© The Authors 2024

1. Introduction

A pulsar is a rapidly rotating and highly magnetised neutron star. When there is misalignment between the spin axis and the dipolar magnetic field axis, with angle α > 0°, and given an adequate configuration with our line of sight¹, a pulsar can be detected by radio telescopes as periodic short bursts of energy. While individual pulses are highly variable for any given pulsar on the time-span of an observation, integrating over a sufficiently large number of rotations and folding temporal observations modulo the spin period P produces a pulse profile (intensity or flux density as a function of the rotation phase ϕ) that, to the zeroth order, appears to be stable over decades of observation (Johnston et al. 2021). The integrated pulse profile of radio pulsars represents one of their most informative traits and is a practical tool to gain insights about key properties determined by time-stable factors such as the geometry or dominance of the strong magnetic field (Lorimer & Kramer 2004).

Linking the pulsar rotational and energetic properties that make up its ‘identity’ (P, spin-down rate Ṗ, estimated characteristic age τ_c ≈ P(2Ṗ), surface magnetic field strength $B\propto \sqrt{P\dot{P}}$, or spin-down energy Ė/Ṗ/P³) to its geometry and exact radio emission pattern remains difficult (Karastergiou et al. 2015). The known pulsar population displays pulse profiles of diverse morphology², where each pulsar seems to have its own multi-frequency profile signature. If we exclude imprints on the profile morphology of propagation effects, such as pulse broadening impeded by the interstellar medium, profiles generally consist of a small number of Gaussian-like shaped components (Kramer et al. 1994); profiles sometimes exhibit a slimmer shape that can be fitted with Lorentz profiles, especially at higher energies (X-rays and γ-rays; e.g. Ferrigno et al. 2023).

Phenomenological studies have considered a form of classification based on the number of distinct components (single, double, triple, multiple/complex; Backer et al. 1976) and their polarimetric properties (Rankin 1983) with discriminants such as the location of the circular polarisation sign change, which is indicative of the crossing between the line of sight and the magnetic pole (e.g. Weisberg et al. 1999; Johnston & Weisberg 2006). Individual components have been studied over varying frequencies in the radio band, allowing the evolution of the beam via tracers, such as the component width, the phase separation between components, and spectral index, to be tracked.

Pulsar models have been developed to explain the source of components in relation to the magnetosphere. Given that P defines the extent of the pulsar light cylinder, the seed plasma responsible for the pulsar emission is considered to take form on the polar cap – the area at the surface of the neutron star containing open field lines. The beamed radiation is emitted by the plasma streaming along the active field lines. In the curvature radiation model (Gangadhara 2004), emission propagates tangentially to dipolar magnetic field lines. In this picture, active field lines are expected to be circularly symmetrical about the magnetic axis, resulting in a hollow cone of emission from the last open field lines with a radius that increases with height from the star surface, and with emission height decreasing closer to the magnetic axis.

When observed at frequencies below about one gigahertz, it is common to observe an evolution of pulse width as a function of frequency: Pulsar profiles tend to be wider at lower frequencies than at higher frequencies. Moreover, in double or multiple profiles, the outer components tend to show a greater outwards shift in phase with a decreasing frequency than those nearing the profile centroid (Mitra & Rankin 2002). This phenomenon is generally interpreted as a given frequency corresponding to a specific beam radius and emission height, and it is colloquially known as radius-to-frequency mapping (RFM; see e.g. Cordes 1978; Thorsett 1991; Mitra & Rankin 2002). This evolution is expressed using the full width at half maximum (W50) or at ten percent of maximum (W10). Some pulsars match different phenomenological classes at different frequencies, where one or more components appear or vanish with varying frequencies (see e.g. Johnston et al. 2008). These frequency-dependent characteristics have helped in determining where in the magnetosphere the radio emission arises.

To explain the symmetrical features observed in many pulsar profiles, the core and hollow cones model (Rankin 1983; Mitra & Rankin 2002) was proposed, comprising a core component of pencil-beam shape positioned along the magnetic axis surrounded by one or more concentric hollow cones. Alternatively, to explain complex and asymmetrical profiles, proposed models include the patchy model (Lyne & Manchester 1988), where emission occurs at random locations within the beam convolved with an annular window function, and the fan beam model (Michel 1987; Dyks et al. 2010; Wang et al. 2014), where the emission region is constrained to elongated streams that follow the paths of a group of active field lines. Recent observations of PSRs J1906+0746 and J1136+1551 (Desvignes et al. 2019; Oswald et al. 2019, respectively) seem to favour the latter model, at least for these pulsars.

The abovementioned Rankin morphology scheme comprises five broad classes: core-single (S_t) pulsars dominated by a central beam emission; core-cone triples (T) with a clear core component sandwiched by two conal outriders; pulsars in class multiple (M) with generally five components, including a central core component and distinct inner and outer conal components; conal-singles (S_d) with emission originating from a highly non-central trajectory across a hollow-conical emission beam; and conal double (D) pulsars with two prominent conal components and broad low-amplitude vestigial core emission in between that tends to become weaker at high frequency. Furthermore, Rankin has discussed a potential quadruple (Q) class where core emission has stopped in M profiles. A similar evolutionary behaviour may occur in T profiles. These classes are not fixed for a given pulsar but are rather part of a morphological continuum where a pulsar class will evolve into other classes. Specifically, there are three typical paths as a function both of radio frequency and orientation with respect to the line of sight: (1) (T?) → D → S_d; (2) S_t → T → D; and (3) S_t → T → M → (Q?). These paths give S_t (core, connected to T) and S_d (cone, connected to D) their name.

By investigating the relation between profile complexity and three physical properties (P, Ė, τ_c), (Karastergiou & Johnston 2007, hereafter KJ2007) visually classified members from an unbiased sample of 283 normal pulsar profiles (P > 20 ms) with a high S/N as belonging to one of three categories: single, double, or complex. The young, fast-spinning (P < 150 ms), highly energetic (Ė > 10³⁵ erg s⁻¹) pulsars from the sample did not display complex profiles, with a 60:40 ratio between single- and double-component profiles. Complex profiles appeared abruptly at P ∼ 150 ms, τ_c ∼ 10⁵ yr, and/or Ė ≲ 10³⁵ erg s⁻¹. Beyond this transition point, the older, slower spinning, less energetic pulsars included complex profiles, with relative fractions between single-, double-, and multiple-component profiles being roughly 45:25:30. Further simulations estimated limits on emission height range for the young and old pulsars, with the young pulsars emitting at a narrow range of high altitudes (∼950–1000 km above the star surface), and the older population emitting anywhere between 20 to 1000 km, with a variable number of active regions at a given height.

Weltevrede & Johnston (2008) investigated morphological differences between low and high Ė pulsars by counting the number of components fitted in total intensity profiles. Results hinted at subtle differences in profile complexity. However, a negligible complexity difference was identified between low and high Ė. For example, pulsars J1034−3224 and J1745−3040 were classified as the most complex by this method and were indeed visually complex; while others such as J1302−6350 were also ranked as complex while being judged as double by eye, though with highly asymmetric components, which required fitting relatively many components. Component asymmetry may hint at a more complex beam composition that is perhaps caused by discrete emission regions at various heights of KJ2007 or originates from geometrical effects such as the impact angle β (e.g. Fig. 5 of KJ2007) and related aberration and retardation of the components (Gupta & Gangadhara 2003), or it may even have an instrumental origin (e.g. an insufficient sampling rate, leading to unresolved short-duration components).

While there has been considerable effort to determine the origin of components seen in individual pulse profiles via directed methods, we note that there has been little effort to systematically compare the shape of pulsar profiles using automated methods. As the full Square Kilometer Array is expected to increase the population of currently known pulsars (∼3400) by more than ten fold, a large fraction of which should already be found during Phase I (Keane 2017; Kramer & Stappers 2015), it is timely to research and develop automated methods for knowledge extraction that require little to no human intervention.

In the following sections, we present ToPP (Topology of Pulsar Profiles), an unsupervised method employing graph theory that can automatically compare and organise a collection of pulsar profiles based on similarities in their morphology and that enables investigation of relations between pulsar profiles and their physical properties.

As a first case study, we applied ToPP to the publicly available data from the European Pulsar Network (EPN) database. From that database, multi-frequency total-intensity profiles can be reliably and relatively homogeneously extracted. Since the start of our work (see Vohl 2021a,b), graph theory has also been used on the P−Ṗ diagram (García et al. 2022; García & Torres 2023) in isolation. Here we focus on the combined set of multi-frequency profiles and the pulsar spin parameters, whose combined information remains largely unexplored. We find discrete evolutionary tracks varying from simple single-component profiles towards a diverse mixture of more complex profiles – interpreted as a mixture of pulsar class evolution and line-of-sight impact angle towards the magnetic axis. We investigate how to utilise the information encoded in the graph topology to predict the class of unknown pulsars and explore relations between morphology, graph locality, and physical parameters.

