Classification of pulsars with Dirichlet process Gaussian mixture model

TL;DR

This paper applies a Dirichlet process Gaussian mixture model to classify pulsar types based on their period and period derivative, revealing insights into their evolutionary relationships and physical properties.

Contribution

The study introduces an unsupervised machine learning approach using DPGMM to classify pulsar populations and analyze their physical characteristics in a unified framework.

Findings

DPGMM effectively clusters pulsar types in parameter space.

Distinct pulsar classes show different characteristic ages and magnetic fields.

Results support models of pulsar evolution and magneto-thermal processes.

Abstract

Young isolated neutron stars (INS) most commonly manifest themselves as rotationally powered pulsars (RPPs) which involve conventional radio pulsars as well as gamma-ray pulsars (GRPs) and rotating radio transients (RRATs). Some other young INS families manifest themselves as anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs) which are commonly accepted as magnetars, i.e. magnetically powered neutron stars with decaying superstrong fields. Yet some other young INS are identified as central compact objects (CCOs) and X-ray dim isolated neutron stars (XDINSs) which are cooling objects powered by their thermal energy. Older pulsars, as a result of a previous long episode of accretion from a companion, manifest themselves as millisecond pulsars and more commonly appear in binary systems. We use Dirichlet process Gaussian mixture model (DPGMM), an unsupervised machine learning algorithm, for analyzing the distribution of these pulsar families in the parameter space of period and period derivative. We compare the average values of the characteristic age, magnetic dipole field strength, surface temperature and transverse velocity of all discovered clusters. We verify that DPGMM is robust and provides hints for inferring relations between different classes of pulsars. We discuss the implications of our findings for the magneto-thermal spin evolution models and fallback discs.

Tables

Table 1. Table 1: The estimated parameters of DPGMM ( m max = 10 subscript 𝑚 max 10 m_{\rm max}=10 , α = 10 − 5 𝛼 superscript 10 5 \alpha=10^{-5} ).   w subscript 𝑤 absent w_{\rm~{}} , 𝝁 𝝁 \boldsymbol{\mu} and 𝚺 𝚺 \boldsymbol{\Sigma} indicate weight, location and shape of the corresponding Gaussian distribution, respectively.

DPGMM Parameters
Cluster	$w_{}$	$𝝁$	$𝚺$
1st	0.0787
2nd	0.0384
3rd	0.3344
4th	0.4111
5th	0.1129
6th	0.0244

Table 2. Table 2: The distribution of the pulsar families over clusters.

	Cluster
Pulsar Family	1st	2nd	3rd	4th	5th	6th	Total
Normal	137	73	703	906	142	21	1982
GRP	38	2	0	15	62	0	117
RRAT	0	0	21	12	1	1	35
AXP	0	0	4	0	0	10	14
SGR	0	0	0	0	0	8	8
XDINS	0	0	6	0	0	1	7
CCO	0	0	3	0	0	0	3
Total	175	75	737	933	205	41	2166

Table 3. Table 3: The distribution binary pulsar families over clusters.

Cluster	Normal	GRPs	Total
1st	103	32	135
2nd	51	0	51
3rd	13	0	13
4th	6	0	6
5th	3	1	4
6th	0	0	0
Total	176	33	209

