PS of Isolated NSs II: from radio-pulsars to magnetars

Population Synthesis of Isolated Neutron Stars with magneto-rotational evolution II: from radio-pulsars to magnetars


Population synthesis studies constitute a powerful method to reconstruct the birth distribution of periods and magnetic fields of the pulsar population. When this method is applied to populations in different wavelengths, it can break the degeneracy in the inferred properties of initial distributions that arises from single-band studies. In this context, we extend previous works to include -ray thermal emitting pulsars within the same evolutionary model as radio-pulsars. We find that the cumulative distribution of the number of X-ray pulsars can be well reproduced by several models that, simultaneously, reproduce the characteristics of the radio-pulsar distribution. However, even considering the most favourable magneto-thermal evolution models with fast field decay, log-normal distributions of the initial magnetic field over-predict the number of visible sources with periods longer than s. We then show that the problem can be solved with different distributions of magnetic field, such as a truncated log-normal distribution, or a binormal distribution with two distinct populations. We use the observational lack of isolated NSs with spin periods s to establish an upper limit to the fraction of magnetars born with G (less than 1%). As future detections keep increasing the magnetar and high-B pulsar statistics, our approach can be used to establish a severe constraint on the maximum magnetic field at birth of NSs.

stars: neutron – pulsars: general – stars: magnetic fields.

1 Introduction

The last couple of decades have seen tremendous advances in our understanding of the properties and evolutionary paths of neutron stars (NSs), largely driven by multi-wavelength detection. Despite this progress, however, several puzzles still remain. In particular, the apparent diversity of the observational manifestations of NSs, which has led over the years to several subclassifications (i.e. rotation-powered pulsars, dim isolated NSs, Anomalous X-ray Pulsars, Soft-Gamma Ray Repeaters, Central Compact Objects, among others) has highlighted the need of understanding possible evolutionary links among different classes, as well as the role played by the birth properties of the NSs.

Studying the properties of NSs as a whole is best done via Monte Carlo simulations aimed at reproducing the galactic population. Since the fraction of pulsars that can be detected close to their birth constitutes a negligible fraction of the total sample, population synthesis studies (PS in the following) generally use the present-day observed properties of pulsars, together with some assumptions about their time evolution, to reconstruct the birth distribution of periods and magnetic fields for the pulsar population.

PS modeling of NSs has a long history, and especially so in the radio band, where the largest sample of pulsars has been detected (e.g. Gunn & Ostriker 1970; Phinney & Blandford 1981; Lyne, Manchester & Taylor 1985; Stollman 1987; Emmering & Chevalier 1989; Narayan & Ostriker 1990; Lorimer et al. 1993; Hartman et al. 1997; Cordes & Chernoff 1998; Arzoumanian, Chernoff & Cordes 2002; Vranesevic et al. 2004; Faucher-Giguère & Kaspi 2006; Ferrario & Wickramasinghe 2006) and in some cases in other wavelengths (Gonthier, Van Guilder & Harding, 2004; Pierbattista et al., 2012). The efforts put over the years into this area of research stem from the fact that these studies allow to constrain the birth properties of NSs, which are in turn intimately related to the physical processes occurring during the supernova (SN) explosion and in the proto-NS. As such, they bear crucial information on the physics of core-collapse SNe, in which most NSs are believed to be formed.

The energy reservoir powering Rotation Powered Pulsars (RPPs) is provided by the loss of rotational energy. RPPs are mostly seen in the radio band ( objects), and we can detect more than one hundred of them also in the -ray and/or -ray frequency range. Part of the rotational energy loss is converted into non-thermal synchrotron and curvature radiation by acceleration of charged particles, which are either extracted from the surface or created by pair production. According to the classical scenario, such acceleration takes place in some regions, called gaps, where a cascade of pairs is supported by the same high-energy photons emitted by the particles, and is thought to produce the coherent radio emission we observe. While radio emission is likely produced above the polar cap (Sturrock, 1971; Ruderman & Sutherland, 1975), high-energy photons are thought to originate in the outer magnetosphere (Cheng, Ho & Ruderman, 1986; Romani, 1996; Muslimov & Harding, 2003) or in the wind region (Pétri, 2012). Of particular interest is the recent FIDO model (Kalapotharakos, Harding & Kazanas, 2014), according to which the gamma-ray emission is produced in regions near the equatorial current sheet. These models can fit the observed values of the radio-lags and reproduce the observed light curve phenomenology and the GeV photon cut-off energies.

The majority of PS studies in the literature have assumed that the magnetic field is fixed at its value at birth, or have adopted simple, parametrized expressions for its decay. Recently, a large theoretical effort has been devoted to model the magnetic field evolution with magnetothermal simulations which follow the coupled evolution of the temperature and the magnetic field in the NS interior (Viganò et al., 2013). We will use results from these simulations to model the magnetic field evolution.

Gullón et al. (2014) (paper I, hereafter) have revisited PS modeling of isolated radio-pulsars incorporating the -field evolution from the state-of-the-art magneto-thermal evolution models, together with recent results on the evolution of the angle between the magnetic and rotational axes from magnetospheric simulations. In an era where the availability of high energy satellites has enlarged considerably our understanding of the NS population, unvealing a variety of observational manifestations of NSs, it is then timely to include these further constraints in the global PS study of the NS population. A first attempt in this direction, including part of the X-ray emitting NS population and magnetic field decay models, was performed by Popov et al. (2010), with the specific goal of uncovering links among the apparently different NS manifestations.