In this article, we thus limited the scope of our application to morphological comparisons of total-intensity profiles over a number of frequencies. We plan to apply the method to other large datasets in the near future, especially those that include homogenous polarisation information (e.g. Posselt et al. 2023), and to collections of fast radio bursts (Petroff et al. 2019; Cordes & Chatterjee 2019; Pleunis et al. 2021).

The next sections of this article are organised as follows: Sect. 2 describes the method and the test data. Section 3 describes the experimental results. Section 4 investigates profile classification using the MST. Finally, Sect. 5 discusses the analysis results, their potential physical implications, lessons learned, and future prospects.

2. Organising pulse profiles

In this section, we describe our methodology for data acquisition, preparation, and analysis.

2.1. Data acquisition from a semi-structured source

The EPN database is a rare public collection of curated multi-frequency pulsar profiles. The database provides access to 2458 integrated profiles for 840 pulsars, taken from 77 individual publication references³ spanning over three decades. Profiles are submitted for inclusion to the database on a voluntary basis. Each pulsar included in the database has one or more profiles taken at specific observing frequencies originating from one or more references. Observations can be visualised interactively in a web browser, and corresponding files can be downloaded.

Profiles are submitted by authors in a range of data formats, which have been curated into semi-structured files of JSON⁴ format. We collect⁵ each pulsar observation appearing on the website. Each file includes a mix of metadata, and a number of curated arrays corresponding to the phase, and total intensity (stokes I). Additionally, some observations also provide linear and circular polarisation information (stokes $L=\sqrt{Q^2+U^2}$ and stokes V). Unfortunately, the number of reasonably reliable, calibrated polarisation profiles is small, and we do not use polarisation information in the current study. For each pulsar, we also collect⁶ measured and derived properties from the ATNF pulsar catalogue (Manchester et al. 2005).

2.2. Data preparation

After collecting information for each pulsar profile from the EPN website, we link each observation to its corresponding file, and store this meta-information (pulsar name, reference, observing frequency, and file location) into tabular data⁷ to ease processing and analysis, and reduce the odds of introducing human errors.

As mentioned in Sect. 1, the pulse profile is composed of a given measured quantity (intensity or flux density) as a function of ϕ sampled at n discrete intervals ${\mathrm{\Delta }}_{\varphi }=\frac{360°}{n}$. Within the EPN sample, profiles are sampled in phase at various rates varying between 26 to 4096 bins as provided by the various authors⁸, with a majority of profiles near 512 bins⁹. We re-grid profiles to a common resolution of 512 bins by applying a spline interpolation using the zoom function from the ndimage module of the Scipy¹⁰ (Virtanen et al. 2020) python package.

We normalise stokes I intensities for all profiles to a common range of [0, 1] to allow comparison between pulsars. Furthermore, we apply a median correction by subtracting the profile median intensity. To compute the S/N of a profile array X, we evaluate the signal as the maximum value of the stokes I profile, and the noise by splitting the 512 phase bins of the stokes I profile into chunks of 8 bins and evaluating the median of medians (μ) and median of standard deviations (σ). The S/N is then computed using the following equation:

$$\begin{array}{c}\hfill S/N=\frac{max(\mathit{X})-\mu }{\sigma }.\end{array}$$(1)

Multi-frequency profiles in our subset have been pre-aligned in phase for any given pulsar in the EPN, with the exception for pulsars J1803−2137 and J1857+0943, which we align by minimising the sum of root-square differences between their respective profiles.

To investigate the profile evolution over frequency, we divided observations into a set of discrete bins of observed frequency ranges in megahertz, hereafter frequency bin, following: [400, 700), [700, 1000), [1000, 1500), and [1500, 2000). We then created a subset of pulsars that meet the following selection criteria: (1) Observations include at least one observation for each frequency bin; (2) each observation has a peak S/N greater than 20; (3) the profile at the lowest frequency bin does not display an obvious scatter broadening tail.

For pulsars where more than one observation is available for a given frequency bin, the one with the highest S/N is selected. From the available 840 pulsars, a subset of 90 meets our set of criteria and remains for further analysis. We show the period and period derivative distribution of these selected pulsars in Fig. 1 and list them in Table A.1.

Fig. 1. P−Ṗ diagram for selected pulsars. Marker colours were defined by mapping the red, green, and blue (R, G, B) channels to the period, spin-down energy, and magnetic field strength (P, Ė, B), respectively, using inverse hyperbolic sine scaling. Additionally, we highlight pulsars in binary systems, spin-powered pulsars with pulsed emission from radio to infrared or higher frequencies, association to a supernova remnant, and x-ray and γ-ray sources.

2.3. Graph

2.3.1. Definitions

A graph G = (V, E) is a structure composed of abstract mathematical objects forming the vertex set V out of which object pairs can be connected forming the edge set E (Busacker & Saaty 1965). We note the number of instances in a set S – its size – as |S|. We represent an edge e connecting two vertices u and v as e_uv = {u, v}. An edge is directed if the order by which the two vertices it connects are visited is of importance, and undirected otherwise. The density of G is defined as $D=\frac{2|E|}{|V|·(|V|-1)}$. G is said to be complete if D = 1; that is, if every vertex is directly connected by an edge to every other vertex in V. A weight w can be associated with an edge to consider a measure of distance separating its two associated vertices; in such a case, given each edge in E has an associated w, we also obtain the weight set W. We denote the exclusion of a vertex u from the set V as V \ u. The degree deg(v) of a vertex v in an unweighted graph corresponds to the number of edges incident to v, while the degree in a weighted graph can be represented as the sum of edge weights incident to v. Following the unweighted degree definition, we refer to a vertex u of deg(u) = 1 as a leaf.

A tree is a special kind of connected graph for which there exists a series of edges connecting any two vertex – a path¹¹ – where following consecutive edges in the set (e.g. [e_ab, e_bc, e_cd, …]) does not cycle (e.g. cannot be of the form [e_ab, e_bc, e_ca]). The distance between two vertices corresponds to the number of edges in the path connecting one to the other. For a weighted graph, the distance can include the sum of the visited edges’ weights. The shortest path between two vertices is the set of edges in a graph that minimises their distance. The measure of centrality of a vertex in a graph carries information about the importance of said vertex in the graph and permits identification of the most significant vertices. In particular, we defined the closeness centrality measure (Freeman 1978) as the average of the shortest path length from a given vertex to every other vertices in V. We refer to the vertex with maximal closeness centrality measure as the root. The contraction G/e of an edge e_uv is the operation E \ e that removes the edge connecting u and v, merging the two vertices into a new vertex uv in V that is connected to all original neighbours of both u and v. If u and v are not directly connected, but rather are connected via a path, G/e removes all edges forming the path, and merges the vertices it spans.

2.3.2. Application

To organise pulsars by similarity of their multi-frequency profiles, we construct a complete undirected weighted graph G = {V, E, W} using our pulsar population, where V is the vertex set of size |V|, E is the set of $|E|=(\begin{array}{c}|V|\\ 2\end{array})=\frac{|V|·(|V|-1)}{2}$ edges, and W is the set of |W|=|E| edge weights connecting every edge in E to a unique weight in W. We represent each pulsar from our population as a vertex¹². We represent the comparison between any two pulsars u, v ∈ V as an edge e with an associated weight w corresponding to their degree of similarity, forming the set e_uvw = {u, v, w}.

From G, we computed the minimum spanning tree (MST) – a subset of edges covering all vertices, without any cycle, that globally minimises the sum of the weights – by applying Kruskal’s algorithm (Kruskal 1956), reducing the set of edges to size |V|−1. The shape of the MST carries information regarding the metric used to define the weights, where an elongated shape shows sign of structured variations in the data, while a rather circular graph would hint towards noise (in relation to the measured quantity; Baron & Ménard 2021). To quantify the MST shape, we evaluate its minor and major axes, corresponding to the estimated half width $\widehat{h}$ and the longest path length l respectively.

Given the nature of a tree, the vertex in V with the smallest closeness centrality measure is the least connected vertex v₀ (a leaf), which is expected to be located at one end of the longest path. We defined the distance between any vertex q ∈ V \ v₀ and v₀ as its level L_v0(q) > 0, where L_v0(v₀) = 0. We defined the MST root as v_r. When plotting the MST, we defined the root level L_vr(v_r) = 0 and the level of any other vertex q relative to v_r.