Table 4. Table 4: The distribution of the companion types of binary pulsars over clusters. Companion types — main-sequence (MS), neutron star (NS), carbon-oxygen (or ONeMg) white dwarf (CO), helium white dwarf (HE), and ultra-light companion or planet (UL) — are obtained from Australia Telescope National Facility (ATNF) Pulsar Catalog (http://www.atnf.csiro.au/research/pulsar/psrcat/).

Cluster	MS	NS	CO	HE	UL
1st	9	0	15	77	23
2nd	0	14	17	12	4
3rd	0	2	3	4	2
4th	4	1	1	0	0
5th	2	1	0	1	0
6th	0	0	0	0	0
Total	15	18	36	94	29

Table 5. Table 5: The characteristic age (year) and dipole magnetic field strength (Gauss) of clusters.

	Characteristic Age	Surface Magnetic Flux
Cluster	(yr)	(Gauss)
1st	$4.49\times {10}^9$	$2.52\times {10}^8$
2nd	$1.35\times {10}^9$	$2.96\times {10}^9$
3rd	$3.06\times {10}^7$	$7.04\times {10}^{11}$
4th	$2.70\times {10}^6$	$1.28\times {10}^{12}$
5th	$1.15\times {10}^5$	$2.21\times {10}^{12}$
6th	$6.75\times {10}^4$	$3.35\times {10}^{13}$

Table 6. Table 6: Statistics of transverse velocity v T subscript 𝑣 T v_{\rm T} (km s -1 ) of non-recycled pulsars (the 3 rd , 4 th , 5 th  and 6 th  clusters). Here N v T subscript 𝑁 subscript 𝑣 T N_{v_{\rm T}} stands for the number of pulsars with the measured speed in each cluster. v T ¯ ¯ subscript 𝑣 T \overline{v_{\rm T}} , v T ∗ ¯ ¯ superscript subscript 𝑣 T \overline{v_{\rm T}^{*}} , σ v T subscript 𝜎 subscript 𝑣 T \sigma_{v_{\rm T}} and σ v T ∗ subscript superscript 𝜎 subscript 𝑣 T \sigma^{*}_{v_{\rm T}} stand for mean, weighted mean, standard deviation and weighted standard deviation of the transverse speed, respectively.

	Transverse Velocity $v_T$ (km s-1)
Cluster	$N_{v_T}$	$\overline{v_T}$	$\overline{v_T^{\ast }}$	${\sigma }_{v_T}$	${\sigma }_{v_T}^{\ast }$
3rd	50	335.5	363.3	625.4	677.4
4th	71	470.0	373.6	753.6	570.1
5th	20	302.6	359.1	248.8	265.4
6th	6	252.2	231.5	141.0	79.7

Table 7. Table 7: Statistics of blackbody temperature (eV) of non-recycled pulsars (the 3 rd , 4 th , 5 th  and 6 th  clusters). Here N T bb subscript 𝑁 subscript 𝑇 bb N_{T_{\rm bb}} stands for number of pulsars with the measured temperatures. T bb ¯ ¯ subscript 𝑇 bb \overline{T_{\rm bb}} , T bb ∗ ¯ ¯ superscript subscript 𝑇 bb \overline{T_{\rm bb}^{*}} , σ T bb subscript 𝜎 subscript 𝑇 bb \sigma_{T_{\rm bb}} and σ T bb ∗ subscript superscript 𝜎 subscript 𝑇 bb \sigma^{*}_{T_{\rm bb}} stand for mean, weighted mean, standard deviation and weighted standard deviation of temperature, respectively.

	Blackbody temperature $T_{\mathrm{bb}}$ (eV)
Cluster	$N_{T_{\mathrm{bb}}}$	$\overline{T_{\mathrm{bb}}}$	$\overline{T_{\mathrm{bb}}^{\ast }}$	${\sigma }_{T_{\mathrm{bb}}}$	${\sigma }_{T_{\mathrm{bb}}}^{\ast }$
3rd	14	219.2	138.2	175.0	83.7
4th	3	128.7	126.3	67.7	55.4
5th	14	172.0	157.5	84.4	51.0
6th	19	395.5	334.1	166.5	171.3

Equations

$$p({\mathbf{x}}_{\mathrm{j}},z_{\mathrm{j}}=i)=p({\mathbf{x}}_{\mathrm{j}}|z_{\mathrm{j}}=i)p(z_{\mathrm{j}}=i).$$

$$p({\mathbf{x}}_{\mathrm{j}})=\sum _{i=1}^{m}p({\mathbf{x}}_{\mathrm{j}}|z_{\mathrm{j}}=i)p(z_{\mathrm{j}}=i).$$

$$\mathcal{N}\left({\mathbf{x}}_{\mathrm{j}}|{\mathbf{\mu }}_{\mathrm{i}},{\mathbf{\Sigma }}_{\mathrm{i}}\right)=\frac{\exp \left(-\frac{1}{2}({\mathbf{x}}_{\mathrm{j}}-{\mathbf{\mu }}_{\mathrm{i}}{)}^{\mathrm{T}}{{\mathbf{\Sigma }}_{\mathrm{i}}}^{-1}({\mathbf{x}}_{\mathrm{j}}-{\mathbf{\mu }}_{\mathrm{i}})\right)}{(2\pi {)}^{d/2}|{\mathbf{\Sigma }}_{\mathrm{i}}{|}^{1/2}},$$

$$\{\mathbf{\mu },\mathbf{\Sigma }\}\sim \mathcal{N}\mathcal{I}\mathcal{W}\left({\mathbf{\mu }}_0,\beta ,\mathbf{\psi },\nu \right),$$

$$\mathbf{w}\sim DP(\alpha ,\mathcal{H}),$$

$$z_{\mathrm{j}+1}=\underset{i}{\mathrm{argmax}}p({\mathbf{x}}_{\mathrm{j}+1},z_{\mathrm{j}+1}=i).$$

$$B_{\mathrm{LC}}=\sqrt{\frac{24{\pi }^{4}I\dot{P}}{c^3{P}^5}}.$$

$$\log (B_{\mathrm{LC}})$$

$$\log (L_{\mathrm{sd}})$$

$$B_{\mathrm{d}}={\left(\frac{3c^{3}I}{8{\pi }^2{R}^6}\right)}^{1/2}\sqrt{P\dot{P}}$$

Full text

††pubyear: 2019††pagerange: Classification of pulsars with Dirichlet process Gaussian mixture model–Classification of pulsars with Dirichlet process Gaussian mixture model

1 Introduction

Rotationally powered pulsars (RPPs), the most common manifestation of neutron stars, were discovered first as radio pulsars (Hewish et al., 1968). More than $\sim 2600$ RPPs discovered to date are understood to be strongly magnetized ($B\sim 10^{12}$ G), rapidly rotating (with spin periods $P\sim 0.1$ s) neutron stars spinning down by torques due to magnetic dipole radiation and particle emission. These objects emit most of the energy, tapped from their rotational kinetic energy, in X- and gamma-ray bands. The total radiative luminosity of RPPs is thus below their spin-down power $L_{\rm sd}\equiv-I\Omega\dot{\Omega}$ where $I$ is the moment of inertia of the compact object, $\Omega=2{\mathrm{\pi}}/P$ is the spin angular frequency and $\dot{\Omega}$ is its time derivative (see Manchester, 2017, for a review). The bulk of the RPP population are the classical radio pulsars with periods $P=0.2-2$ s, period derivatives $\dot{P}=10^{-16}-10^{-13}$ s s*-1*, magnetic fields $B\sim 10^{10}-10^{13}$ G and characteristic ages $\tau_{\rm c}\equiv P/2\dot{P}\sim 10^{3}-10^{7}$ years. Some of the RPPs are associated with supernova remnants (SNRs). Young RPPs are usually isolated objects in the sense that they have no binary companion. The pulsations are attributed to the beaming of the emitted radiation and the “lighthouse” effect produced by the spin of the object. Thus, apart from distant and dim objects, there could be many RPPs whose beam does not sweep our line of sight.

A population among the RPPs are the millisecond pulsars (MSP; Backer et al., 1982) which are old but “recycled” pulsars (Alpar et al., 1982); these objects are understood to descent from low-mass X-ray binaries (LMXBs; see Bhattacharya & van den Heuvel, 1991; Tauris & van den Heuvel, 2006, for reviews) where the neutron star spins-up by accreting matter from a companion via a disc. The accreted matter transfers angular momentum to the neutron star thereby spinning it up to millisecond periods and the magnetic field of the neutron star is reduced to values $B<10^{9}$ G in the process.

More than 200 $\gamma$-ray pulsars (GRPs) discovered by the Large Area Telescope onboard Fermi (Fermi-LAT; Atwood et al., 2009) constitute a subclass of RPPs (Abdo et al., 2010a, 2013). All known GRPs have $L_{\rm sd}>10^{33}$ erg s*-1*. About half of the GRPs detected by Fermi are radio-quite. This possibly indicates that the “radio beam” which is narrower than the “gamma beam” does not pass from our line of sight. Approximately half of the discovered GRPs are MSP which was unexpected before the discovery (Abdo et al., 2009).

Another recently discovered family of RPPs is the rotating radio transients (RRATs) identified in the Parkes multibeam survey (McLaughlin et al., 2006). Unlike ordinary pulsars detected by searches in the frequency domain, RRATs show sporadic emission detected through their bright single pulses. These objects show bursts of duration $2-30$ ms with an interval in the range 4 min-3 hr. Their spin periods are in the range $0.4-7$ s. The period derivatives are measured in 3 sources implying magnetic fields $B\sim 10^{12}-10^{14}$ G and characteristic ages $\tau_{\rm c}\sim 0.1-3$ Myr.

In the last two decades many young isolated neutron star (INS) families (see Popov, 2008; Kaspi, 2010; Harding, 2013; Safi-Harb, 2017, for reviews) other than radio pulsars are identified. These objects are usually radio-quiet and have X-ray luminosities exceeding the spin-down power of the compact object. This indicates the availability of energy budgets other than rotational kinetic energy and the possibility of evolutionary paths other than those leading to RPPs.

Of these groups of objects soft gamma ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) are commonly assumed to be magnetars (see Woods & Thompson (2006); Mereghetti (2008); Rea (2014); Turolla et al. (2015); Mereghetti et al. (2015); Kaspi & Beloborodov (2017); Gourgouliatos & Esposito (2018) for reviews) i.e. INSs with large dipole magnetic fields $B\sim 10^{14}-10^{15}$ G as inferred from their rapid spin down $\dot{P}\sim 10^{-13}-10^{-11}$ s s*-1* (Kouveliotou et al., 1998) and slow periods clustered between $P=2-12$ seconds. According to the magnetar model (Duncan & Thompson, 1992; Thompson & Duncan, 1996) the persistent X-ray emission of these objects with X-ray luminosity of $L_{\rm X}\sim 10^{35}-10^{36}$ erg s*-1* is powered by the decay of this strong magnetic field (Paczynski, 1992; Thompson & Duncan, 1995) in excess of the quantum critical limit $B_{\rm QED}=4.4\times 10^{13}$ G. These objects occasionally show super-Eddington outbursts and, very rarely, giant bursts which are addressed in the magnetar model with the breaking of the neutron star crust due to magnetic stresses and reconfiguration of their fields, respectively. The magnetars are young objects as implied by their characteristic ages $\tau_{\rm c}\sim 10^{4}$ years and about half of them being associated with SNRs.