In line with recent theoretical efforts aimed at understanding the various observational manifestations of NSs, with a unifying evolution model (Perna & Pons, 2011; Pons & Perna, 2011; Viganò et al., 2013), our goal in this work is to combine previous studies of the radio-pulsar population (paper I) with the study of the population of -ray pulsars. For young and middle-age NSs ( Myr), the thermal emission comes from the internal heat which is gradually released or, in the case of NSs with large magnetic fields ( G), by Joule dissipation of the electrical currents circulating in the NS crust which enhances the thermal luminosity. We do not include in this study non-thermal emission powered by losses of rotational energy, mainly because of the lack of detailed models for the physical mechanisms. Predictions for the non-thermal luminosities are based on phenomenological correlations with the rotational energy losses, which include more free parameters that are not linked to physical properties. On the contrary, independent, physically motivated cooling models provide the information linking the model parameters with the thermal emission of NSs. We hence focus on this population of NSs and leave the study of rotation powered, non-thermal emission (-ray and -ray pulsars) for future works.

The article is organized as follows: In § 2 we summarise the main results of paper I about the radio-pulsar population and we describe the assumptions made to model the X-ray thermal radiation. In § 3 we discuss the results of our analysis and finally, in § 4 we present our conclusions and final remarks.

Model Decay Envelope
[G] [G] [s] [s] [century]
A No decay Heavy elements
B Slow Heavy elements
C Fast Heavy elements
D Slow Light elements
E Medium Light elements
F Slow (with toroidal field) Light elements

Table 1: Optimal parameters for the radio-pulsar population for each magnetothermal evolution model considered in this work. Columns to indicate the values of the free parameters that minimize the -value computed through a 2D KS test. Column shows the corresponding birthrate, and column 10 shows the average -value of 10 realizations and its standard deviation .
Figure 1: Contour plot of the -value computed through a 2D KS-test in the - plane (a 2D cut of the general 5D parameter space), for models B (top), C (middle) and D (bottom). The purple crosses indicate the minima found by the simulated annealing procedure in each case, and listed in Table 1.

2 Multi-wavelength population synthesis

From an observational point of view, there is an important qualitative difference, regarding detectability biases, between the radio band and the high-energy bands. In the former, relatively complete surveys allow to quantitatively implement the instrumental detectability in the PS study. In -rays, instead, the sample is not complete, and the sky coverage is inhomogeneous. Moreover, the detectability depends on non-trivial issues, such as the presence of a counterpart with known ephemerides in other wavelengths, or whether or not it has been associated to a transient -ray outburst. This implies that the potential number of visible NSs above a certain threshold flux should be larger than the actual number of detected sources.

In -rays we have to account for the effect of absorption by the interstellar medium. The hydrogen column density, , for each synthetic pulsar is:


where is the local hydrogen number density, whose spatial distribution is estimated following Zane et al. (1995), and the integral is performed along the line of sight . The energy-dependent absorption is calculated following the models of Balucinska-Church & McCammon (1992), which include photoelectric cross sections of elements in the keV range. This allows us to calculate the detected (absorbed) flux on Earth , assuming that radiation is isotropically distributed.

2.1 The radio pulsar population

For completeness, we begin our discussion by reviewing the main results of paper I, to which we refer for technical details. An outline of the procedure to generate a NS population is sketched in the following way:

  • A large number of NSs with uniformly distributed random ages are generated in the framework of a given magneto-thermal evolution model (we explore six different models), which determines the evolution of the magnetic field. The birth location distribution, kick velocity, and evolution in the galactic gravitational potential are described in paper I.

  • Gaussian distributions are assumed for both the initial spin period and the logarithm of the initial magnetic field, . The central values (, ) and widths (, ) are left as free parameters. The radio-luminosity (at 1400 MHz) obeys the following equation:


    where is the spin period, is the period time derivative, mJy kpc and is a random correction chosen from a zero-centered gaussian of (as in Faucher-Giguère & Kaspi 2006), to take into account the large observed dispersion. The index is a free parameter.

  • The NSs are classified as radio-visible according to the same flux and beam visibility criteria of the radio survey used in paper I. The generation of pulsars stops when the number of radio-visible pulsars is equal to the number of radio pulsars in the considered observational sample. This fixes our normalisation for the total number of sources or, equivalently, the birth rate for each synthetic population.

  • The set of free parameters (, , , and ) is determined for each model by minimizing the associated -value of a generalized two-dimensional Kolmogorov-Smirnov test. The statistic is calculated from the four natural quadrants that define each point in the plane of the observed sample. The fraction of data from both samples (observed and simulated) in the radio band, is calculated on each quadrant, and the -value is defined as the maximum difference between observed and simulated fractions over all points and quadrants (for a more technical description we refer the reader to Fasano & Franceschini 1987; Press et al. 1993).