We computed $\widehat{h}$ as the average degree for all levels L_v0(x) ∀x ∈ V. The longest path along edges forms an ordered sequence of pulsar profiles. We found this longest sequence by applying two consecutive breadth-first searches (Dijkstra 1959). First, we started from any node of the graph, and we applied a breadth-first search, which leads to a leaf u – which is one end of the longest path. Then, starting from this leaf, we proceeded with a second iteration of the breadth-first search, traversing the MST and leading to another leaf k. This ordered set of edge-connected vertices {u, v, ..., k} of size l covers the longest path.

We statistically summarised the MST’s topology by evaluating the MST elongation η as the ratio $l/(2\widehat{h})$, which we normalised by the number of vertices in the graph, ${\eta }^{\prime }=\frac{\eta }{|V|}$. If η′ nears unity, it is indicative that the evaluated similarity measure uncovers the presence of a significant trend in the data (Baron & Ménard 2021). In a case where data are categorical rather than continuous with respect to a given metric, or in the presence of multiple distinct trends, unity is not expected.

2.4. Measuring similarity between pulsars

There are many ways in which two profiles can be compared, and these different ways can lead to very different solutions. Several distance measures to compare time series have been discussed in the literature. Esling & Agon (2012) categorise distance measures for time series based on four types: shape, feature, edit, and structure. Shape-based distances compare the overall shape of a time series on its individual values directly. Feature-based distances are applied to obtain a dimensional reduction and a noise reduction. Edit-based distances compares two time series based on the minimum number of operations that are required to transform one series into the other series. Structure-based distances are obtained by comparing higher level structures obtained by modelling or compressing the time series. For example, dimensionality reduction techniques like principal component analysis (PCA) have been utilise for time series classification (Gianniotis 2017). In this work, we will not proceed with a comparative study of distance measures, reserving such an exercise for future work. Rather, we focus our attention on a specific shape-based measure as it naturally maps to the heuristic that individual observed components relate to discrete emission regions and geometry (Sect. 1). Furthermore, the shape-based measure we employ in this work is a demonstrated robust time series classification method (Bagnall et al. 2017).

Shape-based measures comprise two sub-categories: lock-step measures and elastic measures. Lock-step measures directly compare two time series following a one-to-one relation. Examples include the Minkowski distance (L_p-norm of the difference between two vectors of equal length, and where the Euclidean distance corresponds to p = 2, also known as L₂-norm) and the Pearson correlation distance. Lock-step distance measures are sensitive to scale and time shifts. Elastic measures have been developed to overcome these problems, and generally provide better classification accuracy compared to lock-step measures (Górecki 2018). An elastic measure well discussed in the literature is the Dynamic Time Warping (DTW) algorithm (Bellman & Kalaba 1959). Whereas the L₂-norm compares sample x_i with sample y_i for every i, DTW also considers neighbouring samples to allow for shifts and stretches (Niennattrakul & Ratanamahatana 2007) – which in our context allows us to abstract the effect of RFM and compare profiles’ shapes invariant of their widths. It is a similarity measure that approximately models the time-axis fluctuation (or phase-axis in our context) with a nonlinear warping function.

We briefly describe the algorithm. DTW is a pattern matching algorithm part of the ‘dynamic programming’ family of algorithms. The algorithm first builds an (n × m) local cost matrix (LCM), with element (i, j) being the L₂-norm between x_i and y_j. The second step is to construct a warping path W = w₁, w₂, …, w_K, where max(n, m)≤K ≤ m + n − 1. This path traverses the LCM under three constraints. The boundary constraint enforces that the path must begin and end on the diagonal corners of the LCM: w₁ = (1, 1) and w_K = (n, m). The continuity constraint imposes that steps in the path are taken only among adjacent elements in the matrix, including diagonal adjacent elements. Finally, the monotonicity constraint assures that subsequent steps in the path must be monotonically spaced in time. The resulting total distance for path W is the sum of individual LCM elements that the path traverses. For an overview the DTW algorithm and its many variations and constraints, we refer the reader to Giorgino (2009), which also describe the dtw package used in this work, as well as details regarding constraints (namely the step pattern, window type, and window size) used to compute the warping path described below. Figure 2 exemplifies how distance is computed by DTW, with the warping curve on the LCM (left panel) and warping tracks (right panel) computed from the query and reference profiles, exemplified by PSR B0950+08 (J0953+0755) and J1735–0724 (B1732-07) respectively.

Fig. 2. Depiction of distance measurement by the dynamic time warp algorithm. In the left panels, pulse profiles are shown as query and reference profiles corresponding to the x and y axes of the local cost matrix where the ‘warping path’ is computed. In the right panels, the same two profiles now show a subset of warping tracks highlighting matched samples along the warping path. The top row shows a comparison of a profile from PSR B0950+08 (J0953+0755) with itself leading a path of least resistance (following the heat map ‘cost valley’) corresponding to a perfect diagonal in the local cost matrix. The bottom row shows a comparison of PSR B0950+08 (J0953+0755) as a query with J1803–2137 (B1800–21) as a reference. The path of least resistance steps off the perfect diagonal, leading to a higher distance measure. The heat map also highlights the band of size K = 256 used to constrain our search.

Given that a few profiles in the set present an interpulse (a secondary pulse separated by about 180° from the main pulse) located near ϕ ∼ −0.5, the interpulse shape would be compared to the main pulse of another pulsar and confuse our results. To mitigate this problem, we applied an asymmetric step pattern and constrained the path search to a band of size K = 256 around the main diagonal (as first described by Sakoe & Chiba 1978), which corresponds to half a rotation. While the symmetric step pattern returns the average distance for both profiles, the asymmetric case returns the distance separating a query profile with respect to a reference profile. We therefore computed the asymmetric distance by varying the query and reference with either pulsars in a pair and selected the maximal distance (upper bound) of the two to avoid a reduced distance caused by partial matching. We did not normalise the distance, as all profiles have the same length n.

During our early experiment phase, we found that for some lower S/N profiles, when DTW is applied directly to pre-processed profiles (as described in Sect. 2.2), it is sensitive to intensity fluctuations (e.g. noisy profiles being considered similar). Given this observation and the fact that the EPN profiles do not include the original time series used to compute the profiles, we fit each profile using a variational Bayesian estimation of a Gaussian mixture¹³. Using the Dirichlet process as a weight concentration prior and starting with 30 mixture components and a weight concentration prior¹⁴ of 10⁴, we computed an automatic selection of active components, where the number of components to be fitted is itself modelled as a probability distribution (Ay et al. 2020). Finally, we computed the DTW distance between models. We included a symmetry constraint in the distance measure to consider profiles with similar features but flipped in phase (e.g. bright leading component versus bright following component). We considered this symmetry freedom well justified by the fact that the order in which the components appear is determined by the pulsar spin direction on the sky as seen from Earth, which is independent of the emission mechanism. We show both the pre-processed data and models in the following sections.

3. Experimental results

We evaluate similarity between pulsars in our subset population by organising pulsars by pulse shape similarity. As mentioned in Sect. 1, profiles can evolve over frequency. To take profile morphology frequency (non-)evolution into account when comparing pulsars, we construct an MST by averaging w over multiple frequency bins. We begin by discussing the longest sequence, expected to visualise the main evolutionary trend of the pulsar profile population. We then investigate the full MST.

3.1. Longest sequence

Figure 3 shows the longest sequence within the MST, with frequency bins ordered from low frequencies (higher altitude within the magnetosphere) to high frequencies (closer to the neutron star’s surface). Each row correspond to a single pulsar with its name on the left and morphological class, if available in Rankin (1990) or Mitra & Rankin (2011), as detailed in Sect. 4. Next shown is the model profile (black curve) and the pre-processed original profile data (blue dotted curve). The pulsar located at the root of the MST (J0614+2229) is marked in bold. The two branches forming the sequence display length asymmetry, with 8 neighbours on one side of the root, and 10 neighbours on the other. Profiles vary from single profiles with little frequency evolution around the root to more complex profiles at each end of the sequence. Some pulsars show constant morphology over all frequency bins, and others show components appear, fade, or blend in with other components from one frequency bin to the next.

Fig. 3. Longest sequence optimising the MST on averaging distance from all frequency bins. The sequence includes 21% of the original set, shown in Fig. 4. The symmetry constraint (Sect. 2.4) implies that neighbouring profiles in the sequence can be considered similar when inverted in phase. All sources that are mirrored with respect to the root pulsar J0614+2229 are marked with an asterisk before their name, and we display inverted profiles at all frequencies. We note that the apparent shifts in the 400−700 MHz bin between the J0437−4715 and J2145−0750 profile are not, by themselves, a cause for concern, as the profiles displayed here have no underlying intrinsic physical alignment (in contrast to, e.g., a multi-frequency plot). Regardless, the black colouring of the warp between these two profiles indicates a relatively large distance.