X-ray dim isolated neutron stars (XDINSs), or sometimes called “magnificent seven”, are the 7 nearby neutron stars identified through their thermal X-ray emission (Özel, 2013; Potekhin et al., 2015; Mereghetti, 2011) with luminosities of order $L_{\rm X}\sim 10^{30}-10^{32}$ erg s*-1* (see Haberl, 2007; Kaplan, 2008; Turolla, 2009, for reviews). They have a period range similar to the AXP/SGR family (Hambaryan et al., 2017), but are typically older, with characteristic ages $\tau_{\rm c}\sim 10^{5}-10^{6}$ years and kinematic ages of a few $10^{6}$ years (Tetzlaff et al., 2010, 2011, 2012). They have inferred dipole magnetic fields of $B\sim 10^{13}$ G, an order of magnitude lower than magnetars and an order of magnitude higher than conventional pulsars, but their surface magnetic fields as inferred from the narrow absorption features are 7 (Borghese et al., 2015) and 5 (Borghese et al., 2017) times larger in the case of RX J0720.4–3125 and RX J1308.6+2127, respectively.

Yet another young NS family is the central compact objects (CCOs; see De Luca, 2017, for a review) in supernova remnants (SNRs). These $\sim 10$ objects show no sign of RPP activity. Their X-ray spectra is dominated by the thermal emission showing some similarities with magnetars yet they are two orders of magnitude less luminous ($\sim 10^{33}-10^{34}$ erg s*-1*). The periods measured from 3 CCOs are in the range $0.1-0.4$ s and the measured period derivatives imply that the dipole fields of these objects are in the range $B\sim 10^{10}-10^{11}$ G (Gotthelf & Halpern, 2007, 2009; Halpern & Gotthelf, 2011), and hence sometimes are called “anti-magnetars” (Halpern & Gotthelf, 2010; Gotthelf et al., 2013). Yet there is evidence that these objects have much stronger “hidden” magnetic fields (Viganò & Pons, 2012; Viganò et al., 2013; Torres-Forné et al., 2016) as implied by the highly anisotropic emission (Shabaltas & Lai, 2011) leading to the observed high pulsed fraction. According to the “field burial scenario” (Muslimov & Page, 1995; Young & Chanmugam, 1995; Geppert et al., 1999; Ho, 2011; Bernal et al., 2012; Igoshev et al., 2016) this is due to an initial fallback accretion episode the nascent neutron star had suffered soon after the SN explosion. An exception among the CCO family is 1E 161348–5055 in SNR RCW 103 which has an unusually long period of 6.7 hours (de Luca, 2008) and has recently shown magnetar-like bursts (Rea et al., 2016).

Recent discoveries are blurring the borders of the classes: (i) The existence of high magnetic field RPPs (Pivovaroff et al. (2000); see Ng & Kaspi (2011) for a review) that have magnetic field strengths comparable and, in some cases, exceeding that of some magnetars and some of which has even shown magnetar-like X-ray bursts (Gavriil et al., 2008; Archibald et al., 2016; Göğüş et al., 2016); (ii) radio detection from some magnetars (Camilo et al., 2006, 2007; Levin et al., 2010); (iii) identification of low-magnetic field magnetars (Rea et al., 2010) i.e. INS showing magnetars bursts but with ordinary inferred dipole field strengths (see Turolla & Esposito, 2013, for a review). These discoveries indicate that the dipole field strength is not the single parameter leading to the different manifestations favouring the early suggestions (Gavriil et al., 2002; Ekşi & Alpar, 2003; Ertan & Alpar, 2003; McLaughlin et al., 2003) that what causes the magnetar activity could be in the higher multipoles (Alpar et al., 2011; Rodríguez Castillo et al., 2016).