The results for the 6 evolutionary models studied are summarized in Table 1. Here we considered three models with iron envelopes (A, B and C, already presented in paper I), differing in the timescale of the magnetic field decay. We have also used three new models (D, E and F) with light element envelopes. The last one also has an additional strong toroidal magnetic field. Although the envelope composition barely affects the evolution of rotational properties ( and ), it has a relevant impact on the thermal luminosity (see next section).

In general, we find that the differences in the initial period distribution between different models are not significant. Indeed, any broad distribution of initial periods in the range s can be reconciled with the data, but narrow distributions peaked at short spin periods ( ms) are ruled out. The value is more dependent on the initial magnetic field distribution, but this is not strongly constrained because its central value and width are correlated showing an interval of possible solutions: the larger the central value, the wider the distribution. This holds true for models with and without field decay, and with different assumptions for the magnetospheric torque. This correlation is visible in Fig. 1, where we show contour plots of the -value in the plane. We note this is a 2D plot of a section of a 5D parameter space. For conciseness, we only show plots for models B,C, and D (other models look very similar).

The dark blue regions, where the -value is low, define the range of allowed parameters. The purple crosses mark the different solutions for which explicit values for the parameters have been given in Table 1. The main conclusion from the analysis of only the radio pulsar population is that, independently of the underlying physical model, one can find acceptable parameter sets to fit the radio population reasonably well. Thus the radio-pulsar population alone cannot help much to discriminate among different evolutionary models. The next question is whether or not the analysis of pulsar distributions in other bands can break the degeneracy.

Figure 2: Left: Blackbody spectrum with 1 keV (dashed line), and reprocessed spectra using the RCS model with and (solid color lines). Unabsorbed fluxes with a typical galactic value of cm are plotted. Middle: Ratio between the absorbed flux of pure blackbody model to that of the RCS model as a function of the blackbody flux . Right: same as the middle panel, as a function of the logarithm of the magnetic field. The results show only the NSs with detectable thermal emission in the 0.1-10 keV band for model E, with fluxes above erg s cm.

2.2 Thermally emitting NSs

The observational sample

We consider all NSs for which the thermal emission is dominant, or clearly contributing to the total power in the soft X-ray band, and whose origin is attributed to residual cooling from birth. The selected sample we use has been described in Viganò et al. (2013) and it is periodically updated in our website3, includes thirteen closeby RPPs (only those with a clear thermal component), seven X-ray dim Isolated NSs (XINSs), and seventeen magnetars.

The XINSs are a small group of nearby (within about 200-300 parsecs), thermally emitting, isolated NSs. They are radio quiet but relatively bright in the X-ray band, and five of them have periods similar to magnetars. Moreover, their inferred magnetic fields are between the upper end of the radio pulsar population and the magnetar candidates. At variance with the majority of X-ray emitting pulsars, XINSs spectra are rather perfect blackbodies (kT0.1 keV, occasionally with some broad lines), which make them the key objects for research on the NS equation of state, and the cooling properties of the NS crust and core.

The “magnetars” (Olausen & Kaspi, 2014) are a small group of X-ray pulsars with spin periods between s, whose properties cannot be explained within the common scenario for pulsars. Their inferred magnetic fields are  G. The decay and instabilities of such large magnetic fields are thought to power their emission (Duncan & Thompson, 1992; Thompson & Duncan, 1993). Their powerful persistent X-ray emission (much larger than the rotational energy loss) is well modelled by a blackbody with  keV with an additional power-law component (), usually interpreted as the product of Resonant Compton Scattering (RCS) of the seed thermal photons in the dense magnetosphere (Thompson, Lyutikov & Kulkarni, 2002). Hence even this apparent non-thermal component has a thermal origin. Different physical models have succeeded to reproduce the effect of the RCS to fit the magnetars X-ray spectra (Rea et al., 2008; Zane et al., 2009; Beloborodov, 2013). Given the transient nature of these objects (they show flares and sometimes long radiative X-ray outbursts), we have used in this work only X-ray luminosities derived during their quiescent state (available for the 17 selected objects).

Magneto-thermal evolution and emission model

In NSs endowed with strong magnetic fields, the evolution of the internal and surface temperatures and the magnetic field are closely linked. The thermal luminosity is strongly dependent on the time-dependent magnetic field strength and topology, while the evolution of the magnetic field depends on the local values of the magnetic diffusivity and the electron relaxation time, which strongly depend on the temperature. Viganò et al. (2013) performed 2D magneto-thermal simulations of NSs, and studied the evolution of magnetic field strength and thermal luminosity as a function of time. Using this code, we produce a number of tables with the effective temperature and value of the dipolar magnetic field as a function of the age of the NS. Among the many model parameters (initial magnetic field topology, mass, radius, envelope composition, etc), we will focus on the parameters that more strongly influence the observable quantities (luminosity, , ) at a given age. These are the envelope composition (light vs. heavy elements), the initial magnetic field strength at the pole, , and the impurity parameter in the inner-crust , a measure of the average charge distribution of the nuclei present in the crust. We refer to Viganò et al. (2013) for more details.