The vertex distribution within the tree realisation directly relates to the MST being the result of globally minimising Σw. Consecutive profiles x and y in the sequence at a same frequency bin are connected through discrete steps of their warping path marked as vertical lines connecting point x_i and point y_j, with transparency level varying with profile’s intensity for visibility, and coloured by their atomic distance, where the averaged distance w over the four bins is used to build the MST.

The longest sequence¹⁵ comprises 21% of the pulsars in the set, with a normalised elongation η′ = 0.02 obtained from l = 11.4 and $\widehat{h}=4.7$. Given multiple branches in the MST, the longest sequence has to be composed of the two longest branches, disjoint at the most central vertex v_r. The sequence indeed shows a shape evolution distributed between two leaves of the tree, namely J2145−0750 and J0528+2200, located at both ends of the sequence. The root has a single component and is classified as S_t in the literature.

Given our distance measure, the fact that most profiles are not part of the main sequence indicates a simple fact about the data: ordering the data distribution cannot simply evolve linearly, else η′→1. As η′ is low, the tree must have multiple branches, with each branch gradually evolving from single profile (near the root) towards more complex morphology types near the leaves. The assumption is that each branch must correspond to a morphology evolution track that differs from that found in other branches. We investigate these variations in Sect. 3.2.

3.2. MST

Figure 4 shows the full MST, displaying the global distribution of multi-frequency profiles, with the lowest frequency bin at the top and highest frequency bin at the bottom of each panel. Here, all 90 pulsars are organised into branches starting from the root v_r (as described in Sect. 2.3) at level 0. Each vertex is represented as a box that includes the pulse profile, showing a phase range covering W10, plus an extra 0.05 on each side¹⁶ to highlight structure in the main pulse, with the model (solid line) and the pre-processed data (dotted line). Each child is located at level L_r(parent)+1, with the DTW distance between it and its parent indicated as the colour of the connecting segment beneath it (colour bar). Pulsars that are member of the longest sequence, as presented in Fig. 3, are marked with a border. We annotate each branch starting at level 1 with a numerical value (1 to 6), noting branches 1 and 6 are actually leaves.

Fig. 4. Minimum spanning tree obtained using profiles by averaging distance at all four frequency bins. Each box shows the profiles of a given pulsar with the rotational phase showing the ten percent width (W10) plus an extra 0.05 on each side. The lowest frequency is at the top, and the highest is at the bottom of each panel. The root (most central vertex) can be seen at level 0, and a child is located at level L(parent) + 1. The weight w (distance obtained by computing the DTW) between a child and its parent is shown as the colour of the line below each child. Pulsars forming the longest sequence (Fig. 3) are highlighted with orange borders. The digital version of this figure can be zoomed and panned to inspect individual details. Pulsars can be searched by name, which combined with an appropriate zoom, can ease localisation of them in the MST.

Figure 5 presents two additional versions of the MST. In the upper panel, boxes colours are set based on the Rankin scheme with classes taken from Rankin (1990) and Mitra & Rankin (2011), respectively. Boxes without a classification in either reference are then classified given locality (more details in Sect. 4). The lower panel presents boxes coloured via a mapping of P, Ė, and B to the red, green, and blue (RGB) channels respectively. Hence, the longer the period, the redder a point is. Similarly, the higher Ė is, the greener; and finally the higher the magnetic field strength, the bluer. This mapping provides a visual cue that allows for comparison of the boxes in the MST with Fig. 1. For example, boxes in shades of pink and purple are slow pulsars with low Ė; the pinker the stronger B. Similarly, light blue and cyan indicates slightly shorter periods and high Ė; and so on. Finally, recycled pulsars are all in shades of green, given their short periods and low B.

Fig. 5. Same MST as in Fig. 4. Upper panel: Rankin classification scheme. (See Sect. 4 for details.) Grey borders correspond to pulsars without a clear prior classification in the reference papers; classified here using the method described in Sect. 4. To simplify classification interpretation, instead of using a colour bar to display DTW distance between neighbouring vertices, we display w directly below each branch. Lower panel: Relation between pulsar positions in the tree and in P−Ṗ space. The box colour results from mapping P, Ė, and B (Fig. 1) to the red, green, and blue (RGB) channels, respectively.

We can see from the MST structure that profiles over a range of frequencies comprise more complexity than in the case of considering a single frequency bin. In addition of profile morphology variation on a per frequency basis, variations of morphology (or not) across frequency bins for individual pulsar and across pulsars become discriminant. As we discussed in Sect. 3.1, each branch forms a specific morphology evolution track, starting from the simple profile v_r from which various morphological types evolve. The general trend seems to be linked to the degree of symmetry in the single ‘seed’ profile (at vertex level 1) that then morphs into other profile kinds.

The lower panel of Fig. 5 provides hints that the relation between profile morphology and spin parameters is not purely random. Indeed, pulsars with similar shades (and therefore similar properties) tend to cluster, with a number of sub-branches showing parameters evolution. For example, this is the case for the sequence between J1823–3106¹⁷, J2354+6155¹⁸, J1136+1551¹⁹, J2048-1616²⁰, and J0528+2200²¹, evolving from blue to purple to pink with monotonically increasing periods, monotonically decreasing spin-down energies, and generally increasing magnetic field strengths. Similarly, all recycled pulsars are found near or at the leaves in branch 4.

The two leaves at level 1 are occupied by PSRs J0742−2822 and J2313+4253, respectively, classified in the literature as S_t and S_d?. Both stars display profile frequency evolution with a blended trailing component with increasing amplitude towards higher frequencies. PSR J0742−2822 may therefore be a S_t/T given the circular polarisation sign shift at around 600 MHz, and J2313+4253 a conal single given its lack of sign shift in circular polarisation (e.g. Gould & Lyne 1998). This represents a case where polarisation would be crucial to properly distinguish profile class, as the shape and frequency evolution led to classifying J2313+4253 as S_t – supporting our desire to include this in future work (see Sect. 5). Out of the four remaining branches, branches 3 and 4 comprise the bulk of the pulsars in the set; branch 5 being the third in importance, and branch 2 having a single child.

Branch 2. (two vertices) is made of PSRs J1935+1616 and J2018+2839. These two apparently similar profiles according to DTW are also classified as core and conal singles respectively. As our dataset did not allow us to dependably utilise polarisation information, our algorithm is sometimes not yet capable of differentiating between these two based on shape alone. While J1935+1616 has a strong sign shift in circular polarisation over a wide range of frequencies (Gould & Lyne 1998, though displaying a second blended leading component in the 700–1000 MHz bin)²², only the circular polarisation profile of J2018+2839 at 408 MHz²³ shows a sign shift below the trailing core peak. However, this is not seen at the three highest frequency bins of our sample, which rather show an antipodal profile (bump near pulse centroid) – hence a S_d classification.

Branch 3. (42 vertices) starts with vertex pulsar J1903+0135, whose relatively simple profiles that narrow and widen as we move from high to low frequencies. At level 2, the rightmost short branch shows a similar wide to narrow to wide evolution, with the addition of outriders that are more prominent in the highest frequencies. As we move up this branch, one of the two outrider components becomes more prominent than the other, and becomes more and more resolved at higher frequencies for the leaf (J1847−0402). The level 2 leftmost branch instead shows is similar to this last branch, with the subtle difference of showing signs of more than one component within the outrider components (e.g. J2002+4050 and J0332+5434). While these are classified as T, one could equally argue that these are M.

The main central branch at level 2 further splits into four sub-branches of diverse variation patterns at level 3. In this branch, the general trend seems to be a bright narrow component of constant (relative) amplitude, with one or more components gaining amplitude towards lower frequencies. The level 3 leftmost branch morphs into M profiles, with two outer components, and a central component of lesser amplitude, with its amplitude increasing towards lower frequencies. The next branch from the left contains single profiles, including S_t and S_d. The third branch only contains T and M profiles, with one stable component and other varying over frequencies.