In fact there are other observations suggesting that AXP/SGRs have “low” dipole fields: (i) They are not detected in the Fermi/LAT observations (Şaşmaz Muş & Göğüş, 2010; Abdo et al., 2010b) in the GeV range though, according to the outer gap model (Zhang & Cheng, 1997), they are expected (Cheng & Zhang, 2001) to emit high energy gamma-rays should they have super-strong magnetic dipole fields $\gtrsim 10^{14}$ G (Tong et al., 2011). (ii) their transverse velocities are measured (Helfand et al., 2007; Deller et al., 2012; Tendulkar et al., 2012, 2013) to be $200\pm 100$ km s*-1*, similar to the velocities of RPPs though they are expected to have exceptionally large space velocities, $\sim 1000$ km s*-1*, due to getting stronger kicks (Duncan & Thompson, 1992; Thompson & Duncan, 1993) via the rocket propulsion effect should they have super-strong dipole fields. (iii) They do not commonly show the signature of the expected strong energy injection in the SNRs due to hosting a neutron star with a birth period of milliseconds (Vink & Kuiper, 2006; Martin et al., 2014; Borkowski & Reynolds, 2017) (but see Torres (2017)) that could produce super-strong fields by dynamo action (Duncan & Thompson, 1992; Thompson & Duncan, 1993).

Given the transitivity among the INS families the question naturally arises whether the classification of INS above is robust. Could we expect, for example, that more of the RPPs to show magnetar-like behaviour, or that more magnetars to appear in the radio band? What parameters could be leading to the different manifestations for objects with similar period and period derivatives? The magneto-thermal evolutionary theory (e.g. Viganò & Pons, 2012; Viganò et al., 2013) addresses this diversity by the presence of toroidal magnetic fields hidden in the crust as an extra parameter shaping the lives of pulsars. The fallback disc model (Alpar, 2001; Chatterjee et al., 2000) invokes the mass and specific angular momentum (Ertan et al., 2007, 2009) of a putative supernova fallback disc as two parameters leading to the diversity.

A well-known unsupervised machine learning algorithm, Gaussian mixture model (GMM; Press et al., 2007), had been used by Lee et al. (2012) to analyze the distribution of pulsars on the $P-\dot{P}$ diagram and identified six Gaussian clusters (2 for millisecond—recycled—pulsars and 4 for young pulsars). Igoshev & Popov (2013) showed that GMM is over-sensitive to the data and does not demonstrate robust clustering performance. There are more advanced variations of GMM such as Dirichlet process Gaussian mixture model (DPGMM). This model has the advantage over GMM that the number of distributions are determined automatically though it is not used as widely as GMM especially by the astronomy community. In addition its generalization performance is better than GMM, as Bayesian methods are more robust against overfitting (Rasmussen, 2000; Heller, 2008; Witten et al., 2016). We have recently showed that DPGMM is better than GMM in both clustering and classifying pulsars in the $P-\dot{P}$ parameter space (Ay et al., 2019).

The organization of our paper is as follows. In § 2 we review DPGMM employed in this work. In § 3 we present our results and finally, in § 4, we discuss the implications of our findings.

2 Method

2.1 Dirichlet process Gaussian mixture model

Here we will briefly sketch the method employed in this work. The detailed information for DPGMM can be found in (Shin et al., 2009; Del Pozzo et al., 2018).

In finite mixture models, it is assumed that the data consists of a certain number of clusters in which each cluster is generated by a probability distribution. Let $\boldsymbol{X}^{\rm T}=\{\boldsymbol{x_{\rm j}}\}_{\rm j=1}^{\rm N}$ be the collection of $d$-dimensional observed data with $N$ number of instances and assume further that there are $m$ distributions that generate $\boldsymbol{X}$. Because it is not known by which distribution a certain data point is generated, it is useful to define a latent data $\boldsymbol{Z}$ for distribution assignments; $\boldsymbol{Z}^{\rm T}=\{z_{\rm j}\}_{\rm j=1}^{\rm N}$ where $z_{\rm j}=i$ means that $\boldsymbol{x_{\rm j}}$ is generated by the i th distribution, $j\in\{1,2,\ldots,N\}$ and $i\in\{1,2,\ldots,m\}$.

The joint probability of a data point and i th distribution is defined as

[TABLE]

The joint probability value indicates how likely a data point belongs to a cluster. Because $\boldsymbol{x_{\rm j}}$ must be populated by one of the $m$ distributions, we can obtain marginal probability of $\boldsymbol{x_{\rm j}}$ from equation (1) as follows

[TABLE]

The marginal probability value can be considered as the probability that a data point belongs to any cluster. On the right side of the equation, the second term is the probability of the i th distribution which is also known as the weight of the i th cluster such that $w_{\rm i}\equiv p(z_{\rm j}=i)$ where $\sum_{i=1}^{m}{w_{\rm i}}=1$. Besides, if we assume that each cluster is distributed normally, then the first term on the right side of the equation is defined as Gaussian distribution;

[TABLE]

where $\boldsymbol{\mu_{\rm i}}$ is $d$-length mean vector and $\boldsymbol{\Sigma_{\rm i}}$ is $d\times d$ covariance matrix of the i th Gaussian distribution. In this case, the finite mixture model is called as Gaussian mixture model (GMM).

In Bayesian approach, GMM parameters $\boldsymbol{\theta}=\{\boldsymbol{w},\boldsymbol{\mu},\boldsymbol{\Sigma}\}$ are also modeled with probability distributions. The choice of distribution families is based on conjugate prior principle (Raiffa & Schlaifer, 1961) for the computational simplicity. Accordingly, the prior distributions for the parameters of Gaussian distributions are selected as Normal-Inverse-Wishart distribution;

[TABLE]