The magnetic field strength determines the rotational evolution of the NS. In addition, in highly magnetised NSs, the additional energy reservoir provided by the magnetic field dissipation may further increase the surface temperature, and extend the duration of the stage in which neutron stars are bright enough to be seen as X-ray pulsars. was found to be a crucial parameter to understand the clustering of spin periods in isolated X-ray pulsars and why there is an upper limit of isolated X-ray pulsars at about 12 s (Pons, Viganò & Rea 2013). The composition of the envelope is also of particular relevance. Generally speaking, low-B NSs with light element envelopes have significantly larger temperatures and luminosities during the first yr of their life, and lower temperature/luminosity afterwards, once they enter the photon cooling era (Page et al., 2004; Yakovlev & Pethick, 2004; Potekhin, Chabrier & Yakovlev, 2007; Page et al., 2011).

These three parameters characterize the evolution of the magnetic field strength and temperature as the star ages, and therefore the trajectories that NSs follow in the - diagram. Having all these tabulated quantities computed via simulations, for each star generated in the synthetic sample we calculate the present value of magnetic field, , , temperature and luminosity. We assume blackbody emission, and integrate over the whole star surface to calculate the luminosity. For NSs with high magnetic fields, we also add the effect of resonant Compton scattering (see next subsection) in the spectrum. Among all NSs generated (not only those seen as radio-pulsars), we select those for which the absorbed X-ray flux at the Earth location is above erg s cm, as potentially visible X-ray pulsars (for present or future instruments).

Resonant Compton Scattering

The thermal emission is usually assumed to be blackbody-like. However, it is well known that magnetar spectra show a non-thermal tail, whose origin cannot be attributed to the rotational energy of the star. The most popular interpretation is that magnetospheric plasma boosts up thermal photons coming from the surface. For these reasons, it is important to account for the possible effects that spectral distortions in the form of a hard tail would have on the detectability of NSs in X-rays.

We adopt a simplified model for the resonant Compton scattering (RCS, Lyutikov & Gavriil 2006) of seed photons (assumed to have a pure blackbody spectral energy distribution) in NS magnetospheres, and successfully applied to model spectra of several magnetars (Rea et al., 2008). Although it is a 1D, plane-parallel model, far from describing the non-trivial interaction between plasma and radiation in a realistic 3D magnetosphere, it effectively fits X-ray data. RCS models are parametrized by two parameters: the resonant optical depth along a ray , which is related to the plasma density, and the thermal velocity of electrons in units of the speed of light, (see Lyutikov & Gavriil 2006 for more details).

In Fig. 2 (left) we show the comparison between a blackbody of keV, and the spectra predicted by the RCS model for and different values of . We show the spectra as they would be seen when absorbed by the interstellar medium with cm, a typical value for magnetars in the galaxy. Since the interstellar absorption is important below keV, an efficient RCS could enhance the NSs detectability, because it promotes low energy photons to higher energies.

Rea et al. (2008) fit the X-ray spectra in quiescence of several magnetars. We found a correlation between and the RCS parameters and , which can be reproduced by the following analytical fits:


We implement this analytical fit to modify the spectra of the randomly generated NSs in the population synthesis code. RCS does not change the original blackbody spectrum in NSs with low magnetic field, but it generates a hard tail in magnetars. In order to show the effect of the RCS in a typical PS realisation, in the middle panel of Fig. LABEL:fig:rcs2 we show the ratio of the RCS flux to the pure blackbody flux in the keV band, as a function of the blackbody flux. This figure clearly shows the main effect of the RCS correction: a flux enhancement by a factor of a few (up to a factor of 20) is applied to a significant number of sources, which has a visible imprint in the statistics of the number of detectable NSs. This effect is more evident for the dimmest sources: since RCS affects the high-energy tail, we can only obtain a very large enhancement for sources with a large , and these are always dim. The right panel of Fig. 2 shows the same ratio as a function of the initial field strength. It shows that only for the highest fields, in the magnetar range, a significant flux enhancement in the [0.1-10] keV band is produced.

Figure 3: Log N - Log diagrams (top) and period distributions (bottom) for the thermally emitting X-ray pulsars for models A to F and the parameters listed in Table 1. In the bottom panel we show the histogram of sources with erg s cm. In colour lines we show the population synthesis results while the black line shows the observed sources.

3 Results

We begin by analyzing the solutions that best fit the radio-pulsar distribution (Table 1). In Fig. 3 we compare the accumulated number of detected X-ray sources with flux (hereafter - diagram) for each magnetothermal evolution model. The sample of sources is not complete, since a full all-sky survey is not available for the present instruments, and many sources, in particular magnetars, have been discovered while they are in outburst. We can assume that we are actually observing all of the few brightest sources, and we gradually miss more and more sources as the threshold flux decreases. From a quick inspection of the - diagrams, we can immediately conclude that none of the models of Table 1 predicts a number of thermally emitting X-ray sources compatible with the observations. However, as shown in Fig. 1, there is a quite large region in the parameter space with similar statistical significance for the -value. We can then explore whether other parameter combinations (always within the region allowed by the radio-pulsar population analysis) can improve the - results. Since we obtain a systematic lack of bright sources for all models, we have to modify the initial distribution to include more NSs born with high magnetic fields, which can be done if we move up-right in the dark blue region of Fig. 1.