Branch 4. (33 vertices) begins with J0629+2415, a pulsar that shows an asymmetric profile at all frequency bins, with a component fading towards lower frequency. Pulsars following the main branch that is part of the longest sequence tend to conserve this general trend, although sometimes inverted in phase. Still in this branch, morphology begin to vary sharply at level 4 (J2326+6113). From there, sub-tracks form. The leftmost branch of level 5 onwards is inhabited by wide, multi-component profiles, with components being more or less resolved at varying frequency bins. It is interesting to note that while J0814+7429 (B0809+74) is classified as a conal single in the literature, more components are visible by eye. This agrees with the D → S_d track of B0809+74 in Rankin (1983), and much of the discussion on this pulsar since (e.g., van Leeuwen et al. 2002). Here we thus have a number of profiles where a conal single profile at one frequency evolves towards either a triple or multiple. The rightmost branch of level 5 has the pulsar J1932+2220 with one of the most peculiar morphology evolution in the set, with profiles widths varying from fairly wide at our highest frequency bin, then narrowing until the 700–1000 MHz bin, and then widening again at the lowest frequency bin. Potentially we here see a transition from intervening broadening (dispersion, scattering) at lower frequencies – J1932+2220 is relatively high DM – to intrinsic broadening at higher frequencies. This evolution pattern is only found again in J0534+2200. Other pulsars in this branch form a mix of fairly complex profiles evolution over frequency, and include all the recycled pulsars in our sample (green boxes in the lower panel of Fig. 5). Profiles show individual components amplitudes increasing towards higher frequencies. Moreover, a number of cases have components shifting in phase, alike aberration and retardation (Gupta & Gangadhara 2003).

Finally, profiles in the rightmost branch starting at level 2 have components being resolved at the lowest frequencies (top of a panel), which blend at higher frequencies, perhaps due to extra components becoming apparent and creating a plateau (as most obviously seen in J1900−2600).

Branch 5. (eight vertices) starts with J0139+5814 which is formed of single profiles at all bins, with minor signs of an extra components at higher frequencies. It also broadens towards lower frequencies. Pulsars in this branch also show an asymmetry in their main component, and generally bifurcate towards having at least two distinct components at lower frequencies, with the faintest component gaining in amplitude as we move towards higher frequencies. Here again, components tend to shift in phase as a function of frequency.

Finally, we can utilise the MST to estimate morphological complexity since the root is occupied by a simple single profile, and leaves tend to be occupied by more complex profiles. In Fig. 6, we colour each point by the DTW weight w_r, i between the root v_r and any other vertex v_i. Here, pulsars with w_r, i ≲ 0.5 are generally composed of single-component profiles. At higher w_r, i ≳ 0.5, profiles become more complex. The most simple profiles lie in the slow pulsar cloud (P > 0.1s, 10⁻¹⁶ ≲ Ṗ ≲ 10⁻¹³). More complex profiles w_r, i ≳ 1.5 are all found at extremes (low and high Ṗ, fast P). Pulsars with w_ri ≳ 2 are found at all P, but are found below Ė ∼ 10³⁵ erg s⁻¹. We show cumulative distributions of physical parameters as a function of w_r, i ranges in Fig. 7. There are hints that more complex profiles tend to lower Ė before diverging towards the millisecond period regime, which is in agreement with KJ2007. Slightly more noticeable even in Fig. 6 is that profile complexity w_r, i increases with characteristic age τ. The Pearson correlation coefficient between w_r, i and τ between is 0.34, not negligible, but low.

Fig. 6. MST (left) and P−Ṗ diagram (right) coloured with DTW distance wr, i between a given pulsar i and the root pulsar r.

Fig. 7. Cumulative distributions (left: likelihood; right: counts) of physical parameters as a function of profile complexity (wr, i, Fig. 6). All values were taken from the ATNF pulsar catalogue using parameters P0, P1, EDOT, BSURF, Dist, and ZZ, respectively.

4. Classifying new profiles from nearest neighbours

4.1. Motivation and methodology

Previously, newly discovered pulsars would generally be classified by eye into the Rankin (1983) scheme as one of core single (S_t), conal single (S_d), double (D), triple (T), or five-component/multiple (M). In this scenario, an astronomer compares the profile and behaviour of the new source to that of known sources. In this section, we investigate the case of automating – and potentially improving and unbiasing – the human-based classification through the use of the MST. In particular, each source in our sample that is not yet categorised in the literature is compared to its nearest neighbours in the tree. From that we derive its class.

To carry out the experiment, we set the following rules. For each unclassified pulsar, we: (1) Set its class given the majority class of its direct neighbours; (2) if there is a tie among neighbours, assign the class of the neighbour with the smallest w; and (3) if there remain unclassified pulsars after having visited all vertices in the MST, repeat steps 1 and 2, taking into account the newly predicted classes, until all pulsars have been classified.

As a test case, we apply this classification methodology to pulsars in our set that are unclassified in Rankin (1990), Mitra & Rankin (2011). Rankin (1990) classes are based on observations taken directly by the author at 325 MHz, 1.4 GHz, and 1.6 GHz, and by Rankin et al. (1989, 1.4 GHz with Arecibo) plus Xilouris et al. (1991, 1.7 GHz with Effelsberg). Mitra & Rankin (2011) classes are based on observations taken at 325 MHz and 1.4 GHz. As discussed earlier, we are aware that our current dataset necessarily consists of total-intensity profiles only, while the quoted studies also incorporate polarisation information from bespoke, targeted observations. We argue, however, that the current test case is the most appropriate first step for a population-based approach. For instance, having profiles at four frequency bins should provide a proxy to components’ width spectral evolution, which could in principle differentiate between S_d and S_t.

Since a given profile morphology can evolve at different frequencies for a single pulsar, it is not straightforward to simply utilise automatically the classes gathered from the literature to our test set. As we mentioned in Sect. 1, Rankin hypothesised three typical evolution paths as a function both of radio frequency (ν_lower → ν_higher) and orientation with respect to the line of sight: (1) (T?) → D → S_d; (2) S_t → T → D; and (3) S_t → T → M → (Q?). There is, however, no existing human-labelled classification of pulsars on such tracks; only in individual classes. Therefore, the classes we use for this exercise are necessarily partial and not necessarily representative of which track a profile falls into. Moreover, one can imagine that track 1 could be combined to either 2 or 3, lending only the two following tracks: (2′) S_t → T → D → S_d and (3′) S_t → T → M → (Q?), with track 2′ having a more obtuse β than that of track 3′. Given these limitations of using a fixed class for multi-frequency profiles, in addition to the classification obtained using MST shown in Fig. 5, we also computed classes on a per-frequency-bin basis²⁴. Results are presented and summarised in Table 1.

4.2. Results

We discuss our results in terms of internal consistency, by contrasting the profiles’ classification with our own qualitative classification opinions, and by comparing against the literature.

4.2.1. General classification performance

We first discuss the classes obtained on a single frequency bin, as shown in Table 1. We see that in 25 out of 28 cases, a majority of single-frequency classification agrees with one another. We consider this 89% agreement as solid evidence that the classification is consistent. In half of the cases, the majority is 3 out of 4. In 3 cases it is unanimous. It appears the algorithm is robust in its classification.

In most cases (20 out of 28), the classification obtained when using the average over all frequency bins (Sect. 3) in the construction of the MST agrees with the single-frequency MST majority. We interpret this as further affirmation of the use of our algorithm. We will discuss the outliers in more depth in the next subsection.

4.2.2. Qualitative discussion of individual complex classification cases

While the method generally agrees, three cases cannot reach a majority from their single-frequency classification, as they are tied 2–2. PSR J1832–0827 is split between S_t and M, J1901+0331 (B1859+03) is split between S_t and S_d, and J1946+1805 between S_d and T. With the exception of the classic indecisive single case that would require polarisation information, the other cases do show signs of both classes. For instance, J1832–0827 appears by eye as a partial cone (conal components only visible on one side of the core component), with a sharp narrow component that is likely the core component, and leading conal outriders. We note that the difference between T and M isn’t always obvious by eye, which may simply reflect the presence of blended components.

In a few cases, the multi-frequency classification differs from the single-frequency majority:

J0406+6138. has an M majority, and a T for the multi-frequency MST. In this case, individual classes vary between M, D and T. Therefore, it is likely that we are observing the pulsar at relatively high β angle. This pulsar has an exceptionally steep spectral index (Krishnakumar et al. 2019), which is generally associated with core emission (see, e.g., Rankin 1983).

J1722–3207. has a T majority, and is classified as S_t by the multi-frequency MST. This is likely due to the clear presence of a leading narrow component that hardly vary over frequency, while trailing and blended conal components slowly fade as we get to higher frequencies. Ilie et al. (2019) find a steep polarisation-angle (PA) swing for this source. That is generally suggestive of core emission, agreeing with the classification.

J1823–3106. has a T majority, and is classified as S_d by the multi-frequency MST. In this case, the triple class is hard to understand given the models appear as single peak profiles. A closer look at FB2 does show a leading faint conal component that likely gave a slight asymmetry to the global morphology, or simply scatter broadening. One could therefore argue that this profile is on the T → S_d track. However, the steep PA swing in Ilie et al. (2019) does not immediately strengthen that case.