where $\{\mathcal{\boldsymbol{\mu_{\rm 0}}},\beta,\boldsymbol{\psi},\nu\}$ are mean, scaling factor, scale matrix and degrees of freedom. Furthermore, if we assume that the cluster weights are Dirichlet process (DP; Ferguson, 1973) distributed

[TABLE]

then the mixture model is called as Dirichlet process Gaussian mixture model (DPGMM). In the equation, $\alpha$ is the concentration parameter and $\mathcal{H}$ is the base distribution for DP. Because a realization from DP is a infinite length probability vector, DPGMM is a non-parametric, i.e. infinite, mixture model.

In this study, we employed scikit-learn (Pedregosa et al., 2011) implementation of DPGMM that uses variational inference (Attias, 2000; Bishop, 2006; Blei & Jordan, 2006) for learning model parameters. We also benefited from many open source software for scientific computations, analyzes and visualizations, such as scipy (Virtanen et al., 2019), numpy (Oliphant, 2006), pandas (McKinney et al., 2010), matplotlib (Hunter, 2007), seaborn (Waskom et al., 2017) and scikit-image (van der Walt et al., 2014). In addition, Spyder111https://www.spyder-ide.org/ and IPython (Pérez & Granger, 2007) were used as development environments.

2.2 Data Collection

We examine pulsar distribution on $P-\dot{P}$ parameter space in logarithmic scale. In order to increase data samples as much as possible, various neutron star catalogs are combined. Most of the pulsar data are obtained from Australia Telescope National Facility (ATNF) Pulsar Catalog222http://www.atnf.csiro.au/research/pulsar/psrcat/ that contains up to 2658 instances with numerous source types. In addition, 30 instances of SGRs and AXPs (either confirmed or candidates) from McGill Online Magnetar Catalog333http://www.physics.mcgill.ca/p̃ulsar/magnetar/main.html (Olausen & Kaspi, 2014) and 40 instances of other thermally emitting neutron stars such as CCOs and XDINSs from (Viganò et al., 2013)444http://www.neutronstarcooling.info/ are collected. Additionally, 107 RRATs from RRATOLOG555http://astro.phys.wvu.edu/rratalog/ (McLaughlin et al., 2006; Keane et al., 2010). and 117 GRPs from The Second Fermi-LAT catalogue666https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermil2psr.html (Abdo et al., 2013) are included. Finally, a recently discovered (Tan et al., 2018) RPP with an exceptional long period is also included to the data set. As a result, we obtained an amount of 2166 INSs with the observed spin parameters by combining all of these data sources.

3 Results

3.1 Application of DPGMM

DPGMM has a relatively large number of hyper-parameters to be fixed. For example, although DPGMM is an infinite mixture model, scikit-learn implementation of DPGMM requires the hyper-parameter $m_{\rm max}$ that corresponds to the maximum number of clusters to be set to an appropriate value as it would not be practical to consider infinite number of clusters during fitting. We have set the value of $m_{\rm max}$ as 10 as we did not expect to find more than 10 clusters given the number of pulsar families and the 6 clusters already determined by GMM (Lee et al., 2012).

Another hyper-parameter is the concentration parameter $\alpha$ that determines how sparse distributions are allowed to be and thus may effect the number of clusters the model can find. For a given $\alpha$, the optimum number of clusters that best represent the data among the $m_{\rm max}$ clusters and their parameters are determined by DPGMM during the fitting process. We have seen that the value of $\alpha$ in the literature varies significantly: e.g. Haines & Xiang (2013) set $\alpha=0.01$ while Shin et al. (2009) set $\alpha=1$. In another example, Chen et al. (2015) have searched optimum value for $\alpha$ in a range between $0.001$ and $100$, and observed that the results are relatively same independent of different $\alpha$ values.

We created DPGMMs for different $\alpha$ values within a wide range from $10^{-10}$ to $10^{10}$. As a result, we have observed that DPGMM for this wide range of $\alpha$ values discovered either 6 or 7 clusters depending on whether $\alpha$ is lower or higher than $10^{4}$, respectively. We have also observed that the estimated model parameters are remarkably independent of $\alpha$ for all models that discovered 6 clusters ($\alpha\leq 10^{4}$). The same situation is also true for all models that discovered 7 clusters ($\alpha>10^{4}$). In short, we have obtained two different stable models; one having 6 clusters and the other 7 clusters. The recommended value of $\alpha$ is much less than 1 and for such values our model finds 6 clusters. We have employed Bayesian information criteria (Schwarz, 1978) and Akaike information criteria (Akaike, 1974) in a former study (Ay et al., 2019) to determine the optimum values of the clusters in which we found 6 clusters. This result also is consistent with the number of clusters found by Lee et al. (2012) who used multidimensional Kolmogorov-Smirnov test. We thus selected the model with 6 clusters as the optimal one and our results with 7 clusters are not shown here. In Table 1, the estimated parameters of DPGMM created with $m_{\rm max}=10$ and $\alpha=10^{-5}$ are shown.

3.2 Parameter space classification

In Figure 1, all clusters and decision boundaries on $P-\dot{P}$  parameter space revealed by DPGMM are illustrated. Bivariate Gaussian curves are represented by 2$\sigma$ confidence ellipses and each cluster region is enumerated and colored uniquely in here. The mixture models are inductive models so that after fitting a model to a data set, it is able to predict the cluster in which a ‘novel’ instance should belong. Formally, the cluster label of a novel instance $\boldsymbol{x_{\rm j+1}}$ is determined as

[TABLE]

If we classify the entire parameter space with this approach, the separating line between cluster regions naturally arise as shown in Figure 1 and in Figure 2. These lines are called as decision boundaries such that the cluster assignment by the model is different for the one side of the line and the other side.

In the case of mixture of Gaussian distributions, there is a possibility of being more than one region for a single cluster on classified parameter space. To illustrate, there are two regions for the 6th cluster in Figure 2; one is on the top left, and the other is on the bottom right. It is significant that all pulsars belonging to the 6th cluster are in the top left region while the bottom one does not cover any pulsars at all. Thus, it can be said that instead of the 6th cluster, the bottom region should be merged with the neighbor region corresponding to the 2nd cluster. Therefore, the separating line between the 2nd cluster and the bottom region of the 6th cluster seems incorrect. For this reason, we examined the dependability of boundary lines in this study.