We have proceeded as follows: first, we search for a pair of values of , in the dark blue region of Fig. 1 than can reproduce the observed distribution of X-ray pulsars for fluxes above erg s cm (the other 3 parameters do not affect the X-ray population). Once this is done, we fix the values of these two parameters and repeat the minimization procedure to find new values of , , and that best fit the radio-pulsar distribution. In Table 2 and Fig. 4 (top) we show the results for the parameters that best reproduce the - diagram for all models except model A, for which field decay was not included, and therefore the luminosity is not affected by the magnetic field. In other words, the standard cooling of a neutron star without including any additional heating by magnetic field decay cannot be reconciled with the X-ray data. For the other models, we find acceptable solutions with similar statistical significance for the radio-pulsar fits. The -value obtained with the 2D KS test shown in the last column is an average of 10 realizations, and we can see that they are compatible with the results of Table 1 within 1- or 2-.

Figure 4: Same as Fig. 3 for the parameters listed in Table 2.
[G] [G] [s] [s] [century]
Table 2: Optimal parameters for each magneto-thermal evolution model (defined in Table 1) that fit the observed - distribution. The parameter (see eq. 5) is also shown.

Below erg s cm, the observed number of sources saturates due to the sensitivity limit of the present instruments, combined with the interstellar medium absorption, which only allows us to observe dim sources if they are nearby. Thus, we conclude that the degeneracy shown in Fig. 1, where a large region of the parameter space was allowed by the radio-pulsar data, is severely reduced when the thermally emitting X-ray pulsars are included in the analysis. Interestingly, for different models we obtain similar solutions for the initial magnetic field distribution, in the narrow region , and .

However, to reproduce the - diagram is a necessary, but not sufficient condition to test a model. We now turn to explore in more detail the particular form of the period distributions. We remind again that the comparison is made with an incomplete observational data set, because the number of observed sources is only a fraction of the potentially detectable number of NSs. In Fig. 4 (bottom) we show spin period histograms comparing the observed sample (black solid line) with synthetic models.

An important result can already be seen in these distributions. On one hand, one needs a significant fraction of magnetars/high-B field NSs in the initial distribution to account for the number of X-ray sources already present in our observed sample, but on the other hand, if that fraction is too large, the model over-predicts the number of X-ray pulsars with long periods, even assuming the existence of a highly resistive layer in the inner crust of NSs (Pons, Viganò & Rea, 2013). We can adjust the maximum period reached by a NS by varying the impurity parameter in the innermost part of the crust. Increasing the value of results in the decrease of the maximum period, but still not enough to successfully reproduce the observed upper limit of 12 s and, at the same time, explain the observed number of very bright sources. The strong constraint resulting from combining the need of a few very bright sources with the lack of isolated X-ray pulsars with periods longer than 12 s turns out to be a very restrictive condition. The next relevant question is whether or not unknown observational selection effects are causing this discrepancy.

To explore how the period distribution varies with flux, let us focus on model D (the other models are qualitatively similar). In Fig. 5 we compare the observed distribution to the synthetic data with a different cut-off in the absorbed flux (, , and erg s cm). An interesting trend is the apparent bimodal distribution of periods, similar to the observed distribution, although with a different scale and position. There is a partial lack of NSs with s, as seen in the observations, despite the initial period distribution being a single Gaussian. However, even with the highest flux cut-off ( erg s cm) there are always too many visible X-ray pulsars with long periods. Thus, this discrepancy cannot be attributed only to selection effects.

Figure 5: Same as Fig. 4 (bottom) for model D considering different flux thresholds: (blue), (yellow) and (red). All fluxes are in units of erg s cm.

Since it is impossible to a priori take into account (theoretically) unknown selection effects in -rays, we propose in the following a phenomenological approach. A close inspection of Fig. 4 shows that the observational sample appears to be complete (given all the assumptions made in paper and within statistical fluctuations) for fluxes erg s cm. We have also empirically found that the ratio of the number of observed sources to the number of potentially visible sources in the synthetic sample scales linearly with the flux in the range erg s cm. If we assume that the probability to detect an X-ray pulsar with a given flux is simply


where is the absorbed flux in units of erg s cm, we can obtain for each model the value of that reproduces the diagram. It turns out that for all models (eight column in Table 2). This naive approach is useful to estimate the number of missed sources by uncontrolled selection effects.

We can now generate a random observed sample from the synthetically generated sample by imposing that the probability to detect the thermal emission of a pulsar with a given flux is given by the above expression. The corresponding - diagrams and period distributions after such selection is applied are shown in Fig. 6. Although most pulsars with high periods are now not observed, we still find a non-negligible number of bright enough sources with long periods. Our predictions range from 4 to 20 X-ray pulsars to be observed with s, depending on the model. Even considering that we are dealing with low-number statistics, this appears to be a real problem.

Figure 6: Same as Fig. 4, but applying the phenomenological detection probability formula (eq. 5).
[G] [G] [s] [s] [century]

Table 3: Optimal parameters for models D, E and F with a truncated distribution of initial magnetic fields (G).

3.1 Beyond Gaussian distributions of initial magnetic fields.