J2225+6535. has a D majority, and is classified as T by the multi-frequency MST. The well-separated components do appear as double components. However, the asymmetry of the components enables one to argue that extra components are present. There are here some discrepancy between the model and the data, and it is possible that with a better model, the profile would have been classified as M. Wang & Han (2016) find a flat PA.

J2354+6155 (B2351+61). has a T majority, and is classified as D by the multi-frequency MST. Here again, the profile visually appears as a partial conal profile, with a main core component and one (or more, blended) trailing conal outrider(s). From polarisation information, Karastergiou & Johnston (2004) indeed conclude the emission is conal.

4.2.3. Core versus conal single

In a number of cases where multiple peaks are visible to the naked eye (Fig. 5), the classification we obtain tends towards one of the two single classes: wide profiles are often classified as conal single, while those with a narrow component are often classified as core single.

It is well known that witnessing a single peak provides no help in differentiating between S_t and S_d, either requiring observations over a wide frequency range to determine if a single profile developed outriders at high frequency or bifurcated at low frequency; or polarimetry to identify the characteristic circular polarisation frequently exhibited by the core species or linear polarisation associated with the conal variety (Rankin 1990). Therefore, it is expected that the classification performed in this work would confuse S_t and S_d in certain cases.

For example, J1822–2256 (B1819−22; Davies et al. 1972) is classified as core single in Table 1. In their subpulse-drift paper Basu & Mitra (2018) find the source is consistent with being S_d. Their conclusion is based mostly on polarisation data. However, the zero-crossing in Stokes V²⁵ indicates that the line of sight crosses the pole and, hence, that we likely are looking at a core component. This is not the case in profiles at all frequencies, and hence, the single flavours may vary between cone and core components at different frequencies.

It may well be that characteristics such as the pulse width or asymmetry automatically enable one to differentiate them. We note that in the multi-frequency MST, some unclassified pulsars are surrounded by classified pulsars, and their predicted classes appear to be sensible. This is the case for J1832−0827 (B1829−08) for example – predicted as S_t, and shares many similarities with J1759−2205 that is also classified as S_t. That is interesting, because from sub-pulse drift or polarisation studies, little could previously be determined about the geometry of J1832−0827 (e.g., Weltevrede et al. 2006; Johnston et al. 2023, respectively).

Lastly, a select few of unclassified pulsars were found in one sequence within the MST, all of which derive their predicted conal-single S_d class from the same root source: J0630-2834. This is series J0953+0755 (B0950+08) → J1713+0747 → J0014+4746 (B0011+47) and J0437–4715 → J2145–0750. We note that B0950+08 was already tentatively classed as conal in Lyne & Manchester (1988), which is encouraging, but the source is complicated, as discussed in Karastergiou et al. (2003) and Bilous et al. (2022), for example. Three of the four subsequent sources in the sequence are millisecond pulsars (MSPs; see next section). In those, classification is generally elusive (as discussed in e.g. Rankin et al. 2017). For J0014+4746, Rankin (2022) tentatively concluded that the emission is conal, in agreement with our automatic classification.

Overall, Table 1 displays both the successes and current limitations of our method. The algorithm is generally self-consistent, and a large majority of its findings are corroborated by studies on the individual sources at hand, as described above. The limitations primarily follow from our fixed seed classes, which are used to classify other pulsars. Arguably this is a limitation of the classification scheme in general. Ideally, our conclusions in the ‘M–F’ column should be an evolution track, not single classification. However, we do not have at this stage the tools to automatically recognise such a track; and there are basically no systematic human-labelled data on these tracks either. We will discuss these results in more detail in Sect. 5, with an outlook to further improve our method, in future work.

4.2.4. The MSPs

It is striking that ToPP clusters all MSPs together. This is best visible as the green-coloured set in the top-right part of Fig. 5 (bottom panel). This was achieved purely from the profiles, without spin information. While we cannot rule out a chance alignment, given the relatively low number of profiles, we hypothesise that the actual underlying reason is twofold, and as follows. First, the profiles are placed in a part of the graph that is generally inhabited by core and cone singles. None of the MSP profile resemble conal doubles, for example, or five-component profiles. The sampling time of our MSP profiles was high enough to rule out sampling smearing as hiding these features. Thus, the MSP profiles seem to actually be self-similar. Secondly, the MSPs are grouped at high vertex levels. While they are smoothly connected to the normal population, the distances between these profiles is generally large. Jointly those two causes group the MSPs together.

The fact that there is a smooth transition between the regular pulsars and the MSPs signifies that one physical processes is fundamental to both. Much single-pulse behaviour that was previously though to only apply to regular pulsars has since been observed in MSPs too: subpulse-drift (Edwards & Stappers 2003) and mode changing (Mahajan et al. 2018), for example. We argue that classification and interpretation schemes, too, should strive to include both recycled and non-recycled pulsars.

5. General discussion and future work

We evaluated ToPP to automatically seek for morphological types, their evolution over frequency, and their global classes. For the first time, we present here a visual summary of 90 pulsars in a condensed format. All previous publications to date, to best of our knowledge, presented pulsars individually. Beyond the effort to utilise Rankin classes and predict classes from unclassified pulsars (Sect. 4), our condensed, self-organised method makes visually intelligible a number of trends in spectral evolution that only become obvious when considering a collection of pulsars. Utilising multiple frequency naturally leads to an MST that is slightly more complex to interpret than in the single-frequency-bin case. In the multi-frequency case, while we see individual branches display similar profile types, e.g. symmetrical triple, we also find profiles being arranged by frequency evolution patterns. For example, a few pulsars show narrowing of profiles in FB2 and/or FB3 range before widening at FB1 and FB4 – the most obvious cases being J1932+2220 and J0534+2200. In the case of J1932+2220, all profiles were taken from Gould & Lyne (1998). The dispersion smearing ranges and the pulse-widths for this pulsar reported in their paper at different frequency bands are listed in Table 2:

With the exception of the profile observed at 925 MHz, the dispersive smearing could be contributing around 50% or more to the pulse width (which could be up to 100% at 1410 MHz). Therefore, the observed intriguing trend of the pulse width over frequency for this pulsar might simply be determined by the instrumental setup, and not be intrinsic to the pulsar emission regions. Similar effects might also be present in other relatively high dispersion measure pulsars with intrinsically narrow profile widths.

The trends we see in the MST can be due to multiple factors. Firstly are geometric parameters such as β or ζ. Clearly these can determine how many, if any, core and cone components cross our line of sight. Secondly, there may be intrinsic sub-families of pulsars, based, for example, on relatively unchanging characteristics such as the magnetic field strength B (magnetars – regular pulsars – MSPs). Finally, the state of the emitting system may evolve when the pulsar period decreases, and this cosmic accelerator (Ruderman & Sutherland 1975; van Leeuwen & Timokhin 2012) operates on an ever lower spin-down energy budget. In this case, the pair-formation cascade over the polar cap may be less efficient at lower Ė (Timokhin & Harding 2019), possibly somehow causing the change from core to conal dominated profiles (see Olszanski et al. 2022, for a discussion).

These three general factors are not easy to disentangle. For example, if we compare the overall profile shape of J1932+2220 to that of its direct neighbour J1807–0847 in Fig. 4, both having similar profile envelopes but slightly different periods, period derivatives, spin-down energies and magnetic field strengths. It is plausible that conal components similar to those observed in J1807–0847 are also present in J1932+2220 but blended in phase with the central component due to higher B, Ė and faster P (Fig. 8). A similar argument can be put forward for PSRs J0452−1759 and J1900−2600, which are also direct neighbours in the MST. Future similarity investigations of profile shape and physical parameters that include polarisation information and high phase resolution should allow for further investigation into these questions.

Fig. 8. Two example cases of neighbouring pulsars raising the question whether conal components are visible due to a geometrical effect or to the pulsar’s state. Spin parameters were taken from the ATNF pulsar catalogue: J1932+2220 (branch 4, level 5): P ∼ 0.14 s, Ṗ ∼ 10−13.2, Ė ∼ 1035.9 erg s−1, B ∼ 1012.5 G; J1807–0847 (branch 4, level 6): P ∼ 0.16 s, Ṗ ∼ 10−16.5, Ė ∼ 1032.4 erg s−1, B ∼ 1010.8 G; J0452–1759 (branch 4, level 3): P ∼ 0.54 s, Ṗ ∼ 10−14.2, Ė ∼ 1033.1 erg s−1, B ∼ 1012.3 G; J1900–2600 (branch 4, level 4): P ∼ 0.61 s, Ṗ ∼ 10−15.7, Ė ∼ 1031.6 erg s−1, B ∼ 1011.6 G.