As mentioned before, the joint probability defined in equation (1) and used in equation (6) for cluster assignments shows that how likely a data point $\boldsymbol{x_{\rm j}}$ is a member of i th cluster. Its value depends on the cluster weight and the multivariate Gaussian distribution defined in equation (3) where the exponent term in the numerator is actually half the negative of the square of Mahalanobis distance (Mahalanobis, 1936; De Maesschalck et al., 2000) that is a distance measure in terms of standard deviation. Roughly speaking, it can be said that the cluster assignment for a data point is determined by the weight of clusters and the Mahalanobis distance between the data point and cluster means. The cluster assignment is predominately determined by the weight of the clusters for nearby data points. On the other hand, the effect of Mahalanobis distance on the decision making of the model increases for the receding data points. As the Mahalanobis distance is not a deterministic measure, the cluster assignments may be incoherent for the distant regions. To illustrate, the bottom right region in Figure 2 is separated from the 2nd cluster and assigned to the 6th cluster by the model even though all pulsars belonging to the 6th cluster are in top left region. Therefore, we defined the dependability of the cluster assignments as the marginal probability values specified in equation (2) that can be considered as how likely a data point is a member of any cluster, as stated before. If the marginal probability of $\boldsymbol{x_{\rm j}}$ is lower, then the dependability of the cluster assignment for this data point is not reliable as well. In this way, we showed the degree of the dependability of border lines in Figure 1 and Figure 2 based on the reliability of data points over border lines, so that the border line is darker where it is more reliable, and it is fainter where it is less reliable. Accordingly, it is obvious that the color of the border line between the 2nd and the 6th clusters in Figure 2 is very faint that means it is not reliable as expected. In this way, we tried to eliminate the handicap of parameter space classification in the case of mixture of Gaussian distributions.

As we have discussed before, we have introduced 7 families of pulsars namely rotationally powered pulsars (RPPs) that involve millisecond pulsars, rotating radio transients (RRATs), gamma-ray pulsars (GRPs) as well as conventional radio pulsars, magnetars (AXPs, SGRs), X-ray dim INS (XDINSs) and central compact objects (CCOs). As seen in Figure 1, DPGMM identifies six clusters for the value of the concentration parameter $\alpha=10^{-5}$. In this section, we analyse whether these six clusters are linked to any parameters of pulsars other than $P$ and $\dot{P}$ (e.g. space velocity or surface temperature) that could attribute further meaning to the different clusters.

3.3 Pulsar distribution over the clusters

The distribution of pulsar types over DPGMM clusters is given in Table 2. Accordingly, the most of the ‘non-recycled’ RPPs detected in the radio band are located in the 3rd, 4th and 5th clusters. The 1st and the 2nd regions contain mostly millisecond pulsars. It is significant that millisecond pulsars are separated into two subclasses by DPGMM as in the case of GMM (Lee et al., 2012). In § 4, we discuss the possible origin of this division. The 2nd region is very large and extends beyond the pulsar death line. The 3rd region contains most of the RRATs and XDINSs, all CCOs, low-B magnetars (J1647–4552, J0418+5732, J2301+5852 and J1822–1604) and those RPPs close to the death line. There is again a 4th region populated by RPPs with a narrow $\dot{P}$ range. About $10\%$ of GRPs overflow to this region. The 5th region mostly contains ‘non-recycled’ GRPs while the recycled ones are in the 1st cluster. The 6th region is predominantly occupied by ’High-B‘ AXPs and SGRs excluding the low-B magnetars. Moreover, most of the pulsars are grouped in the 3rd and 4th clusters that are located at the center of the diagram. There are fewer pulsars in other clusters relatively.

3.4 Effect of binary companion for recycled pulsars

The distribution of binary pulsars over clusters is given in Table 3. Accordingly, the binary pulsars are heavily grouped in the 1st and the 2nd clusters where millisecond pulsars ‘live’. This, of course, is understood through the “recycling scenario” (Alpar et al., 1982) that address the very existence of these objects by spin-up of the neutron star through long-term accretion of matter and angular momentum from a binary companion. It is interesting, however, that both GMM and DPGMM found the millisecond pulsars are distributed into two clusters. What could be the underlying astrophysical distinction between these two clusters?

Almost all millisecond GRPs are located in the 1st cluster as seen in Figure 1 and in Table 3. This is obviously due to the higher rotational power, $L_{\rm sd}=4{\mathrm{\pi}}^{2}I\dot{P}/P^{3}$, of this cluster as well as the magnetic field strength at the light cylinder radius $R_{\rm LC}=cP/2{\mathrm{\pi}}$ given by

[TABLE]

This can more readily be inferred from Figure 2 where DPGMM classification of $B_{\rm LC}-L_{\rm sd}$ parameter space is illustrated.

As $\{L_{\rm sd},B_{\rm LC}\}$ are the linear functions of $\{P,\dot{P}\}$ in logarithmic scale such that

[TABLE]

then the data matrix $\boldsymbol{Y}=[\log\boldsymbol{B}_{\rm LC}\,\,\log\boldsymbol{L}_{\rm sd}]$ is also distributed with the mixture of Gaussians as is the original data $\boldsymbol{X}=[\log\boldsymbol{P}\,\,\log\boldsymbol{\dot{P}}]$. This is true for any data matrix $\boldsymbol{Y}$ which is linearly related with $\boldsymbol{X}$. Therefore, an analysis in the $B_{\rm LC}-L_{\rm sd}$ parameter space is equivalent to the analysis in the $P-\dot{P}$ diagram. For this reason, we applied the same DPGMM ($m_{\rm max}=10$, $\alpha=10^{-5}$) to the $B_{\rm LC}-L_{\rm sd}$ parameter space. The obtained pulsar classification in this plane is almost same as the classification in $P-\dot{P}$  parameter space as should be. As seen in Figure 2, $\gamma$-ray emission can not be the underlying astrophysical cause for the existence of two separate clusters for millisecond pulsars; it will rather be a consequence.