All previous results (and most studies in the literature) rely on the assumption that the initial magnetic field distribution is log-normal. As the initial magnetic field increases, old NSs reach longer periods and therefore the narrower beam size produces a strong selection of radio-pulsars. Therefore, the low field part of the distribution is what we can practically constrain using the radio-pulsar population, while the high field part of the distribution remains unconstrained, because the distribution of high field radio-pulsars is highly degenerate with other parameters (beaming angle, radio-luminosity correlation, etc.). We have seen that the optimal parameters for radio-pulsars underpredicts the number of X-ray sources with high fluxes. We hence need to increase the width of the Gaussian to generate more bright NSs with high magnetic fields, but the problem is that the Gaussian becomes too wide, and we run into an overprediction of X-ray pulsars with periods s, which has not been observed. The easiest, and perhaps more obvious, way out is to relax the assumption of an initial log-normal magnetic field distribution. For example, if we consider model D, we can visualize the distribution functions for the best fit to the radio-pulsar population (green/solid) and X-ray populations (brown/dot-dashed) in Fig. 7. The excess of long period X-ray pulsars is originated by the high tail of the latter. We now discuss two possibilities, to reconcile the synthetic models with the data.

Figure 7: distributions for model D that best fit the radio-pulsar (green/solid line) and X-ray pulsar (brown/dot-dashed) populations. In red/three-dot dashed and blue/dashed we plot the results for truncated and bimodal distributions discussed in Section 3.1, respectively.
Figure 8: Same as Fig. 6 assuming a log-normal distribution with cutoff at G, and the parameters of Table 3.

A truncated magnetic field distribution.

If the origin of the neutron star magnetic field is attributed to magneto-hydrodynamic instabilities (e.g., any dynamo process related to rotation or convection), a nonlinear saturation is expected when the amplified magnetic field becomes dynamically relevant to suppress the instability. This may happen at typical fields of the order of G, depending on the particular mechanism (Obergaulinger et al., 2009; Obergaulinger, Janka & Aloy, 2014). Thus, we can simply introduce a cut-off value to the log-normal distribution (red/three-dot dashed curve in Fig. 7). In Fig. 8, we can see that the - diagram for the X-ray population can also be successfully reproduced by a truncated distribution with G and increased and (Table 3). The period distribution is significantly modified, and now there are no visible sources with s, as shown in Fig. 8.

(%) (%) [century]

Table 4: Results for models D, E and F with a bimodal distribution. We show the relative weight of the first population (, same parameters as in Table 1) and second population of highly-magnetized stars (, flat distribution in the range ).

A second population of magnetised NSs.

The alternative possibility that we have explored is the existence of a second population of magnetised NSs, whose origin is different from that of standard radio-pulsars. The reason for this second population could be attributed to two different evolutionary paths of massive stars if they are isolated or in binary systems. Since about 80% of massive stars are in binaries, and assuming that half of them go through a common envelope phase or some evolutionary stage that produces a magnetized stellar core, one could expect about 40% of the total pulsar population being born with a binary massive progenitors, possibly leading to the formation of more magnetized NSs (Spruit, 2008; Langer, 2012; Clark et al., 2014). We have investigated what happens if, in addition to the population generated with the parameters of Table 3 by fitting only the radio pulsar population (hereafter denoted by NS1), we add a second population (denoted by NS2), uniformly distributed in the range and normalized to account for the total number of NSs observed (blue/dashed curve in Fig. 7). The results are shown in Fig. 9, where again we can observe that both the - diagram and the period distribution are satisfactorily reproduced. We have verified that this second component does not modify the radio pulsar population, since their magnetic fields are so large that the vast majority of NS2 sources remains undetectable as radio-pulsars. The total birthrate, however, will increase accordingly, still marginally compatible with the observations. We also note that one can obtain similar results by assuming other distributions for this second population, for example a narrow log-normal distribution centered in G.

Figure 9: Same as Fig. 6 with a bimodal distribution of described in §3.1.2.

3.2 The upper limit to the birth rate of very strongly magnetized NSs.

Before concluding our analysis on the period distribution, a simple calculation can stress once more the conflict between having a significant fraction of NSs with magnetic fields at birth above a certain value, and the lack of observed isolated NSs with periods s. To quantify this incompatibility, we begin by obtaining the probability of observing a star, born with a given initial magnetic field , as an X-ray pulsar with a period longer than 12 s. The results for three different values of the initial magnetic field are given in Table 5, for the representative model D (similar values are obtained for other models). In this calculation, we assumed a uniform distribution of ages in the range 0-560 Myrs (with constant birth rate), and we generate a large number of NSs with the same magnetic field. Then, we can obtain the probability as the fraction of synthetic visible sources with s to the total number of stars generated.

Next we address the following question: what is the maximum number of stars born with compatible (within statistical fluctuations) with the absence of observed neutron stars with periods longer than 12 s? This number can be estimated taking into account that if is the number of stars born with a given magnetic field, the probability of observing none of them with periods longer than 12 s is given by , according to the Poisson distribution, and the probability that we will observe at least one star with  s is .

We can use this value to reject the null hypothesis in the following sense: for a given model, we can calculate the value of that corresponds to a probability of , and we know that in 99% of the realizations with generated stars, we will not detect any source with  s. Therefore, we can conclude that if a model predicts stars, the fact that none has been observed rejects the model at the 99% confidence level (statistical fluctuations may result in a no-detection only in of the cases).

