Gaussian mixture model and population synthesis of radio pulsars
aa r X i v : . [ a s t r o - ph . H E ] J un Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 10 October 2018 (MN LaTEX style file v2.2)
Gaussian mixture models and the population synthesis of radiopulsars
A.P. Igoshev ⋆ and S.B. Popov † Sobolev Institute of Astronomy, Saint Petersburg State University, Universitetsky prospekt 28, 198504, Saint Petersburg, Russia Sternberg Astronomical Institute, Lomonosov Moscow State University, Universitetsky prospekt 13, 119991, Moscow, Russia
ABSTRACT
Recently, Lee et al. used Gaussian mixture models (GMM) to study the radio pulsar popula-tion. In the distribution of normal pulsars in the P – ˙ P plane, they found four clusters. Wedevelop this approach further and apply it to different synthetic pulsar populations in orderto determine whether the method can effectively select groups of sources that are physicallydifferent. We check several combinations of initial conditions as well as the models of pulsarevolution and the selection effects. We find that in the case of rapidly evolving objects, theGMM is oversensitive to parameter variations and does not produce stable results. We con-clude that the method does not help much to identify the sub-populations with different initialparameters or/and evolutionary paths. For the same reason, the GMM does not discriminateeffectively between theoretical population synthesis models of normal radio pulsars. Key words: methods: statistical, pulsars: general, X-rays: binaries, gamma-ray burst: general
While radio pulsars are known for more than 45 years(Hewish et al. 1968), many aspects of their initial properties andthe evolution are not well understood. The evolution of a pul-sar is a long process and cannot be studied directly. Fortunately,about 2000 pulsars at different stages of their evolution are cur-rently being observed (Manchester et al. 2005). Therefore, sincethe early seventies, different methods of statistical analysis of theradio pulsar population and its evolution have been proposed (seeVivekanand & Narayan 1981; Vranesevi´c & Melrose 2011 and ref-erences therein). The latest attempt was made by Lee et al. (2012)(hereafter L12).These authors use a Gaussian mixture model (GMM) to iden-tify groups of pulsars in the P – ˙ P plane. In this approach, itis assumed that the density of objects inside a group (a cluster)can be fitted with a two-dimensional normal distribution (in linearor logarithmic scale). Here, we develop further the approach usedby L12 . We analyse the distributions in the space of periods andperiod derivatives, both on logarithmic scale. In other words, thesources belonging to every cluster have log-gaussian distribution.In the plots below, these groups of sources are represented by - σ contours, and we call them ellipses, which actually stands for “agroup selected by the GMM with log-gaussian distribution of den-sity of objects in both coordinates”.We are interested in understanding the origin of the clusters ⋆ E-mail: [email protected] † E-mail: [email protected] We thank Kejia Lee for providing the original code used in the study L12 found by L12 for normal radio pulsars, as well as in checking ifthe GMM can be used to differentiate between different populationsynthesis calculations. We present the analysis of several toy mod-els with well-understood properties, provide an additional analysisof the observed population of the radio pulsars, and finally, applythe cluster analysis to study the realistic synthetic sets of the dataproduced by the advanced population synthesis models.The paper is organized as follows. In the next section webriefly describe the GMM method. In Section 3, the GMM methodis applied to different pulsar populations — both, synthesized andobserved. Our population synthesis models and calculations arealso presented in the third Section. In Section 4, we discuss theeffectiveness of the GMM method to falsify theoretical populationsynthesis models and present two examples of the application ofthis method to other sets of sources. Finally, we present our con-clusions.
The GMM is an iterative method which allows one to find clus-ters of objects that follow the Gaussian distribution. Each set of theclusters describes the data with some probability. We are interestedin finding the one which gives a higher probability for a smallernumber of clusters.The GMM as an unsupervised learning algorithm has demon-strated its efficiency in different areas where it used to ex-tract empirical knowledge from data samples (Kim & Kang 2007;Reynolds et al. 2000). Note that Gaussians on which the method is c (cid:13) A. Igoshev and S. Popov based are not basis functions. Nevertheless, the central limit theo-rem may be used to explain our motivation to use this method. In-deed, a position of a pulsar in the P – ˙ P plane is a two-dimensionalindependent random variable because it is determined by the effec-tive magnetic field and the age. The logarithm of the effective mag-netic field is in the range [10 , , and the logarithm of the age isin the range [3 , . This defines a rather narrow variance range.The central limit theorem states that a set of the large number ofindependent random variables where every variable belongs to adistribution with similar mean and variance has limiting distribu-tion function approaching the normal distribution. This is our mainmotivation to use Gaussians to study the pulsar distribution.All calculations start with only one cluster. We then probe se-quentially ∼
20 different initial positions of a cluster which coverthe entire data range. It is found that the size of the initial clus-ter does not change the results, while the position of the initialcluster is crucial. At the next step, we maximize the likelihoodby varying the model parameters. The details of the procedure oflikelihood maximization, which is called below the Expectation-Maximization algorithm (Press et al. 2007) are described in detailin L12. This procedure is run for all initial guesses for one cluster.As a result, we have a few (typically 4-5 out of 20) sets of clus-ter parameters with significantly different properties. For each ofthem we generate 10 samples of data points, and for each such sam-ple we carry out the two-dimensional Kolmogorov-Smirnov (K-S)test. The most probable set of cluster parameters is used later asa fixed initial guess for the first cluster. Now, we add the secondcluster. For the second cluster, we also make 20 initial guesses. TheExpectation-Maximization algorithm is run, and finally, after ap-plication of the K-S test, we fix the initial guess for both clusters(note, that the initial guess for the first cluster could change whenwe analyse a two-cluster model). Following this algorithm we findall sets of cluster parameters sequentially from one cluster to somemaximum number of clusters.Up to this moment we do not exclude any significantly distinctsets of clusters. When all final sets of clusters (in our study up to4-5 clusters in each set) are found, we again generate 20 samplesof data points.In our numerical experiments, we found that if we use setswith more than 5 clusters then we obtain a large number of differentsolutions with equally high probabilities (this is called overfittingin L12). On the other hand, when we limit the number of clustersto 4-5, we typically obtain no more than one or two solutions withcomparably high probabilities.We compare these samples to each other by means of the K-Stest, and then study the distribution of D ’s (this is the Kolmogorov-Smirnov statistic).Finally, we calculate the probability that the data are drawnfrom the synthetic sample which corresponds to the chosen set ofclusters, and the analysed data (observed or obtained by a popula-tion synthesis model) are drawn from the same parent distributionbased on the mean value of D for the analyzed data.As the final result we select the set of clusters which has thehighest probability (higher than ∼
90 per cent) and which containsthe smallest number of clusters. As it is mentioned in L12, the cu-mulative distribution for the two-dimensional K-S test is not welldefined. Due to this, it was suggested to compute a cumulative dis-tribution of D for every model, and then we rely on this numerical Here it is reasonable to consider logarithms of quantities because we ap-ply the GMM in the logarithmic scale. statistics when we calculate the probability. Based on our experi-ence such numerical cumulative distribution does not provide bet-ter resolution than about a few percent. We say that the probabilityis 99.9 per cent when D lies to the left from the begining (the firstdata point) of the cumulative distribution, i.e. in the region whereour simulation does not have enough resolution. In this section we present the cluster analysis using the Gaussianmixture model for several synthetic and observed populations ofnormal radio pulsars.First, we study the same population as in L12 to be certain thatwe can reproduce their results. Once we have made sure that thisis indeed the case, we apply the method to other populations: theobserved and the synthetic one.In all cases, for a given set of the data points in the P — ˙ P plane (observed or calculated) several runs of the Expectation-Maximization algorithm are made for different sets of initialguesses (see details in L12). These runs produce different sets ofthe Gaussian ellipses (the number of ellipses can also be different)which describe the data with some probability defined according tothe multi-dimensional Kolmogorov-Smirnov test. Typically, we tryto use as small a number of ellipses as possible to describe the datawith a reasonable probability, and we demonstrate only one plotwith the largest probability for the given number of ellipses.When a distribution of pulsars is analysed, the GMM methodprovides multiple solutions, and to choose the most appropriate oneit is necessary to use some formal criteria. We use the K-S test prob-ability as such formal parameter. Some sets of ellipses which do nothave as high probability as the presented ones, potentially can re-flect physically linked groups of pulsars better. However, since apriori we do not have enough information, it is possible that wecan miss such good solutions. L12 made their analysis for a mixed population where, in addi-tion to normal radio pulsars, sources like magnetars and close-bycooling radio quiet neutron stars (NSs) were also considered. Weexclude those objects which have been initially identified not as ra-dio pulsars, but as sources of a different nature (magnetars, rotatingradio transients, cooling near-by NSs, etc.). Results are shown inFigure 1. We present two most probable realisations based on thesame data set.As in L12, four clusters are necessary to describe the data.First, let us note that both pictures are different from the one inL12. As we removed from the data most of the sources with large P and ˙ P (upper right part of the diagram: all magnetars and close-by cooling NSs) the ellipse corresponding to this part of the dia-gram in L12 is changed. However, this is not the only modification.This is especially obvious in the left panel of Figure 1, which has ahigher probability: in this plot one of the ellipses is shifted towardsthe region with smaller ˙ P with respect to the figure in L12. Onlythe ellipse lying along the death line saves its position. In our opin-ion, this argues against the effectiveness of the method to identifyphysically or evolutionary related groups of normal radio pulsars.We also tried to check what happens if the set of the pulsarsconsidered is slightly modified. We randomly excluded 10 per centof the sources (Figure 2). Again, we see that the results are changed. c (cid:13) , 000–000 MM and population synthesis of radio pulsars -20-18-16-14-12-10-3 -2 -1 0 1 2 l og P ˙ , s e c ond s / s e c ond s log P, seconds -20-18-16-14-12-10-3 -2 -1 0 1 2 l og P ˙ , s e c ond s / s e c ond s log P, seconds Figure 1. P – ˙ P diagram for normal pulsars with ellipses derived with the GMM. No magnetars, rotating radio transients, cooling near-by NSs, etc. areincluded into the data set. The probability that the model and the data are drawn from the same parent distribution is 98 per cent for the model in the left panel,and 95 per cent for the model in the right panel. -20-18-16-14-12-10-3 -2 -1 0 1 2 Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 3.
The Parkes multibeam survey pulsars and the related clustersidentified with the GMM. The number of data points is 905. The probabilitythat the model and the data are drawn from the same parent distribution is94 per cent.
In one (left) plot we obtained close to a carbon-copy of the distri-bution from L12, though magnetars are excluded here. In another,which has a very similar probability, high-B pulsars do not forma separate group (note, that the data are the same, only the initialcombination of the ellipses is different).Finally, we take only 905 pulsars detected in the Parkes multi-beam survey (Manchester et al. 2001) to get rid of some of the se-lection effects (Figure 3). The picture is now simplified and be-comes more stable. One three ellipses are required. One is mostprobably related to the existence of the death line. To see if we canfind interpretations for the others we run several population synthe-sis models with different assumptions.
One of the main tasks of this study is to try to understand the originof the clusters found for the observed populations. To do this in asystematic manner we perform several runs for different syntheticpopulations based on well-understood assumptions with regards totheir evolution, the spatial distribution, etc. We use population syn- thesis models of different complexity to see if an addition of a newoption (for example, the field decay or the Galactic spiral structure)produces a new cluster or even changes the whole picture.Let us briefly describe the basic features of our popula-tion synthesis models. Typically, pulsars are born in four spiralarms (unless the opposite is explicitly stated). The spiral armsare parametrized according to Vall´ee (2005). The initial spin pe-riods are taken from the Gaussian distribution with h P i = 0 . s and σ p = 0 . s. The initial magnetic fields are described bythe Gaussian distribution in log-scale: h log B / [G] i = 12 . and σ B = 0 . (if opposite is not explicitly stated). Pulsars are bornwith constant rate, the oldest have ages of years.The motion of every pulsar in the Galaxy is determinedby its birth place, the kick velocity, and the Galactic gravita-tional potential. In all models we use the potential proposed byKuijken & Gilmore (1989). The kick velocity distribution is cho-sen in the form: p ( v l ) = 12 h v l i exp (cid:20) − | v l |h v l i (cid:21) . (1)Here h v l i = 180 km s − (Faucher-Gigu`ere & Kaspi 2006); | v l | denotes the absolute value of v l . Every component of the ve-locity vector is determined according to the probability definedby Equation (1). We neglect the Shklovskii effect as well as thechanges of the relative position of the Sun and the spiral arms.A pulsar is “detected” if its luminosity exceeds the limit S min for the Parkes multibeam survey (Manchester et al. 2001). The ra-dio luminosity is calculated as (Faucher-Gigu`ere & Kaspi 2006): log L = log( L P ǫ P ˙ P ǫ ˙ P ) + L corr , (2)with L = 0 .
18 mJy kpc , ǫ P = − . , ǫ ˙ P = 0 . . L corr isa normally distributed random function with a zero average and σ L corr = 0 . . P is the spin period, and ˙ P is the period deriva-tive in units − s s − . The beaming fraction is calculated as inFaucher-Gigu`ere & Kaspi (2006). A pulsar is detectable only abovethe death-line (Ruderman & Sutherland 1975; Rawley et al. 1986): BP > .
12 10 G s − , (3)where B is the equatorial magnetic field.In all runs the number of objects used in the analysis is equalto 905, i.e. equal to the number of the data-points from the Parkes c (cid:13) , 000–000 A. Igoshev and S. Popov -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 2.
The same as Figure 1, but 10 per cent of the objects are randomly removed. The probability that the model and the data are drawn from the sameparent distribution is 94 per cent for the model in the left panel, and 95 per cent for the model in the right panel.
Table 1.
List of synthetic models.Name Fig. Electron density Prob. RemarksModel I. 4 NE2001 0.99 –Model II. 5 DM = 15 D DM = 15 D DM = 15 D R Gal
Model VI. 9 NE2001 0.94 bimodal in P Model VII. 10 NE2001 0.97 bimodal in B Model VIII. 11 NE2001 0.98 Popov et al. (2010),no decayModel IX. 11 NE2001 0.96 Popov et al. (2010),decay survey (see Fig. 3). The analysis based on the GMM appears tobe dependent on the number of pulsars considered. Therefore, wedecide to use the same number of pulsars in all generated samplesas it is in the Parkes multibeam survey (see 3.1).All models are listed in the Table 1, and short comments aregiven for each.
We start with the simplest models with the standard magneto-dipolelosses with a constant angle between the spin and the magneticdipole axis. The magnetic field is constant, too.In the first model that we present, the electron density is cal-culated according to the NE2001 (Cordes & Lazio 2002). The datacan be well described with four clusters, see Figure 4. Three clus-ters are rather similar to those we see in the fit for the Parkes multi-beam survey data (Figure 3). The left cluster, containing just fourpulsars, breaks this similarity. Note that this cluster is statisticallysignificant. Indeed, the GMM finds no similar models with justthree clusters. And the most probable model with three clusters isdrawn from the same parent distribution as the data with the prob-ability 92 per cent (compare with 99 per cent for 4 clusters).The results presented in Figure 4 reproduce all main featuresseen in the plot for the Parkes data (Figure 3). However, slightlydifferent runs of the population synthesis code, where the integra-tion of the individual pulsar histories lasted not for yrs, butless, for example yrs, produce very different results in terms -20-18-16-14-12-10-3 -2 -1 0 1 2 Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 4.
Model I. The magnetic field is constant. Pulsars are born in spi-ral arms. The electron density is calculated with the NE2001 model. Theprobability that the model and the data are drawn from the same parentdistribution is 99 per cent. of the identified clusters. In our opinion, this demonstrates that evenfor the simple models with well understood ingredients the methoddoes not produce a stable set of easily interpretable clusters.We can simplify the Model I by using a less complex modelfor the electron density distribution in the Galaxy. In the Model II,the dispersion measure is calculated with a simple relation: DM = 15 D. (4)Here D is the distance in kpc. The minimum brightness isdependent on the interstellar scattering time, τ scat . If the DM istaken in the form (4), then τ scat is calculated according to Equation(7) from Lorimer et al. (2006).If we compare Models I and II, the results are changed sig-nificantly. In the simpler model, only two clusters are necessary todescribe the data, see Figure 5. We conclude that the effects con-nected with the fluctuations of the electron density are strong, andinfluence the number and the distribution of ellipses quite signifi-cantly.In addition, we made calculations for the modified simplestModel IV when the electron density distribution is calculated ac-cording to NE2001. We do not present the figure, but the results c (cid:13) , 000–000 MM and population synthesis of radio pulsars -20-18-16-14-12-10-3 -2 -1 0 1 2 Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 5.
Model II. The same as Model I in Fig.4, but instead of the NE2001model we use a simple linear relation between distance and dispersion mea-sure (Eq. 4). The probability that the model and the data are drawn from thesame parent distribution is 89 per cent. are changed significantly. This again confirms that the influence ofelectron density distribution on the picture of gaussian clusters isessential.In principle, the identification of just two clusters in a set ofthe observational data can lead to a conclusion of some dichotomyeither in the initial properties or in the evolutionary laws. Here,obviously, the pulsars are born from a single mode population. Theevolution also does not contain any process which can result in adichotomy of sources. Therefore, the origin of such dichotomy ispuzzling. We will simplify the model even more in Section 3.2.3in an attempt to clarify it.
Now we add to our scenario the magnetic field decay. The firstmodel of this type (Model III) is similar to the one shown in Fig-ure 4, but the field decays exponentially. For the illustrative pur-poses we choose a simple model: B ( t ) = B exp (cid:20) − tτ mag (cid:21) . (5)Here τ mag = 5 10 years. We see (Figure 6), that the set of ellipsesis very different from the Model I. In most cases the pulsars withlarger ˙ P are joined in one Gaussian with the pulsars from the mainpopulation, so it is difficult to identify if there is a magnetic fielddecay in this model, looking at the distribution of the ellipses. Let us make further simplifications of the Model II presented inFigure 5. Now we exclude the spiral structure. This is the simplestmodel we study here (Model IV), see Fig. 7. The picture looks onlyslightly different.Another modification related to the spatial distribution is veryartificial. We want to study how the finite size of the Galaxy in-fluences the picture. Therefore, we increase the size of the Galaxy, R Gal , by a factor of three, and increase the pulsar production raterespectively. The results for the Model V are presented in Figure 8.Now the distribution is described with three clusters and two ofthem are elongated along the death-line. -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 6.
Model III. Similar to Model I, but the magnetic field decays ex-ponentially. The probability that the model and the data are drawn from thesame parent distribution is 99.9 per cent. -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 7.
Model IV. Similar to Model II, but the birth places of the pulsarsare not related to the spiral arms. The probability that the model and thedata are drawn from the same parent distribution is 96 per cent. -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 8.
Model V. Similar to Model II, but the size of the Galaxy and thetotal birth rate are increased. The probability that the model and the data aredrawn from the same parent distribution is 92 per cent.c (cid:13) , 000–000
A. Igoshev and S. Popov -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 9.
Model VI. The initial spin period distribution consists of twoGaussians: h P i = 0 . s, σ P = 0 . s, and h P i = 0 . s, σ P = 0 . s.The probability that the model and the data are drawn from the same parentdistribution is 94 per cent. In this subsection we present two models with modified initial dis-tributions of the spin periods and the magnetic fields. Instead of us-ing a single Gaussian for each of these parameters, in one case weuse two well-separated Gaussian distributions for the initial spinperiods; in another we use two Gaussians for the magnetic fields(which are assumed to be constant during evolution).In our Model VI the initial spin period distribution consists oftwo Gaussians: h P i = 0 . s, σ P = 0 . s, and h P i = 0 . s, σ P = 0 . s.Three clusters are identified, see Figure 9. They overlap suf-ficiently. We suppose that for the GMM it is difficult to distin-guish pulsars from different sub-populations because of their sig-nificant mixing. After a short time (about – years) pulsarswith shorter periods are braked enough and can be confused withyounger pulsars from the second sub-population with longer initialperiods.The picture of the clusters distribution for the case of the bi-modal distribution of the initial magnetic field could look simpleras pulsars with stronger magnetic fields are higher in the P – ˙ P plane, and they are never mixed with low-magnetized pulsars with-out the field decay. Nevertheless, in the case of two Gaussians forthe initial magnetic field we have to make them very well sepa-rated to identify these groups using the GMM-code. When the firstGaussian is centred on log B / [G] = 12 . and the second — on log B / [G] = 13 . (both with σ B = 0 . ), the GMM is not able todistinguish them. Even for log B / [G] = 12 . and log B / [G] =13 . (both with σ B = 0 . ) the method does not work well. Onlywhen we take log B / [G] = 11 . and log B / [G] = 12 . with σ B = 0 . (Model VII) the two groups are clearly described bydifferent ellipses. However, in this case it was also easily visibleby eye, see Figure 10. Note, that the additional clusters are alsonecessary to describe the data, similar to Model I. -20-18-16-14-12-10-3 -2 -1 0 1 2 Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 10.
Model VII. The initial magnetic field distribution is bimodal: log B / [G] = 11 . , σ B = 0 . and log B / [G] = 12 . , σ B = 0 . .The probability that the model and the data are drawn from the same parentdistribution is 97 per cent. It is interesting to apply the GMM to the results of the advancedmodels of the population synthesis. We use results of the calcula-tions using the code of Popov et al. (2010) .We study two data sets: one with and one without the mag-netic field decay. Both are fitted to reproduce the properties of theobserved population. The results are presented in Figure 11. In thecase of the constant field (Model VIII, left panel), the results canbe described with just two clusters, as in the case of the most sim-ple Model IV (Figure 7), but the ellipses have different properties.If the field is allowed to decay (Model IX), then we need threeclusters. However, the ellipses now are clearly determined by theoutlying points. The method does not distinguish the group of theinitially higher magnetized pulsars which experienced significantfield decay. Note, that in the original calculations no death-line wasused, only the selection related to the flux, so no ellipse stretchedalong the death-line is seen. Note however, if we exclude pointsbehind the death line, then the pictures is not changed significantly. Based on the significant amount of numerical experiments we canconclude that the GMM method is not effective for the studies ofnormal radio pulsars. Let us try to speculate why this is the case.The main initial distributions of pulsars in most of our modelsis defined by two Gaussians: one for the spin periods and the sec-ond — for the magnetic fields (in log-scale). However, the GMMmethod is applied to analyze the distributions in the P – ˙ P not inthe B – P plane. For the standard magneto-dipole braking ( n = 3 ) ˙ P can be expressed as: ˙ P = αB P . (6) We thank prof. J.A. Pons for providing these data sets. Note, that theexact realisations used here are different from those presented in the paperby Popov et al. (2010). c (cid:13) , 000–000
MM and population synthesis of radio pulsars -20-18-16-14-12-10-3 -2 -1 0 1 2 Log P ˙ , s e c ond s / s e c ond s Log P, seconds -20-18-16-14-12-10-3 -2 -1 0 1 2
Log P ˙ , s e c ond s / s e c ond s Log P, seconds
Figure 11.
Models VIII (left panel) and IX (right panel). The data was calculated with the population synthesis model used in Popov et al. (2010). In the leftpanel pulsars evolve with constant magnetic fields. The K-S test probability for this model is 98 per cent. In the right panel results for decaying fields arepresented. For this model the K-S test probability is 96 per cent.
Here α = 8 π R / (3 Ic ) = 10 − cm s g − , R is the NS radius, I is the moment of inertia, and c is the speed of light. Therefore,the value of ˙ P does not follow the Gaussian distribution (in log-scale or not) from the very beginning of our simulations. If theevolution does not follow the magneto-dipole formula the argumentis also true for the models we used. In particular, it explains why theGMM method fails to describe with a single cluster the ensemblesof pulsars when the maximum duration of the evolutionary trackfollowed was short. For example, we performed the calculationsfor the maximum track duration equal to ∼ years. And inthis case the results cannot be explained by a single cluster.Simple replacement of ˙ P by B does not improve the situa-tion. Pulsars are rapidly evolving objects, and in the P – ˙ P planewe can find objects at different stages of their evolution. The tracksof pulsars in the P – ˙ P plane in the simple case of the constanteffective magnetic field are straight lines. The length of the line isdetermined by the age, the initial spin period and the initial mag-netic field. Therefore, the longer lines are in the upper part of the P – ˙ P plot. During the evolution, the Gaussian distribution of pul-sars in the P – B plane shifts and rotates. In addition to the alreadyevolved pulsars, there are new ones being born constantly. Conse-quently, the distribution would not follow a Gaussian distributioneven in the case of P – B plot.L12 discussed the problem of robustness of the method in theirstudy, too. The authors show that if 3 per cent of all pulsars are ran-domly removed from the dataset, no significant changes appear inthe structure of the cluster distribution. In our opinion, 3 per centis a small number. New observations routinely bring many tens,or even hundreds of newly discovered objects. So, 10 per cent isa more realistic number to study the stability of the method. Itis shown above that a 10 per cent modification of the number ofpulsars changes the results significantly. In addition, as we havedemonstrated, a systematic exclusion of even a small numer of ob-jects (magnetars etc. in the studied case) also changes the set ofthe ellipses significantly. This brings us to the conclusion that themethod is not very effective. In our opinion, the GMM is not very effective in distinguishingphysically or evolutionary related groups when applied to rapidly evolving populations observed with significant selection effects.That is why we decided to apply the method to more “stable” setsof data. First, we decided to use the GMM-code to study the well-known distribution of the gamma-ray bursts (GRBs) in theduration-hardness diagram. It is well established that there areat least two populations: the long soft and the short hard bursts(Klebesadel 1992; Kouveliotou et al. 1993). The duration and hard-ness of GRBs data are obtained from the BATSE 4B Gamma-RayBursts Catalog (Paciesas et al. 1999). We applied the GMM methodto the distributions of bursts by logarithm of duration ( t ) and bylogarithm of hardness. The hardness is defined as S = F −
300 keV F −
100 kev , (7)where F −
300 keV is the flux in the energy range 100-300 keV,and F −
100 kev is the flux in the range 50-100 keV.Naively, one expects that the method will easily describe thedistribution with two clusters which correspond to the two GRBtypes. But this is not the case! In Figure 12 we see that the methodrequires four clusters. Depending on the initial guess, a differentconfiguration of ellipses with nearly equal likelihood can appear,and in some of them the ellipses overlap.For another test we choose the high-mass X-ray binaries, inparticular the Be/X-ray systems. For these, the existence of twotypes have been established by Knigge et al. (2011). We want tocheck, whether the Expectation-Maximization algorithm can alsoidentify this dichotomy. The data on Be/X-ray binary systems istaken from the catalog by Raguzova & Popov (2005) .The GMM method was applied to the distribution of Be/X-ray pulsars by logarithm of the spin period and the logarithm ofthe orbital period. The results are shown in Figure 13. Indeed, twoGaussians are enough to describe the data. However, we do not seethat the method separates sources into long (both, the orbital andthe spin) and the short periods. Instead, the short period sourcesare united with those with the longest periods. We have to note In L12 it was demonstrated that the method can be successfully used forthe millisecond pulsars, which evolve much more slowly compared to thenormal radio pulsars. The catalog is available on-line athttp://xray.sai.msu.ru/ ∼ raguzova/BeXcat/c (cid:13) , 000–000 A. Igoshev and S. Popov -1.5-1-0.5 0 0.5 1 1.5 2-2 -1 0 1 2 3 H a r dne ss Log t , seconds -1.5-1-0.5 0 0.5 1 1.5 2-2 -1 0 1 2 3 H a r dne ss Log t , seconds Figure 12.
Gamma-ray bursts analysis. Two realizations of the Expectation-Maximization method are shown. In both cases, the well-known bimodal distri-bution in the duration-hardness cannot be described by two Gaussian clusters. The probability that the presented model and the data are drawn from the sameparent distribution is 99 per cent for the model in the left panel, and 97 per cent for the model in the right panel.
Log P o r b , da ys Log P spin , seconds 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 3.5 4
Log P o r b , da ys Log P spin , seconds
Figure 13.
Analysis of the P spin – ˙ P orb distribution of Be/X-ray binaries with the GMM method. The probability that the model and data are drawn from thesame parent distribution is 99.9 per cent for both models. here, that the objects with the longest periods that we use were notincluded in the sample in Knigge et al. (2011). In our study, the GMM method was at first applied to the distribu-tions of the known normal radio pulsars in the P – ˙ P plane. It isfound that the exclusion of magnetars, the thermally emitting NSs,etc. modifies the set of the clusters as compared to those foundby L12. The pulsars detected in the Parkes multibeam survey maybe described with just three clusters. It was also found that the re-sults of the Gaussian cluster finding procedure are not robust (orthe GMM method is oversensitive and the small changes in the datacause significant changes in the results). The relative position of theellipses were changed when we excluded random 10 per cent of thepulsars from the observational data set.In the second part of our study we generated ensembles ofpulsars using population synthesis models of different complexity.First, we find that the GMM method is strictly dependent on thetotal number of pulsars in the analyzed ensemble. The choice ofthe electron density model also has a strong influence on the clusterdistribution. The spiral structure of our Galaxy has a smaller effect. Such features as the bimodal distribution of the initial parametersare hardly recognized by the GMM method for the realistic choicesof parameters. The magnetic field decay changes the distribution ofthe clusters. Typically, if the field is decaying it is necessary to usemore clusters to describe the data.We conclude that the GMM is not effective to test models ofthe normal radio pulsar evolution because of the method’s over-sensitivity. ACKNOWLEDGMENTS
We thank Kejia Lee for numerous useful comments during thewhole work on this study and on the text of the manuscript. Specialthanks to Vasily Belokurov, who carefully read the paper whichhelped to improve its style. We also thank the referee for usefulremarks and suggestions. The work of S.P. was supported by theRFBR grant 12-02-00186. The work of A.P. was supported by SaintPetersburg University grant 6.38.73.2011. c (cid:13) , 000–000 MM and population synthesis of radio pulsars REFERENCES
Cordes J. M., Lazio T. J. W., 2002, arXiv: astro-ph/0207156Faucher-Gigu`ere C.-A., Kaspi V. M., 2006, ApJ, 643, 332Hewish A., Bell S. J., Pilkington J. D. H., Scott P. F., Collins R. A.,1968, Nature, 217, 709Kim S. C., Kang T. J., 2007, Pattern recognition, 40, 1207Klebesadel R. W., 1992, in Ho C., Epstein R. I., Fenimore E. E.,eds, The durations of gamma-ray bursts Gamma-ray bursts - ob-servations, analyses and theories. Cambridge University Press,Cambridge, pp 161–168Knigge C., Coe M. J., Podsiadlowski P., 2011, Nature, 479, 372Kouveliotou C., Meegan C. A., Fishman G. J., Bhat N. P., BriggsM. S., Koshut T. M., Paciesas W. S., Pendleton G. N., 1993, ApJ,413, L101Kuijken K., Gilmore G., 1989, MNRAS, 239, 651Lee K. J., Guillemot L., Yue Y. L., Kramer M., Champion D. J.,2012, MNRAS, 424, 2832Lorimer D. R., Faulkner A. J., Lyne A. G., Manchester R. N.,Kramer M., McLaughlin M. A., Hobbs G., Possenti A., StairsI. H., Camilo F., Burgay M., D’Amico N., Corongiu A., Craw-ford F., 2006, MNRAS, 372, 777Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993Manchester R. N., Lyne A. G., Camilo F., Bell J. F., Kaspi V. M.,D’Amico N., McKay N. P. F., Crawford F., Stairs I. H., PossentiA., Kramer M., Sheppard D. C., 2001, MNRAS, 328, 17Paciesas W. S., Meegan C. A., Pendleton G. N., 1999, ApJS, 122,465Popov S. B., Pons J. A., Miralles J. A., Boldin P. A., Posselt B.,2010, MNRAS, 401, 2675Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P.,2007, Numerical recipes: the art of scientific computing. Cam-bridge University Press, CambridgeRaguzova N. V., Popov S. B., 2005, Astronomical and Astrophys-ical Transactions, 24, 151Rawley L. A., Taylor J. H., Davis M. M., 1986, Nature, 319, 383Reynolds D. A., Quatieri T. F., Dunn R. B., 2000, Digital signalprocessing, 10, 19Ruderman M. A., Sutherland P. G., 1975, ApJ, 196, 51Vall´ee J. P., 2005, AJ, 130, 569Vivekanand M., Narayan R., 1981, Journal of Astrophysics andAstronomy, 2, 315Vranesevi´c N., Melrose D. B., 2011, MNRAS, 410, 2363 c (cid:13)000