The companion types of some binary pulsars have been detected and the distribution of these companions over the clusters is shown in Table 4. Accordingly, there is a remarkable distinction between the properties of the companion objects in the two groups indicating to different astrophysical origins and recycling history. The 1st cluster is populated by millisecond pulsars with main sequence (MS), helium white dwarf (HE), carbon-oxygen white dwarf (CO) and ultra-light (UL) companions while lacking any neutron star companions. The 2nd cluster is populated by those with neutron star (NS) and white dwarf companions (CO and He), while having only 4 UL and no MS companions. Lee et al. (2012) attributes the identification of two separate clusters to the chemical composition of the companions. This implies that millisecond pulsars in the 2nd cluster systematically had larger mass companions that lived shorter. Their larger periods (hence the smaller rotational power and lack of gamma-emission) can be attributed to the shorter recycling history they suffered. It is remarkable that DPGMM distinguishes these two groups only through their spin parameter distributions.

3.5 Characteristic age and magnetic field

We calculated the characteristic age, $\tau_{\rm c}\equiv P/2\dot{P}$, and magnetic dipole field strength as inferred under the assumption of a orthogonal rotating magnetic dipole model

[TABLE]

of each cluster where we assumed $I=10^{45}\,{\rm g\,cm^{2}}$ and $R=10^{6}\,{\rm cm}$ for the moment of inertia and the radius of neutron stars, respectively. These parameters calculated based on the mean value ($\boldsymbol{\mu_{\rm i}}=\{P_{i},\dot{P}_{i}\}$) of corresponding Gaussian distribution. Therefore the results given in Table 5 can be assumed as the average values for each cluster. The table is arranged such that the age of cluster is increased with the cluster order as well as the magnetic dipole strength. Hence, the 1st cluster is the youngest one with the highest dipole field where the 6th cluster is the oldest one with the lowest dipole field. Beside that, the millisecond pulsars in the 2nd cluster have an order of magnitude larger magnetic dipole field strengths compared to those in the 1st cluster. According to our discussion in the previous subsection, stronger fields can be understood as a consequence of the shorter ‘recycling’ episode these systems passed through.

3.6 Transverse velocity

According to the field burial model (Muslimov & Page, 1995; Young & Chanmugam, 1995; Geppert et al., 1999; Ho, 2011), an initial, brief but intense, fallback accretion episode (Chevalier, 1989) following the supernova explosion is important in shaping the magnetic fields of nascent pulsars. The amount of matter accreted determines how much the initial magnetic field is buried under the crust and thus the time-scale of diffusion of the field to the surface. Space velocities of pulsars due to the kick they receive during their birth may thus determine how much mass the neutron star can accrete, as the accretion rate is inversely proportional to the cube of the velocity in Bondi-Hoyle accretion. The field diffusion time-scale as inferred from the measured braking indices of young pulsars and the space velocities are inversely correlated (Güneydaş & Ekşi, 2013; Ekşi, 2017) as expected from the field-burial scenario. The question then naturally arises whether the space velocity of pulsars have any role in the observed diversity of pulsars.

In order to check this possibility we showed the transverse velocity of non-recycled pulsars (i.e. those in the 3rd, 4th, 5th and 6th clusters in the $P-\dot{P}$ diagram) as seen in Figure 3. The figure is colored depending on the magnitude of the velocity so that the color changes with the magnitude from yellow to red sequentially. The transverse velocity statistics of same clusters are given in Table 6. Since the probability that each pulsar is a member of the cluster to which they belong is different, we calculated the weighted statistics of their transverse velocities by using their joint probability values (see equation (1)) with respect to corresponding Gaussian distributions as weights. Accordingly, the weighted standard error is very large due to the low number of space velocity measurements in each group and the results may not be statistically significant. We see, however, that the space velocities of pulsars in the 4th cluster is the largest. We also see that the average speeds of pulsars in the 6th cluster is about half of those in the 4th cluster. The velocity data was collected from the sources in which the spin parameters of pulsars were taken as specified in § 2.2.

3.7 Blackbody temperature

In Figure 4, the blackbody temperature of non-recycled pulsars are illustrated. The statistics of black body temperature of these clusters are given in Table 7. The weighted statistics are computed in the same way described in the previous section. According to these results, the temperature of the 6th cluster where magnetars are frequently hosted in is almost twice as large as temperatures of other clusters. In other words, AXP/SGRs are the hottest pulsar type. However, the results may not be statistically significant due to the small sample size. The blackbody temperatures of these pulsars are acquired from the sources mentioned in § 2.2 and also from Gonzalez et al. (2007); Kaspi & McLaughlin (2004); Webb et al. (2004); Hu et al. (2017); Li et al. (2005); Chang et al. (2011); Caraveo et al. (2010); Marelli et al. (2011); Güver et al. (2012); Posselt et al. (2012).

4 Discussion

We have examined the pulsar population with an unsupervised machine learning algorithm, Dirichlet process Gaussian mixture model (DPGMM). We confirmed the earlier result obtained with Gaussian mixture model (GMM) by Lee et al. (2012) that the millisecond pulsar population has two clusters and that the normal pulsar population consists of four clusters.

We have considered possible hidden parameters (space velocity and surface temperature) as possible underlying cause of the distinct clusters. We have first tried DPGMM on higher dimensional planes, e.g.  3-D planes ($P-\dot{P}-v_{\rm T}$ and $P-\dot{P}-T_{\rm bb}$) and 4-D plane ($P-\dot{P}-v_{\rm T}-T_{\rm bb}$). We have found that the results for the higher dimensions were unreliable due to the lack of the data instances available for the time being. We have then examined the transverse velocity and the blackbody temperature statistics for each cluster. We have found that the 4th cluster has about twice the average space velocity compared to the slowest cluster, the 6th one, (see Figure 3 and Table 6) providing a marginal evidence that space velocities of young neutron stars could play a role in their astrophysical manifestations by controlling the total mass accreted during post-supernova fallback and the depth of field burial (Güneydaş & Ekşi, 2013; Ekşi, 2017). This indicates the pulsars in the 4th cluster suffered less supernova fallback accretion and less field burial. Indeed the average magnetic field of pulsars in the 4th cluster are similar to those in the 5th cluster and are twice of those in the 3rd cluster. If the magnetic fields of pulsars decay as they age within a timescale less than a million years, the difference in the magnetic fields of the 4th and 3rd clusters could be also attributed to the former systems being younger and subject to less field decay. Yet it is unlikely that the magnetic fields of typical pulsars decay in such short timescales favouring the idea that the space velocities could be the underlying reason for the difference in the magnetic fields of the 4th and 3rd clusters.