In Fig. 10 we plot the probability of no-detection of NSs with periods longer than 12 s against the birthrate, for different values of . This can be interpreted as an upper limit to the birthrate of stars with initial magnetic field , for a given confidence level. For a confidence level of 99% and a magnetic field of G we obtain . This value corresponds to a birthrate of objects/century which means, by comparison with the total galactic NS birthrate NSs/century, that less than can be born with magnetic fields G. Of course, this fraction is even smaller if we reduce our confidence level. The maximum birthrate for stars with initial magnetic field greater than G can be times larger than for G for the same confidence level. This is satisfied by both alternative models described in the two subsections above, which cannot be discriminated with current data.


Table 5: Probability of detection of stars whose period is longer than s for samples with fixed initial magnetic field .
Figure 10: Probability of no-detection of NSs with s as a function of the birthrate, for three different initial magnetic fields: G (blue squares), G (green triangles) and G (red diamonds). The symbols show results of simulations with 20 million NSs and the dashed lines correspond to the theoretical Poissonian value .

3.3 distribution.

To complete our discussion, we comment on our results for the distribution of . In Fig. 11 we show the histograms for the models considered, which roughly agree with the observations. However, all models and initial field distributions (log-normal, truncated log-normal, bimodal) show a similar trend: synthetic samples display an excess of sources with mean values s s and a lack of a few sources with s s. This is due to the fact that even the strongest magnetic fields ( G) hardly achieve s s.

The observed sources with very high are magnetars and it is likely that modelling their spin evolution with the simple magnetospheric torque formulae (Spitkovsky, 2006; Beskin, Istomin & Philippov, 2013; Philippov, Tchekhovskoy & Li, 2013) is an oversimplification. For instance, strong particle winds can be non-negligible in the early yr of their lives (Tong et al., 2013) when these sources undergo frequent bursts and flares. As a result of additional torques, the can increase by an order of magnitude for the same magnetic field. This would imply that the magnetic field estimate by the classical magnetodipolar formula is a factor of a few higher than the real value. This effect is only expected to work during a short time and does not change the results of the analysis of the whole NS population, but it may contribute to the anomalous high values of of the few youngest magnetars.

Figure 11: Observed and predicted distributions. Top: models with log-normal distributions. Middle: models with truncated distributions discussed in section 3.1.1. Bottom: models with bimodal distributions discussed in Section 3.1.2. In colour lines we show the PS results while the black line shows the observed sources.

4 Summary

We have extended previous works on population synthesis of isolated neutron stars to study the combination of initial conditions and microphysical parameters required to explain the radio-pulsars and thermally emitting -ray pulsar populations within the same evolutionary model.

Our first finding is that the - diagram of X-ray pulsars can be well fitted by several models that, at the same time, reproduce the radio-pulsar distribution (because the parameter space of initial magnetic field distribution of radio-pulsars has a large region with the same statistical significance). Models with iron envelopes require a higher mean value of the initial magnetic field and a wider distribution than those with light elements.

Our second relevant result is that the obtained sample of observed sources with erg s cm is nearly complete, and we estimate that, with more sensitive instruments, the thermal emission of about one hundred of NSs with fluxes erg s cm is potentially detectable. This number can be increased up to 500 sources if the sensitivity is good enough to cover the whole sky with erg s cm.

However, when examining in detail the period distribution all models with log-normal distributions of the initial magnetic field fail to explain why there are no X-ray pulsars with periods longer than 12 seconds. This sharp upper limit cannot be reproduced, even with models with a very high resistivity in the crust/core interface, if the initial field distribution is assumed to be log-normal. The problem is solved if other distributions are invoked, such as a truncated log-normal distribution, or a binormal distribution with two distinct populations. The most interesting feature, common to the two alternative cases considered, is that both cases require that magnetars cannot be born with a very high magnetic field. For a confidence level of 99%, we obtain an upper limit to the birth rate of NSs with G of objects/century, which represents a tiny fraction of the population (less than 1%).

The existence of an upper limit to the magnetic field at birth in the range G, combined with fast magnetic field dissipation and observational selection effects, explains the lack of isolated pulsars with long periods. Although this has a drawback, the lack of a few pulsars with very high values of , note that we have assumed that NSs spindown only by the rotating magnetospheric torque. In very young, active magnetars, losses of angular momentum carried by particles (wind) can compete with the magnetospheric torque (Tong et al., 2013), increasing the value of for a given magnetic field. A more detailed study of these phenomena, and other modifications of the spin-down evolution is reserved to future works.