In the work we present here, our distance measure does not depend on spectral index or the polarisation profile, fraction, and polarisation position angle. The strength of our approach is in the large number of profiles it compares. While the data provided by the EPN database forms the largest curated and publicly available set, the metadata of its profiles do not unambiguously state if the profile is flux and/or polarisation calibrated. The data reduction underlying the submissions is quite heterogeneous, and profiles were submitted in different formats with different header requirements. We thus concluded that only flux-normalised Stokes-I profiles are robust enough at the moment for determining the distance between these profiles. Future application of our method on sets of more homogeneously analysed pulse profile collections such as the MeerKAT Thousand Pulsar Array (Johnston et al. 2020) or Pulsar Radio Emission Statistics Survey (PRESS, PI: W. van Straten) with the ultra-wide bandwidth receiver at Parkes could benefit from also including spectral index or polarisation information. We could then also extend our method to utilise frequency-dependent profile templates (Pennucci 2019) to compare pulsars rather than a select few frequency bins.

ToPP is a method based in graph theory and optimisation methods, which are inherent concepts in artificial intelligence and machine learning. Given our specific interest in investigating the morphology of pulsar average profiles, we did not aim at curating the most generic method this present work. Rather, we tailored a method that allowed us to peer at features of interest with respect to common emission models such as the core-cone model in order to automatically organise them by morphology and spectral evolution. We curated the algorithm with domain-knowledge (choice of shape-based metric, asymmetric step pattern, and symmetry constraint) rather than making it purely agnostic.

The sequencer algorithm presented by Baron & Ménard (2021) is fully agnostic, aggregating distances measured at different scales, and evaluating multiple distance metrics (L2, Kullback-Leibler Divergence, Monge-Wasserstein or Earth mover distance) to seek for the longest manifold. The sequencer utilises the graph width and length to evaluate how much information is carried by the longest sequence. Branches outside the longest sequence are completely omitted as the algorithm specifically seek for scaling relations. ToPP extends this shortcoming through visualisation of the full MST, and the ability to investigate information from specific branches and sub-branches programmatically. For very large datasets, visualising the full tree may become impractical. To that end, edge contraction (defined in Sect. 2.3) and presented in Vohl (2021b) can provide ways to synthesise the MST to its most informative vertices.

Shan et al. (2015) considered the classification of single and double peak pulsar profiles by applying the fuzzy C-means clustering algorithm in conjunction with wavelet-based feature extraction, showing singularity-based features provided smaller classification error than shape-based extracted features. We note that our method does not specifically use extracted feature to classify profiles. Furthermore, fuzzy C-means clustering requires the number C of clusters to be specified, which is limiting in cases where C is not known a priori. However, this issue could be alleviated with a Dirichlet process approach. Future efforts to generalise ToPP should consider an evaluation of wavelet-based feature extraction as potential alternative to DPGMM prior to building the complete undirected weighted graph.

As opposed to other hierarchical clustering methods such decision trees, which require a training set, ToPP (DTW, DPGMM, MST) is inherently unsupervised. Many other unsupervised clustering methods exist – each with pros and cons with respect to various aspects, such as time complexity, memory complexity, and the ability to deal with noise – and we refer the interested reader to recent reviews (Labbé et al. 2023; Yin et al. 2024).

The various MSTs presented above highlight ToPP’s capacity to organise pulsars as a function of profile morphology (given DTW), leading to the classification results presented in Sect. 4.2. Beyond classification, one could utilise an MST to recover profiles with similar (global²⁶) morphological features as a pulsar of interest. This could be done by visiting close neighbours of the pulsar of interest, possibly limiting the search by setting a threshold on w. Future work could test such a scenario through a dedicated simulation study.

6. Conclusion

The integrated pulse profile of radio pulsars represents one of their most informative traits, carrying information on properties determined by time-stable factors such as the geometry or dominance of the strong magnetic field. Pulsar profiles have enabled theorists to model the configuration of its magnetosphere and emission mechanisms. However, linking the pulsar radio emission properties to its identity and geometry remains difficult. We have presented ToPP, a methodology to automatically compare a population of pulsar average profiles and organise them by similarity. By comparing the shape of modelled profiles, we constructed MSTs, yielding branches inhabited by pulsars displaying similar morphologies and frequency evolution and evolving from simple cases near the MST root towards more complex features near the leaves.

Employing the MST, we investigated the case of predicting the ‘Rankin class’ based on a pulsar’s location within the MST. We showed that the MST can distinguish, in many cases, between core and conal single profiles without any knowledge of spectral index or any polarisation information. We found limitations in using a fixed class when considering profiles over a range of frequencies. Future work should consider using morphological tracks such as the tree tracks introduced by Rankin (1983) as labels when classifying pulsars on multi-frequency MST. Furthermore, future application of our method on sets of more homogeneously gathered and analysed pulse profile collections would allow the polarisation information to be used within the similarity measure.

Given that the MST root vertex (in a sense, the average profile) is composed of a single profile at all of our frequency bins, with little to no frequency evolution, we utilised the root to evaluate profile complexity. After evaluating the link between profile morphological complexity and pulsar identity (P, Ṗ, B, Ė; Figs. 6 and 7), we noted that the more complex profiles tend to a lower Ė, in line with results from KJ2007, and a higher characteristic age τ. Overall, the succinct MST provides a way to identify frequency evolution trends, and these trends can be of geometric or energetic nature.

With ToPP, we have only just started to scratch the surface of the condensed information provided by graph representations. As the number of known pulsars grows into the many thousands, such structured representations are essential for both finding the population trends and identifying interesting outliers.

Acknowledgments

DV acknowledges funding from the Netherlands eScience Center grant AA-ALERT (027.015.G09). JvL and YM were supported by the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199 (‘ALERT’). JvL further acknowledges funding from Vici research project ‘ARGO’ (number 639.043.815), and from CORTEX (NWA.1160.18.316), under the research programme NWA-ORC, both financed by the Dutch Research Council (NWO). We thank Michael Keith for help with accessing the EPN database. DV thanks Anna Bilous, Jason Hessels, Simon Johnston, Aditya Parthasarathy, Andrzej Szary, and Joel Weisberg for useful discussions. In memory of François Laviolette. Software: Astropy (Astropy Collaboration 2013); astropy:2018; astropy:2022 version 5.2, dtw-python Giorgino (2009) version 1.3.1, Matplotlib (Hunter 2007) version 3.6.2, Numpy (Harris et al. 2020) version 1.24.1, Scipy (Virtanen et al. 2020) version 1.10.0, Pandas (Pandas Development Team 2020) version 1.5.2, psrqpy (Pitkin 2018) version 1.2.4, scikit-learn (Pedregosa et al. 2011) version 1.2.0.

References

Appendix A: Pulsar set

This section presents information about the pulsar set utilised in our study, including pulsar names, physical parameters, and the original reference for each profile (Table A.1). Physical parameters and uncertainties were taken from the ATNF pulsar catalogue (Manchester et al. 2005). Distances in the ATNF pulsar catalogue are based on the YMW16 electron-density model (Yao et al. 2017); with relative errors within the uncertainty limits of independently measured distances on which the model is based. While relative errors have a factor of 0.9 as an upper limit, the mean and median have factors of 0.42 and 0.05 respectively. Hence, we choose to use a fiducial factor of 0.25. In addition, we present the original phase sampling distributions in Fig. A.1.

Fig. A.1. Distribution of phase sampling rate (n) among profiles in our subset population prior to re-gridding to a common n = 512 bins through spline interpolation. Each histogram constructed using Scott’s rule. See text for details on sample selection and preparation.

Appendix B: MST variants

This section presents MST variants discussed in Sections 3 and 4. In particular, Fig. B.1 presents Fig. 4 using the full phase information. Figures B.2 to B.5 present MSTs obtained using individual frequency bins and their corresponding classification results listed in Table 1.

Fig. B.1. Minimum spanning tree from Fig. 4 (lower panel) but shown here with full phase information, ϕ ∈ [ − 0.5, 0.5], where cases with an interpulse can now be seen.

Fig. B.2. Classification using profiles in the 400–700 MHz bin.

Fig. B.3. Classification using profiles in the 700–1000 MHz bin.

Fig. B.4. Classification using profiles in the 1000–1500 MHz bin.

Fig. B.5. Classification using profiles in the 1500–2000 MHz bin.

Appendix C: Spin parameters of pulsar subsets

In this section, we investigate the difference between the two main branches of the MST presented in Fig. 4, namely branches 3 and 4.