We note that the very high space velocities of CCOs imply that they are not likely to suffer very strong post-supernova fallback accretion. This does not neatly fit into the ‘field burial scenario’ (Muslimov & Page, 1995; Geppert et al., 1999) as the small dipole fields of these objects, in this model, are attributed to very intense fallback accretion. This may indicate the role of the initial magnetic field and spin period as other parameters controlling the amount of the accreted fallback mass: The low magnetic fields and relatively slow initial periods of CCOs could allow them to accrete large amount of mass to bury their fields in spite of their large space velocities. Another solution high-speed-low-field puzzle could be that their dipole fields are in fact very strong but their rotation and spin axis are aligned to a very high precision such that the cluster of their dipole moment perpendicular to the rotation axis inferred from spin-down is so small. This idea, however, would work only in the vacuum dipole model i.e. if a corotating plasma does not exist around CCOs; Magnetohydrodynamic (MHD) dipole model (Spitkovsky, 2006) require non-vanishing spin-down even in the aligned case. This idea is worth considering given that only three of the CCOs show pulsations while there is evidence for strong magnetic fields (Shakura et al., 2012).

It is remarkable that the speed of objects in the 6th group is ‘typical’. This group involves AXP/SGRs which according to the magnetar model should have large space velocities, $\sim 1000$ km s*-1* due to the “rocket effect” (Duncan & Thompson, 1992; Thompson & Duncan, 1993). This may be because it is their super-strong quadrupole fields, rather than the dipole fields, that distinguishes these objects from the rest of the isolated pulsars. The higher spin-down rates of these objects could then be attributed to additional torques, contributing to the magnetic dipole torque, due to winds (Harding et al., 1999; Thompson et al., 2000) or supernova fallback discs (Wang et al., 2006; Ertan et al., 2007).

The surface temperatures of pulsars in the 6th cluster are found to be relatively large (see Figure 4 and Table 7). This is an expected result as this cluster involves AXP/SGRs which, according to magnetar model, have decaying fields. This also is expected according to the fallback disc model as a result of accretion (Alpar, 2001; Chatterjee et al., 2000). A critical discrimination between the two models would be to compare the surface temperatures of HBPSR’s which presumably are not accreting objects, with those of AXP/SGRs. The upper limits on the temperatures of high magnetic field pulsars PSR B0154+61 (Gonzalez et al., 2004), J1814-1744 and J1847-0130 (Keane et al., 2013) imply that they are not significantly hotter than pulsars with lower magnetic fields thus raising the question why the fields of these objects do not decay and contribute to the X-ray luminosity. Yet the temperature of B1509-58 (Hu et al., 2017) lies between the temperatures of AXP/SGRs and typical RPPs, and the temperature of PSR J1119–6127 (Gonzalez et al., 2005) is among the highest measured from typical RPPs. The observed pulsed fraction of this source is also higher than those of typical RPPs (Gonzalez et al., 2005) possibly indicating the role of strong fields and unlikely to be related to the presence of a putative fallback disc. Also, the recent analysis (Çerri-Serim et al., 2019) shows that the noise strength and X-ray luminosity magnetars lacks the correlation known to exist in accreting X-ray pulsars thus suggesting that they are not accreting.

If cluster regions do indeed have anything to do with the evolutionary connections among pulsar families, the 5th and the 6th regions may be considered as two different origins for young neutron stars on the $P-\dot{P}$ diagram. A nascent neutron star may be born as a Crab-like GRP ($P_{0}\sim 10$ ms, $\dot{P}_{0}\sim 10^{-12}$ s s*-1*) or as a magnetar ($P_{0}\sim 0.5$ s, $\dot{P}_{0}\sim 10^{-8}$ s s*-1*) assuming the field is not produced by the initial dynamo effect, but flux conservation (Ferrario & Wickramasinghe, 2006) depending on its spin parameters, respectively. The classification presented here clearly relates low-B magnetars with XDINSs and RRATs and implies they could descend from magnetars by field decay. Yet statistics suggests magnetars are much less common objects compared to XDINSs and RRATs. This evolutionary scenario could, however, be validated only if magnetar stage is very short or there are many undiscovered transient magnetars. The recent identification of the high energy component in XDINS (Yoneyama et al., 2019) indeed suggests that magnetars and XDINSs are evolutionary linked.

There are, of course, strong selection effects on the representation of pulsars on the $P-\dot{P}$ diagram. CCOs, RRATs and XDINSs should be much more common in the galaxy than they are detected. Yet it is not possible to do the similar analysis with a ‘flux limited sample’ at present as such filtering of data would result with much less instances of data than the model we employed, DPGMM, requires. Given that many types of pulsars are lower-represented in our sample, it may be argued that some of the clusters discovered by DPGMM are artefacts of describing an intrinsically non-Gaussian distribution. Yet, DPGMM classification we present here broadly coincides with the existing astrophysical classification despite it is an unsupervised model.

In summary, we have shown that DPGMM, an unsupervised machine learning model, identifies 6 pulsar classes: two for recycled millisecond pulsars corresponding to different accretion histories; one roughly coinciding XDINSs, low-B magnetars and RRATs, but involves many RPPs as well; one involving the bulk of RPPs and some RRATs; one involving the very young Crab-like RPPs including the bulk of the GRPs; and finally one that involves the AXP/SGRs and HBPSRs.

Acknowledgements

KYE acknowledges support from TUBITAK, the national foundation of science in Turkey, with the project number 118F028. We thank M. Ali Alpar and Erbil Gügercinoğlu for useful discussion.