This work was supported in part by the grants AYA2013-42184-P and Prometeu/2014/69, and by the New Compstar COST action MP1304. MG is supported by the fellowship BES-2011-049123. DV and NR acknowledges support from grants AYA2012-39303 and SGR2014-1073. NR is additionally supported by an NWO Vidi award. We gratefully acknowledge useful discussions and suggestions from Peter Gonthier.


  1. pagerange: Population Synthesis of Isolated Neutron Stars with magneto-rotational evolution II: from radio-pulsars to magnetarsReferences
  2. pubyear:


  1. Arzoumanian Z., Chernoff D. F., Cordes J. M., 2002, ApJ, 568, 289
  2. Balucinska-Church M., McCammon D., 1992, ApJ, 400, 699
  3. Beloborodov A. M., 2013, ApJ, 762, 13
  4. Beskin V. S., Istomin Y. N., Philippov A. A., 2013, Physics Uspekhi, 56, 164
  5. Cheng K. S., Ho C., Ruderman M., 1986, ApJ, 300, 500
  6. Clark J. S., Ritchie B. W., Najarro F., Langer N., Negueruela I., 2014, A&A, 565, A90
  7. Cordes J. M., Chernoff D. F., 1998, ApJ, 505, 315
  8. Duncan R. C., Thompson C., 1992, ApJL, 392, L9
  9. Emmering R. T., Chevalier R. A., 1989, ApJ, 345, 931
  10. Fasano G., Franceschini A., 1987, MNRAS, 225, 155
  11. Faucher-Giguère C.-A., Kaspi V. M., 2006, ApJ, 643, 332
  12. Ferrario L., Wickramasinghe D., 2006, MNRAS, 367, 1323
  13. Gonthier P. L., Van Guilder R., Harding A. K., 2004, ApJ, 604, 775
  14. Gullón M., Miralles J. A., Viganò D., Pons J. A., 2014, MNRAS, 443, 1891
  15. Gunn J. E., Ostriker J. P., 1970, ApJ, 160, 979
  16. Hartman J. W., Bhattacharya D., Wijers R., Verbunt F., 1997, A&A, 322, 477
  17. Kalapotharakos C., Harding A. K., Kazanas D., 2014, ApJ, 793, 97
  18. Langer N., 2012, Annu. Rev. Astro. Astrophys., 50, 107
  19. Lorimer D. R., Bailes M., Dewey R. J., Harrison P. A., 1993, MNRAS, 263, 403
  20. Lyne A. G., Manchester R. N., Taylor J. H., 1985, MNRAS, 213, 613
  21. Lyutikov M., Gavriil F. P., 2006, MNRAS, 368, 690
  22. Muslimov A. G., Harding A. K., 2003, ApJ, 588, 430
  23. Narayan R., Ostriker J. P., 1990, ApJ, 352, 222
  24. Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2009, A&A, 498, 241
  25. Obergaulinger M., Janka H.-T., Aloy M. A., 2014, MNRAS, 445, 3169
  26. Olausen S. A., Kaspi V. M., 2014, ApJS, 212, 6
  27. Page D., Lattimer J. M., Prakash M., Steiner A. W., 2004, ApJS, 155, 623
  28. Page D., Prakash M., Lattimer J. M., Steiner A. W., 2011, Physical Review Letters, 106, 081101
  29. Perna R., Pons J. A., 2011, ApJL, 727, L51
  30. Pétri J., 2012, MNRAS, 424, 2023
  31. Philippov A., Tchekhovskoy A., Li J. G., 2013, arXiv:astro-ph/1311.1513
  32. Phinney E. S., Blandford R. D., 1981, MNRAS, 194, 137
  33. Pierbattista M., Grenier I. A., Harding A. K., Gonthier P. L., 2012, A&A, 545, A42
  34. Pons J. A., Perna R., 2011, ApJ, 741, 123
  35. Pons J. A., Viganò D., Rea N., 2013, Nature Physics, 9, 431
  36. Popov S. B., Pons J. A., Miralles J. A., Boldin P. A., Posselt B., 2010, MNRAS, 401, 2675
  37. Potekhin A. Y., Chabrier G., Yakovlev D. G., 2007, AP&SS, 308, 353
  38. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1993, Numerical Recipes in FORTRAN; The Art of Scientific Computing, 2nd edn. Cambridge University Press, New York, NY, USA
  39. Rea N., Zane S., Turolla R., Lyutikov M., Götz D., 2008, ApJ, 686, 1245
  40. Romani R. W., 1996, ApJ, 470, 469
  41. Ruderman M. A., Sutherland P. G., 1975, ApJ, 196, 51
  42. Spitkovsky A., 2006, ApJ, 648, L51
  43. Spruit H. C., 2008, in American Institute of Physics Conference Series, Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More, Bassa C., Wang Z., Cumming A., Kaspi V. M., eds., pp. 391–398
  44. Stollman G. M., 1987, A&A, 178, 143
  45. Sturrock P. A., 1971, ApJ, 164, 529
  46. Thompson C., Duncan R. C., 1993, ApJ, 408, 194
  47. Thompson C., Lyutikov M., Kulkarni S. R., 2002, ApJ, 574, 332
  48. Tong H., Xu R. X., Song L. M., Qiao G. J., 2013, ApJ, 768, 144
  49. Viganò D., Rea N., Pons J. A., Perna R., Aguilera D. N., Miralles J. A., 2013, MNRAS, 434, 123
  50. Vranesevic N. et al., 2004, ApJL, 617, L139
  51. Yakovlev D. G., Pethick C. J., 2004, Annu. Rev. Astro. Astrophys., 42, 169
  52. Zane S., Rea N., Turolla R., Nobili L., 2009, MNRAS, 398, 1403
  53. Zane S., Turolla R., Zampieri L., Colpi M., Treves A., 1995, ApJ, 451, 739
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description