Firstly, Fig. C.1 presents the evolution of spin parameters P, Ė and B as a function of the pulsar order within the longest sequence discussed in Section 3.1. There is a general trend for pulsars period varying from slow pulsars in branch 3 (left to the root PSR J0614+2229) towards faster pulsars in branch 4. Similarly, there appear to be a correlation between Ė and sequence order, while B appears near constant throughout, with the exclusion of the three millisecond pulsars at the right end. To evaluate the monotonicity of the relationship between each physical parameter and the sequence index, we indicate the Spearman correlation coefficient r and related p-value in each panel. We evaluate the correlation for two scenarios: normal pulsars (P > 0.1 ms) and for all pulsars in the sequence. Results show that these parameters are partially monotonically correlated.

Secondly, columns (1) and (2) of Fig. C.1 present cumulative distributions for P, Ṗ, Ė, B, and Age, where column (1) includes all its pulsars, and (2) includes only pulsars with P ≥ 0.1 ms. We perform two sample Kolmogorov-Smirnov (KS) tests between the samples of both branches. While distributions in terms of P, Ṗ, B, and age appear as significant when considering all pulsars, when p-values are approximately less than or equal to 0.01, only a marginal difference can be found in the spin period when considering normal pulsars only. Finally, column (3) compares two sub-branches of branch 3: one including primarily single-component profiles and the other having more complex profiles. Based on their respective cumulative distributions, these two sub-branches appear to be drawn from the same distribution.

Fig. C.1. Pulsar spin parameters ordered as a function of the longest sequence (Fig. 3). Pulsar names (J2000.0) are accompanied by their MST vertex level in parentheses.

Fig. C.2. Cumulative distributions of spin parameters for (1) branches 3 and 4; (2) branches 3 and 4 while excluding millisecond pulsars; and (3) branch 3, level 3, sub-branches 2 and 4. We indicate the p-values obtained by performing a two-sample Kolmogorov-Smirnov test between the samples of both branches compared in each panel. Marginal differences can be found on the spin period between branches 3 and 4 when considering normal pulsars only.

Appendix D: Distance uncertainty

As mentioned in Appendix A, distance to pulsars used in this work are taken from the ATNF pulsar catalogue. These distances are based on YMW16, and we represent distance uncertainty using a fiducial factor of 0.25. Here, we want to visualise the effect of uncertainties on the cumulative distributions presented in Fig. 7. Given distances d ∈ D, we recompute the cumulative distribution after generating 10 000 realisations randomly sampled from a normal distribution of mean d, with standard deviation d/4. Results are shown in Fig. D.1.

Fig. D.1. Cumulative distributions for distances from 10,000 realisations per measured distance from our 90 pulsars sampled from a normal distribution around the distance uncertainty.

All Tables

All Figures

Fig. 1. P−Ṗ diagram for selected pulsars. Marker colours were defined by mapping the red, green, and blue (R, G, B) channels to the period, spin-down energy, and magnetic field strength (P, Ė, B), respectively, using inverse hyperbolic sine scaling. Additionally, we highlight pulsars in binary systems, spin-powered pulsars with pulsed emission from radio to infrared or higher frequencies, association to a supernova remnant, and x-ray and γ-ray sources.

In the text

Fig. 2. Depiction of distance measurement by the dynamic time warp algorithm. In the left panels, pulse profiles are shown as query and reference profiles corresponding to the x and y axes of the local cost matrix where the ‘warping path’ is computed. In the right panels, the same two profiles now show a subset of warping tracks highlighting matched samples along the warping path. The top row shows a comparison of a profile from PSR B0950+08 (J0953+0755) with itself leading a path of least resistance (following the heat map ‘cost valley’) corresponding to a perfect diagonal in the local cost matrix. The bottom row shows a comparison of PSR B0950+08 (J0953+0755) as a query with J1803–2137 (B1800–21) as a reference. The path of least resistance steps off the perfect diagonal, leading to a higher distance measure. The heat map also highlights the band of size K = 256 used to constrain our search.

In the text

Fig. 3. Longest sequence optimising the MST on averaging distance from all frequency bins. The sequence includes 21% of the original set, shown in Fig. 4. The symmetry constraint (Sect. 2.4) implies that neighbouring profiles in the sequence can be considered similar when inverted in phase. All sources that are mirrored with respect to the root pulsar J0614+2229 are marked with an asterisk before their name, and we display inverted profiles at all frequencies. We note that the apparent shifts in the 400−700 MHz bin between the J0437−4715 and J2145−0750 profile are not, by themselves, a cause for concern, as the profiles displayed here have no underlying intrinsic physical alignment (in contrast to, e.g., a multi-frequency plot). Regardless, the black colouring of the warp between these two profiles indicates a relatively large distance.

In the text

Fig. 4. Minimum spanning tree obtained using profiles by averaging distance at all four frequency bins. Each box shows the profiles of a given pulsar with the rotational phase showing the ten percent width (W10) plus an extra 0.05 on each side. The lowest frequency is at the top, and the highest is at the bottom of each panel. The root (most central vertex) can be seen at level 0, and a child is located at level L(parent) + 1. The weight w (distance obtained by computing the DTW) between a child and its parent is shown as the colour of the line below each child. Pulsars forming the longest sequence (Fig. 3) are highlighted with orange borders. The digital version of this figure can be zoomed and panned to inspect individual details. Pulsars can be searched by name, which combined with an appropriate zoom, can ease localisation of them in the MST.

In the text

Fig. 5. Same MST as in Fig. 4. Upper panel: Rankin classification scheme. (See Sect. 4 for details.) Grey borders correspond to pulsars without a clear prior classification in the reference papers; classified here using the method described in Sect. 4. To simplify classification interpretation, instead of using a colour bar to display DTW distance between neighbouring vertices, we display w directly below each branch. Lower panel: Relation between pulsar positions in the tree and in P−Ṗ space. The box colour results from mapping P, Ė, and B (Fig. 1) to the red, green, and blue (RGB) channels, respectively.

In the text

	Fig. 6. MST (left) and P−Ṗ diagram (right) coloured with DTW distance wr, i between a given pulsar i and the root pulsar r.
In the text

	Fig. 7. Cumulative distributions (left: likelihood; right: counts) of physical parameters as a function of profile complexity (wr, i, Fig. 6). All values were taken from the ATNF pulsar catalogue using parameters P0, P1, EDOT, BSURF, Dist, and ZZ, respectively.
In the text

Fig. 8. Two example cases of neighbouring pulsars raising the question whether conal components are visible due to a geometrical effect or to the pulsar’s state. Spin parameters were taken from the ATNF pulsar catalogue: J1932+2220 (branch 4, level 5): P ∼ 0.14 s, Ṗ ∼ 10−13.2, Ė ∼ 1035.9 erg s−1, B ∼ 1012.5 G; J1807–0847 (branch 4, level 6): P ∼ 0.16 s, Ṗ ∼ 10−16.5, Ė ∼ 1032.4 erg s−1, B ∼ 1010.8 G; J0452–1759 (branch 4, level 3): P ∼ 0.54 s, Ṗ ∼ 10−14.2, Ė ∼ 1033.1 erg s−1, B ∼ 1012.3 G; J1900–2600 (branch 4, level 4): P ∼ 0.61 s, Ṗ ∼ 10−15.7, Ė ∼ 1031.6 erg s−1, B ∼ 1011.6 G.

In the text

	Fig. A.1. Distribution of phase sampling rate (n) among profiles in our subset population prior to re-gridding to a common n = 512 bins through spline interpolation. Each histogram constructed using Scott’s rule. See text for details on sample selection and preparation.
In the text

	Fig. B.1. Minimum spanning tree from Fig. 4 (lower panel) but shown here with full phase information, ϕ ∈ [ − 0.5, 0.5], where cases with an interpulse can now be seen.
In the text

	Fig. B.2. Classification using profiles in the 400–700 MHz bin.
In the text

	Fig. B.3. Classification using profiles in the 700–1000 MHz bin.
In the text

	Fig. B.4. Classification using profiles in the 1000–1500 MHz bin.
In the text

	Fig. B.5. Classification using profiles in the 1500–2000 MHz bin.
In the text

	Fig. C.1. Pulsar spin parameters ordered as a function of the longest sequence (Fig. 3). Pulsar names (J2000.0) are accompanied by their MST vertex level in parentheses.
In the text

Fig. C.2. Cumulative distributions of spin parameters for (1) branches 3 and 4; (2) branches 3 and 4 while excluding millisecond pulsars; and (3) branch 3, level 3, sub-branches 2 and 4. We indicate the p-values obtained by performing a two-sample Kolmogorov-Smirnov test between the samples of both branches compared in each panel. Marginal differences can be found on the spin period between branches 3 and 4 when considering normal pulsars only.

In the text

	Fig. D.1. Cumulative distributions for distances from 10,000 realisations per measured distance from our 90 pulsars sampled from a normal distribution around the distance uncertainty.
In the text