1. Introduction
The observed peculiar velocities of isolated radio pulsars with respect to the Galactic environment are
$\sim 200-500$
km s
$^{-1}$
(e.g. Igoshev, Reference Igoshev2020), which is much faster than
$\sim 15-40$
km s
$^{-1}$
for their high-mass progenitors (e.g. Tetzlaff et al., Reference Tetzlaff, Neuhäuser and Hohle2011; Carretero-Castrillo et al., Reference Carretero-Castrillo, Ribó and Paredes2023). This suggests that neutron stars experience a significant natal kick during their formation in supernova explosions. Unfortunately, the absence of observational techniques to determine the full three-dimensional velocities of neutron stars has impeded the firm establishment of the statistical properties of their Galactic motionFootnote
a
. Therefore, the details of the distribution of pulsar kicks remain somewhat uncertain. Thus, some researchers advocate for a broad unimodal velocity distribution (Faucher-Giguère & Kaspi, Reference Faucher-Giguère and Kaspi2006; Hobbs et al., Reference Hobbs, Lorimer, Lyne and Kramer2005; Chatterjee et al., Reference Chatterjee, Brisken and Vlemmings2009), while others identify two components with typical dispersions of
$\sim 100$
and
$\sim 500$
km s
$^{-1}$
, respectively (Cordes & Chernoff, Reference Cordes and Chernoff1998; Arzoumanian et al., Reference Arzoumanian, Chernoff and Cordes2002; Verbunt et al., Reference Verbunt, Igoshev and Cator2017). The physical origin of such bimodality remains unclear. One widely discussed suggestion is that a significant fraction of isolated neutron stars could originate from disrupted binary systems (Kuranov et al., Reference Kuranov, Popov and Postnov2009; Chmyreva et al., Reference Chmyreva, Beskin and Biryukov2010; Verbunt et al., Reference Verbunt, Igoshev and Cator2017). Indeed, a supernova explosion could affect the orbital velocities of a binary system that survives the event (Beniamini & Piran, Reference Beniamini and Piran2024). However, if the explosion is powerful enough to disrupt the system, the natal kick will overlap the orbital velocity. It has been shown that neutron stars born in disrupted systems likely have identical velocity distributions to those formed from isolated stars (Kuranov et al., Reference Kuranov, Popov and Postnov2009). Therefore, this work does not distinguish between the two origins of pulsars.
At the same time, the kick velocities of neutron stars are thought to be related to their rapid spin-up at birth, as both processes could result from the interaction between the newborn star and the supernova ejecta (e.g. Wang et al., Reference Wang, Lai and Han2007; Janka, Reference Janka2017; Janka et al., Reference Janka, Wongwathanarat and Kramer2022; Coleman & Burrows, Reference Coleman and Burrows2022). Therefore, it makes sense to expect a correlation between the pulsar spin axis and velocity directions. This correlation should exist for the peculiar velocity of objects younger than
$\sim 10$
Myr. However, it breaks down at Galactic dynamical timescales, as discussed by Mandel & Igoshev (Reference Mandel and Igoshev2023).
Studies of pulsar radio emission have shown that their apparent transverse velocities are distributed anisotropically with respect to their spin axes. Evidence of a strong correlation between spin and velocity was found (Johnston et al., Reference Johnston, Hobbs and Vigeland2005; Wang et al., Reference Wang, Lai and Han2006; Noutsos et al., Reference Noutsos, Kramer, Carr and Johnston2012, Reference Noutsos, Schnitzeler, Keane, Kramer and Johnston2013). This analysis was based on applying the Rotating Vector Model (RVM) to the pulsar’s linear polarisation signal. This model assumes that the polarisation position angle follows the projection of a rotating dipole field onto the viewing plane (Radhakrishnan & Cooke, Reference Radhakrishnan and Cooke1969; Lyutikov, Reference Lyutikov2016). However, this approach has two limitations. First, two orthogonal polarisation modes are possible for electromagnetic waves propagating through a magnetic field. These X- and O- (extraordinary and ordinary, respectively) modes may be present in observations of a single pulsar and cannot be distinguished (Stinebring et al., Reference Stinebring, Cordes, Rankin, Weisberg and Boriakoff1984; Gangadhara, Reference Gangadhara1997; Wang et al., Reference Wang, Wang and Han2014). Second, the RVM is insensitive to the spin direction of a radio pulsar; a clockwise rotator exhibits the same polarisation behaviour as a counterclockwise rotator. Together, these degeneracies prevent us from determining whether pulsar velocities are aligned with or perpendicular to their spins. Both scenarios remain possible.
Nevertheless, Yao et al. (Reference Yao, Zhu and Manchester2021) have recently presented the first clear detection of three-dimensional spin-velocity alignment in PSR J0538+2817. This is the only such evidence obtained so far. Therefore, alternative observational verification involving the statistics of a subset of objects, independent of polarisation, is necessary.
This work presents a method of testing spin-velocity alignment in isolated radio pulsars based on geometric considerations. If the spins of pulsars align with their velocities, the transverse components of the latter for objects with different magnetic angles must differ. Magnetically aligned pulsars (with small spin-magnetic axis angles) move almost along the line of sight, resulting in small, weakly scattered transverse velocities. In contrast, magnetically orthogonal pulsars move perpendicular to the line of sight, meaning their transverse velocities are close to their full 3D velocities. Consequently, one would expect different transverse velocity distributions for nearly aligned and orthogonal pulsars.
We tested this idea using existing data on the magnetic inclinations, proper motions and distances of isolated radio pulsars. Additionally, we performed a detailed population synthesis to verify that the observed difference in the velocity distribution was consistent with the properties of a realistic population, and could be detected in observations.
The paper is organised as follows. In Section 2, we compare the transverse velocities of radio pulsars with different magnetic obliquities. Section 3 describes the population synthesis of these objects, with the results presented in Section 4. Section 5 provides a brief discussion, while conclusions are collected in Section 6. Appendix A contains the table with precise details of the pulsar subset under investigation and figures showing the expanded results of the population synthesis undertaken.
2. Pulsars’ velocity distributions
2.1 Magnetically orthogonal and aligned pulsars
In order to compare the transverse velocities of magnetically aligned and orthogonal pulsars, it is necessary to create representative subsets of observed sources for each alignment type. In our previous work (Biryukov & Beskin, Reference Biryukov and Beskin2023), we compiled a catalogue of 77 isolated radio pulsars identified in the literature as having extreme magnetic inclination angles, either small (
$\alpha$
near 0) or large (
$\alpha$
near 90 degrees). The most reliable method of estimating
$\alpha$
involves fitting the rotating vector model (RVM) to the pulsar’s linear polarisation signal. However, achieving high precision in this estimation often proves challenging, primarily due to the low signal-to-noise ratio in the outer “wings” of the average pulse profile. Indeed, polarisation data carry the most information about the magnetic inclination angle in this part of the pulse (Wang et al., Reference Wang, Han and Xu2023; Johnston et al., Reference Johnston, Kramer and Karastergiou2023). Therefore, the magnetic inclination angle of a particular pulsar is usually unknown.
However, the pulsar emission properties could provide valuable information for distinguishing between weakly and strongly oblique rotators without the need for precise angle estimation. These properties include average pulse widths, distinguishable core and cone emission components, the presence of an interpulse and qualitatively different polarisation behaviour in both the main pulse and interpulse. Many authors have used these properties to classify pulsars into two categories (Lyne & Manchester, Reference Lyne and Manchester1988; Rankin, Reference Rankin1990, Reference Rankin1993b; Maciesiak et al., Reference Maciesiak, Gil and Ribeiro2011; Keith et al., Reference Keith, Johnston, Weltevrede and Kramer2010; Malov & Nikitina, Reference Malov and Nikitina2013; Johnston & Kramer, Reference Johnston and Kramer2019). The list of 77 pulsars mentioned above was compiled using these classification criteria. For the present study, we extracted a subset of this list, retaining only objects with determined full proper motions and distances. We found 13 nearly aligned and 25 almost orthogonal rotators that meet these conditions. Figure 1 plots these pulsars on the classical
$P-\dot{P}$
diagram and within galactic coordinates, alongside other single, non-recycled pulsars from the ATNF database (Manchester et al., Reference Manchester, Hobbs, Teoh and Hobbs2005)Footnote b, Footnote c. Despite the limited sample size, the selected pulsars provide a representative cross-section of the overall pulsar population. Details of the selected pulsars are presented in Table 2 in Appendix A. The parameters of these objects were also taken from the ATNF. Half the distances were estimated from dispersion measures, while parallaxes were determined independently for 19 objects.
Pulsars under consideration. Left plot: the
$P-\dot P$
pulsar diagram. Orange stars show weakly aligned rotators (with small magnetic angles). Blue circles show nearly orthogonal ones (with large magnetic angles). Grey dots represent classical (non-recycled) isolated pulsars listed in the ATNF catalogue. Right plot: The same pulsars are plotted in the Galactic sky coordinates.
2.2 Pulsars’ peculiar velocities
Consider the motion of pulsars relative to the galactic inertial reference frame (x, y, z), with its origin at the Solar System Barycentre. The x-axis is directed towards the galactic centre, the y-axis is aligned with the galactic rotation and the z-axis points towards the north galactic pole. If l and b are the pulsar’s galactic coordinates and
$\mu_l$
and
$\mu_b$
are the corresponding proper motions, with d being the distance to the pulsar, then the apparent transverse velocity vector of the pulsar is
$\begin{equation} \mathbf{v}_{\mathrm t} = kd\mu_l\left ( \begin{array}{c} -\sin l \cos b \\ \cos l \cos b \\ 0 \end{array} \right) + kd \mu_b\left ( \begin{array}{c} -\cos l \sin b \\ -\sin l \sin b \\ \cos b \end{array} \right),\end{equation}$
where
$k \approx 4.74$
km s
$^{-1}$
/(pc arcsec yr
$^{-1}$
). Physically, this velocity comprises the Sun’s Galactic motion, the circular velocity of a pulsar’s local standard of rest (LSR) within the Galaxy and the peculiar velocity
${\mathbf v}_{pec,t}$
relative to the LSR. We assume that only the latter component generally contains information about the kick that a pulsar received after a supernova explosion. Therefore, it relates to the kick-velocity correlation. Although this correlation spreads out over a pulsar’s lifetime (Mandel & Igoshev, Reference Mandel and Igoshev2023), it still applies to relatively young objects that have not yet crossed their death line (see also Section 4.1 below).
Therefore, in this work, we focus on the statistics of the projection of the pulsar’s peculiar velocity onto the star’s viewing plane.
$\begin{equation} \mathbf{v}_{\mathrm pec,t} = \mathbf{v}_t + \mathbf{v}_{\odot, t} - \mathbf{v}_{\mathrm circ,t},\end{equation}$
where
$\begin{equation} \mathbf{v}_{\odot, t} = \mathbf{v}_\odot - (\mathbf{v}_\odot \cdot \mathbf n) \mathbf n,\end{equation}$
is the projection of the Solar System’s Barycenter galactic velocity. Unit vector
$\mathbf n$
is directed towards the pulsar:
$\begin{equation} \mathbf n = \left ( \begin{array}{c} \cos l \cos b \\ \sin l \cos b \\ \sin b \end{array} \right).\end{equation}$
In our calculations, we adopt an approximate value
$\begin{equation} \mathbf{v}_\odot = (10, 235, 7.5)\,\text{km s}^{-1},\end{equation}$
which is well consistent with recent estimations (Schönrich et al., Reference Schönrich, Binney and Dehnen2010; Disberg et al., Reference Disberg, Gaspari and Levan2024). Due to the uncertainties in estimating the distance of pulsars from their dispersion measures, we do not need to know this particular velocity with high accuracy. Furthermore, we also use a simplified, fully analytic model of the three-component gravitational potential
$\phi(R, z)$
, as offered by Carlberg & Innanen (Reference Carlberg and Innanen1987) and modified by Kuijken & Gilmore (Reference Kuijken and Gilmore1989), for the circular velocity
$\mathbf{v}_{\mathrm circ,t}$
. Calculating
$\begin{equation} \mathbf{v}_\textrm{ circ}(R) = \sqrt{R \cdot \nabla \phi(R, 0)}\left ( \begin{array}{c} y\\ -x + R_\odot \\ 0 \end{array} \right ), \end{equation}$
we assume it to be aligned with the galactic plane and depending only on the cylindrical radius
$R = \sqrt{x^2 + y^2}$
and the galactocentric distance of the Sun
$R_\odot = 8.5$
kpc, which is consistent with the adopted gravitational potential.
Figure 2. (top) shows the cumulative distributions of the absolute values
$|\mathbf{v}_{\mathrm{pec},t}|$
for both types of pulsars. The difference between the two distributions is clear: orthogonal rotators (shown in blue) have systematically larger and more scattered velocities. This is precisely what one would expect in the case of spin-kick alignment. The Kolmogorov–Smirnov two-sample test rejects the hypothesis that both subsets represent the same parent distribution at a
$p = 0.016$
confidence level, and the Anderson–Darling test detects an even more significant difference:
$p = 0.007$
. However, this difference is not due to the systematics of the distance estimates, as the distributions of the distances for both groups are almost identical, as can be seen in the lower plot of Figure 2.
Top plot: Cumulative distributions of the absolute values of the transverse peculiar velocities for both types of pulsars. The solid orange line represents weakly oblique (“aligned”) rotators, while the dashed blue line represents strongly oblique (“orthogonal”) rotators. There are 13 aligned and 25 orthogonal pulsars were used, whose parameters are listed in Table 2 in Appendix A. The two distributions differ: orthogonal rotators show systematically larger and more scattered velocities. The Kolmogorov–Smirnov two-sample test rejects the hypothesis that both subsets represent the same parent distribution at a
$p_{\mathrm{KS}} = 0.016$
confidence level, and the Anderson–Darling test gives a
$p_{\mathrm{AD}} = 0.007$
. This divergence can be interpreted as an imprint of the pulsar spin-kick alignment. Magnetically aligned pulsars, in particular, tend to move along the line of sight, while orthogonal rotators move perpendicular to it. Bottom plot: Cumulative distributions of estimated distances for both types of pulsars. Their similarity suggests that the difference detected in the top plot is not due to systematic differences in distance measurement.
These distributions can be explained by assuming that aligned pulsars tend to move along the line of sight. Consequently, they exhibit low dispersion in their transverse velocities. This is consistent with the spin-kick alignment hypothesis. However, it is impossible to determine whether the kick is in the same direction as or opposite to the pulsar’s spin axis. This is the first direct, kinematic, population-based evidence of pulsar spin-kick alignment.
We performed a detailed population synthesis of radio pulsars to provide a robust theoretical argument that this difference is indeed expected and detectable. The fundamental difficulty here lies in accurately accounting for the selection bias in estimates of transverse velocity. As the quantities in our subset have been measured using different instruments and techniques, our aim is merely to demonstrate the existence of the described distribution splitting effect and estimate its amplitude.
3. Population synthesis
For reasons of simplicity, in our work, we generally reproduces the well-established algorithm initially described by Faucher-Giguère & Kaspi (Reference Faucher-Giguère and Kaspi2006) (hereafter referred to as FK06) with modifications relevant to the analysis’s aims. This section provides a detailed description of our model and computational setup.
3.1 Pulsar initial Galactic positions
The coordinates and velocities of the synthetic pulsars were considered within a Galactocentric Cartesian right-handed reference frame, with the observer shifted from the centre for
$R_{\odot} = 8.5\text{ kpc}$
and located within the equatorial plane. Following FK06, we simulated the positions of newborn pulsars along the four spiral arms established by Georgelin & Georgelin (Reference Georgelin and Georgelin1976) and quantified by Wainscoat et al. (Reference Wainscoat, Cohen, Volk, Walker and Schwartz1992). All birth locations were then smoothed relative to the centroid of the spiral arms, as described in FK06.
However, the surface density of newborn pulsars was modelled slightly differently, with
$\begin{equation} p(R) = A \left (\frac{R}{R_\odot} \right )^4 \times \exp \left( -6.8 \cdot \frac{R}{R_\odot} \right ), \end{equation}$
where R is the galactocentric cylindrical radius and
$A \approx 71.3$
is a normalisation constant (Yusifov & Küçük, Reference Yusifov and Küçük2004). This distribution describes the localisation of young OB-type stars, which are believed to be neutron star progenitors. In FK06, another distribution was adopted for the same purpose. Specifically, the probability density function of evolved pulsars, p(R), was used, which differs slightly from equation (7) and was also found by Yusifov & Küçük (Reference Yusifov and Küçük2004). However, (7) seems more physically motivated and, as demonstrated below, successfully reproduces observed pulsar statistics.
For pulsar birth places along the galactic z-axis (vertical), the double-sided exponential distribution with
$\langle z_0 \rangle = 50$
pc has been used.
3.2 Pulsar initial velocities
The initial velocities of the synthetic pulsars were calculated as the vector sum of the progenitor’s circular velocity (6) and the isotropic kick velocity
${\mathbf v}_\textrm{ kick}$
. The absolute value of the latter was based on a double-sided exponential distribution
$\begin{equation} p(v_{\mathrm 1D}) = \dfrac{1}{\langle v_{\mathrm 1D} \rangle} \exp \left (-\dfrac{|v_{\mathrm 1D}|}{\langle v_{\mathrm 1D} \rangle} \right) \end{equation}$
with
$\langle v_\textrm{ 1D} \rangle$
= 180 km s
$^{-1}$
as was found by FK06. According to this, the average full three-dimensional velocity of newborn pulsars is
$\langle v_\textrm{ kick, 3D} \rangle$
= 380 km s
$^{-1}$
. The kick velocity was considered to be independent of other pulsar parameters. However, numerical simulations of the NS formation process suggest a potential correlation between the initial spin period and the velocity of the stellar remnant (Burrows et al., Reference Burrows, Wang, Vartanyan and Coleman2024). We will study this effect in further papers.
3.3 Spin-kick alignment
Two scenarios of the spin-kick relationship were investigated: isotropic and statistically aligned.
For each synthetic pulsar, a unit vector,
$\mathbf w$
, representing the orientation of its spin axis, was generated. In the case of an isotropic kick, the vector
$\mathbf w$
was modelled independently of
$\mathbf v_\textrm{ kick}$
. In the case of an aligned kick, however, it was generated so that the angle between the kick velocity vector and the unit vector representing the orientation of the spin axis followed a zero-centred normal distribution with a standard deviation of
$\sigma_{\mathrm kick} = 15^\circ$
. This value is consistent with empirical estimations by Noutsos et al. (Reference Noutsos, Kramer, Carr and Johnston2012). In particular, this value of
$\sigma_{\mathrm kick}$
corresponds to a reasonably small value that allows the distribution of the angle between the pulsar velocity and spin axes projections onto the viewing plane to be reproduced (see also Section 4.1).
3.4 Pulsar physical parameters and their evolution
Pulsar ages
The ages of the synthetic pulsars were generated uniformly within the interval from 0 to 1.2 Gyr, corresponding to a constant galactic pulsar birthrate over the last 1.2 Gyr. The upper bound of this interval is determined by the longest pulsar lifetime within the adopted spin-down and deathline models (see below). The equation of motion for the synthetic pulsars,
$\ddot {\mathbf r} = -\nabla \phi({\mathbf r})$
, has been solved numerically until the prescribed age of the pulsar. During these calculations, we control the total energy of the pulsar (kinetic plus potential) and discard the result if it diverges by more than one per cent over the pulsar’s lifetime.
Spin-down model
Unlike FK06, we adopt a more realistic pulsar spin-down model derived from MHD and PIC simulations (Spitkovsky, Reference Spitkovsky2006; Philippov et al., Reference Philippov, Tchekhovskoy and Li2014). Thus, the evolution of the pulsar spin period, P(t), is described by the equation
$\begin{equation} P(t) \cdot \dfrac{dP(t)}{dt} = \xi B^2 \left[1 + 1.4\sin^2 \alpha(t) \right]. \end{equation}$
Here, B is a constant surface (dipolar) magnetic field;
$\alpha(t)$
is the variable angle between the pulsar magnetic and spin axes; and
$\xi = 4\pi^2 R_{\mathrm NS}^6/Ic^3 = 1.75\times 10^{-39}$
sec/Gs
$^2$
. The neutron star radius
$R_{\mathrm{NS}}$
and the moment of inertia I in this equation were taken for a star with a mass of 1.5
$M_{\odot}$
, assuming the WFF2 equation of state (Wiringa et al., Reference Wiringa, Fiks and Fabrocini1988). Since we are interested in the evolution of classical pulsars with moderate magnetic fields of the order of
$10^{12}$
G, we ignore their possible magnetic field decay. Indeed, due to the proximity of the expected Ohmic dissipation timescale to the typical pulsar lifetime, their magnetic fields do not change significantly during evolution (e.g Igoshev et al., Reference Igoshev, Popov and Hollerbach2021; Aguilera et al., Reference Aguilera, Pons and Miralles2008).
The evolution of the magnetic angle,
$\alpha(t)$
, was modelled consistently with equation (9) as
$\begin{equation} P(t)^2 \cdot \dfrac{d\alpha(t)}{dt} = - \xi B^2 \sin\alpha(t) \cdot \cos\alpha(t). \end{equation}$
The distribution of the initial pulsar periods
$P_0 = P(0)$
, was assumed to be normal with an average
$\langle P_0 \rangle$
and a standard deviation
$\sigma[P_0]$
. These two parameters were estimated during the optimisation procedure. The same is true for the distribution of magnetic fields (described via
$\langle \log B \rangle$
and
$\sigma[\log B]$
, respectively).
At the same time, initial values for the magnetic angle,
$\alpha_0 \in 0..90^{\circ}$
, were taken from the isotropic distribution so that
$\begin{equation} p(\alpha_0) = \dfrac{1}{2} \sin\alpha_0. \end{equation}$
Radio luminosity
The pulsar pseudo-luminosityFootnote d at a frequency of 1.4 GHz was modelled in a manner similar to that described by Gullón et al. (Reference Gullón, Miralles, Viganò and Pons2014):
$\begin{equation} L_{\mathrm{1.4GHz}} = L_0 \cdot 10^{L_{\mathrm{corr}}} \cdot \left( P^{-3} \dot P \right)^{1/2}, \end{equation}$
where the pulsar period and its derivative are expressed in seconds and seconds per second, respectively. According to FK06, the correction parameter
$L_{\mathrm{corr}}$
was assumed to be normally distributed with a zero mean and a standard deviation
$\sigma[L_{\mathrm{corr}}] = 0.8$
. The constant
$L_0$
(mJy kpc
$^2$
) is the final free parameter in our simulations and was assumed to remain constant for all pulsars throughout a run. Synthetic radio flux was then defined as
$\begin{equation} F_{\mathrm{1.4GHz}} = \dfrac{L_{\mathrm{1.4GHz}}}{d^2}, \end{equation}$
where d is the modelled distance to the pulsar from the Solar System Barycentre in kiloparsecs.
Death line
The pulsar death line was adopted in the form
$\begin{equation} \dot P_\textrm{ crit} = (2.82 \cdot 10^{-17}\text{ sec}^2) P^3 \end{equation}$
which was initially established by Rawley et al. (Reference Rawley, Taylor and Davis1986) from the observations of a long-period pulsar. This death line is equivalent to the equation used by FK06 (
$B/P^2 \lt 0.17 \times 10^{12}$
Gs/sec
$^2$
) if one assumes
$B = 3.2\cdot 10^{19} \sqrt{P \dot P}$
Gs – the standard estimation of the pulsar magnetic field. As a purely empirical filter, equation (14) remains close to the classical theoretical prediction
$\dot P_\textrm{ crit} \propto P^{2.25}$
by Ruderman & Sutherland (Reference Ruderman and Sutherland1975).
3.5 Pulsar observational selection
Viewing geometry
We assume that synthetic pulsars emit radiation within the two identical, symmetric pencil-shaped beams, which are directed along the magnetic axis. The radius
$\rho$
of each beam was calculated according to Rankin’s formula for the opening angle of pulsar outer conal emission
$\begin{equation} \rho = 5.7^\circ \text{ }P^{-1/2}, \end{equation}$
where P is in seconds (Rankin, Reference Rankin1993a).
The width of the unscattered pulse
$w_{10}$
(at a tenth of the maximum) has been calculated for every synthetic pulsar. Geometrically, it follows the equation
$\begin{equation} \cos \left ( \frac{w_{10}}{2} \right) = C(\alpha, \rho, \theta) = \frac{\cos\rho - \cos\alpha\cos\theta}{\sin\alpha \sin\theta}, \end{equation}$
where
$\theta \in 0..180^\circ$
is the observer’s obliquity with respect to the NS spin axis:
$\cos\theta = -(\mathbf w \cdot \mathbf n)$
. The pulse width (16) has been calculated for both pulsar beams so that
$C_\textrm{ n} = C(\alpha,\rho,\theta)$
for “north” beam and
$C_\textrm{ s} = C(180^\circ - \alpha,\rho,\theta)$
for the “south” one, respectively.
The decision as to whether the simulated pulsar is “directed” towards the observer at time t was made as follows:
-
• If
$|C_\textrm{ n}| \lt 1$
then pulsar is detectable and
$w_{10} = 2\arccos(C_\textrm{ n})$
independently on the value
$C_\textrm{ s}$
; -
• If
$|C_\textrm{ n}| \gt 1$
and
$|C_\textrm{ s}| \lt 1$
then pulsar is detectable and
$w_{10} = 2\arccos(C_\textrm{ s})$
; -
• If
$|C_\textrm{ n}| \gt 1$
and
$|C_\textrm{ s}| \gt 1$
then pulsar cannot be detected.
Detection threshold and control subset
We have almost replicated the algorithm described in FK06 in order to model telescope sensitivity. Specifically, our aim was to reproduce pulsars that could be detected by the Parkes and Swinburne Multibeam Surveys (Manchester et al., Reference Manchester, Lyne and Camilo2001; Edwards et al., Reference Edwards, Bailes, van Straten and Britton2001). Within these surveys, 1057 isolated, rotation-powered pulsars were detected at a central frequency of 1.4 GHz and found in the ATNF database. However, note that
$W_{10}$
values are known only for 381 of these pulsars, and full proper motion (
$v_t$
) is known for only 106 of them.
In order to estimate survey sensitivity, the background brightness temperature,
$T_{\mathrm{sky}}$
, is required. We calculated it using a slightly different method to that employed by FK06. Specifically, while FK06 used an electronic version of the
$T_{\mathrm{sky}}$
maps obtained by Haslam et al. (Reference Haslam, Klein and Salter1981) at a central frequency of 408 MHz, in our work we adopted the analytical approximation of the same maps found by Narayan (Reference Narayan1987) in the form
$\begin{equation} T_\textrm{ sky}(\textrm{ @408 MHz}) = 25 + \frac{275}{[1 + (l/42)^2]\cdot [1 + (b/3)^2]}\text{ K},\end{equation}$
where l and b are galactic longitude and latitude, respectively, taken in degrees. We then scaled
$T_{\mathrm{sky}}$
to 1.4 GHz according to a power law with a spectral slope of −2.8 (Lawson et al., Reference Lawson, Mayer, Osborne and Parkinson1987). The brightness temperature at the sky coordinates of a synthetic pulsar, as well as the dispersion measure, DM, towards its positionFootnote
e
, were then used to estimate interstellar distortion of the signal and the survey detection threshold,
$F_{1.4, min} \propto T_\textrm{ sky} \sqrt{DM}$
. The observed pulse width,
$W_{10}$
also took into account interstellar dispersion and telescope parameters.
4. Results
4.1 Pulsar population properties
We considered a discrete grid in the space of the free parameters of our model. There are five of these parameters: two describing the distribution of the initial period (
$\langle P_0 \rangle$
and
$\sigma[P_0]$
), two describing magnetic fields (
$\langle \log B \rangle$
and
$\sigma[\log B]$
), and one describing the luminosity-period relationship (
$L_0$
). The grid steps were set to 0.05 dex for the initial spin period and magnetic field distribution parameters, and to 0.125 dex for the luminosity normalisation constant
$L_0$
. The values corresponding to the greatest similarity between the synthetic and observed distributions are listed in Table 1. The number of potentially detectable galactic pulsars and their estimated birth rate are also provided.
These parameters are consistent with those obtained in FK06 and in a more recent analysis by Igoshev et al. (Reference Igoshev, Frantsuzova and Gourgouliatos2022). Interestingly, our model predicts short pulsar lifetimes of 2.5–3 Myr on average, although much older pulsars are undoubtedly possible. This estimation is reasonable, but could be affected by some incompleteness in our understanding of pulsar beam width, the luminosity model, telescope selection effects and the spin-down law (see the discussion below). Within this time interval, however, one could expect a pulsar’s peculiar velocity to still reflect its natal kick velocity. This explains the distribution split observed in Figure 2.
The corresponding synthetic distributions of pulsar observables are shown in Figure 6 of the Appendix A. The solid orange and dashed blue lines represent the aligned and isotropic kick scenarios, respectively. The grey-filled histograms show the observed distributions for 1057 control objects. Following Faucher-Giguère & Kaspi (Reference Faucher-Giguère and Kaspi2006), we plot the classical magnetic field estimation instead of the spin period derivative:
$\begin{equation} B_\textrm{ md} = 3.2\cdot 10^{19} \sqrt{P \dot P}\text{ Gs}, \end{equation}$
where P and
$\dot P$
are taken in seconds and sec/sec respectivelyFootnote
f
. We conclude that the obtained synthetic distributions are close to the observed ones.
Additionally, Figure 7 in Appendix A shows the distributions of further pulsar parameters that are important for analysing their kinematics and viewing geometry. These are the magnetic angle
$\alpha$
, the pulse width at 10 per cent of the maximum (
$W_{10}$
), the proper motion components
$\mu_l \cos b$
and
$\mu_b$
, the transverse velocity
$v_t$
and the distance d.
These parameters do not reproduce the observed distributions very well. However, this is not unexpected, as we ignored any selection bias that could affect estimations of pulsar proper motions, distances and, therefore, transverse velocities. In our model, synthetic pulsars have systematically larger observed velocities of approximately 150 km s
$^{-1}$
. We assume that at least two types of selection effect can account for this difference. The first is systematic errors in estimating the velocities of slow pulsars, which lead to the false detection of small
$v_t$
. The second is the small number of pulsars used to verify the kick velocity distribution (8). Thus, in FK06, only the sample of 34 pulsars observed by Brisken et al. (Reference Brisken, Benson, Goss and Thorsett2002) and Brisken et al. (Reference Brisken, Fruchter, Goss, Herrnstein and Thorsett2003) was used. Ultimately, this slight inconsistency cannot affect our results. However, we can conclude that none of these parameters are sensitive to the specific spin-kick alignment scenario.
The only parameter that exhibits such sensitivity is the angle
$\Psi$
, which is the angle between the projections of the pulsar’s spin and velocity vectors. Its properties have been investigated in detail previously (Johnston et al., Reference Johnston, Hobbs and Vigeland2005; Noutsos et al., Reference Noutsos, Kramer, Carr and Johnston2012, Reference Noutsos, Schnitzeler, Keane, Kramer and Johnston2013). It exhibits a non-uniform distribution, which is the primary statistical evidence for the pulsar spin-kick correlation. The orientation of the spin axis in the sky is associated with the position angle,
$\mathrm{PA}_{\mathrm{pol}}$
, of the linear polarisation of the pulsar’s emission. In particular,
$\mathrm{PA}_{\mathrm{pol}}$
is the position angle at the maximum derivative of the PA swing. However, observations cannot distinguish between clockwise and counterclockwise rotating pulsars. This leads to an uncertainty of 90 degrees in
$\Psi$
. Additionally, we cannot differentiate between the X- and O-modes of pulsar polarisation, which exacerbates this uncertainty. As a result, the observable
$\Psi$
is the smallest angle between the direction of the proper motion
$\mathrm{PA}_{\mathrm{v}}$
and one of the two axes, one of which is given by the position angle
$\mathrm{PA}_{\mathrm{pol}}$
and the other by
$\mathrm{PA}_{\mathrm{pol}} + 90^\circ$
Population synthesis best parameters
Figure 3 shows synthetic distributions of
$\Psi$
obtained in population synthesis and the observed distribution from Noutsos et al. (Reference Noutsos, Kramer, Carr and Johnston2012). We conclude that an isotropic kick cannot reproduce the observations, whereas a spin-aligned kick can. This result is consistent with that obtained by Noutsos et al. (Reference Noutsos, Kramer, Carr and Johnston2012), but is based on an extensive, detailed population synthesis.
Distributions of the angle
$\Psi = $
PA
$_{\mathrm{v}}$
- PA
$_{\mathrm{pol}}$
, representing the divergence between the proper motion vector and the projection of the pulsar spin axis onto the viewing plane. The grey bars show the observed distribution, taken from Noutsos et al. (Reference Noutsos, Kramer, Carr and Johnston2012), based on 54 pulsars. Values are normalised from 0 to 45 degrees due to degeneracy of the Rotating Vector Model regarding the
$X-$
and
$O-$
polarisation modes in pulsar emission. Angle
$\Psi$
is the only parameter sensitive to the spin-velocity correlation. An isotropic kick scenario (shown by the blue dashed line) results in a nearly uniform distribution of
$\Psi$
, which is inconsistent with observations. Conversely, a spin-aligned kick scenario (the solid orange line) produces a more suitable, non-uniform distribution.
4.2 Velocities of aligned and orthogonal rotators
However, the main parameter of interest in the current work is the peculiar transverse velocity of the pulsars. Figure 4 shows the modelled distributions of this quantity, which are similar to those in Figure 2. The solid black lines show the resultant distributions involving every synthetic pulsar, which look almost identical in both spin-kick relationship scenarios. However, when weakly and strongly oblique pulsars are considered separately, a high spin-kick alignment produces a clear distribution split. This can be seen in the right-hand plot in Figure 4. A couple of distributions corresponding to nearly aligned synthetic pulsars with
$\alpha$
less than 25 or 10 degrees are shown by thin orange lines. Similar distributions for almost orthogonal objects (with
$\alpha$
greater than 65 and 80
$^\circ$
) are shown by thin blue lines. These two sets of distributions are separated from each other, as well as from the joint distribution. In particular, aligned pulsars show systematically smaller and less scattered transverse velocities. This is precisely expected in the case of strong spin-kick alignment, which is observed in real pulsars. The real distributions from Figure 2 are also shown on these plots by light grey lines.
These plots illustrate the paper’s primary theoretical outcome. The distributions of the modelled (synthetic) peculiar velocities,
$v_{\mathrm{pec}} = v - v_{\mathrm{LSR}}$
, are shown after projection onto the viewing plane. The left plot contains results for an isotropic kick velocity, and the right plot contains results for a spin-aligned velocity. The solid black lines show the overall distribution of pulsar velocities. However, observational selection effects are taken into account in the same way as for the distributions shown in Figures 6 and 7. The thin orange and blue lines on both plots represent the velocity distributions of pulsars that are nearly aligned (with
$\alpha$
< 25 and 10 degrees) and nearly orthogonal (with
$\alpha$
> 65 and 80 degrees), respectively. A clear effect seen in the simulated data is the splitting of the velocity distribution for these two types of pulsar in the case of spin-kick alignment. Specifically, weakly oblique pulsars have systematically smaller and less dispersed observed transverse velocities than orthogonal ones. At the same time, isotropy of the kick destroys this effect. The observed distributions from Figure 2 are also shown as light grey lines on the plots. Notably, the observed and synthetic distributions are very close to each other in the case of spin-kick alignment, even though this was not the aim of the calculations. In general, this plot shows that spin-kick alignment manifests as different transverse velocities for aligned and orthogonal pulsars in a synthetic pulsar galaxy with realistic underlying properties. Furthermore, the strength of this difference is very close to that observed.
It is notable that the modelled distributions on the right-hand plot agree well with the observed ones. While this coincidence was not the objective of the simulations, it likely reflects the adequacy and realism of the constructed population synthesis model.
Therefore, we conclude that the population synthesis performed supports the idea that the observed difference in the velocities of aligned and orthogonal pulsars is due to spin-kick correlation. Furthermore, kick velocities tend to align with the spin axis’s direction.
5. Discussion
5.1 Realistic rotational evolution of radiopulsars
Although the equations of pulsar spin-down (9) and (10) are confidently constrained from an accurate numerical study, they still remain purely theoretical results. The long-term rotational evolution of real pulsars appears to be more complexFootnote
g
. Thus, the so-called braking indices
$n_{\mathrm{br}} = 2 - P\ddot P/\dot P^2$
, are expected to lie within the range
$3..3.25$
during a pulsar lifetime, assuming a constant magnetic field (Philippov et al., Reference Philippov, Tchekhovskoy and Li2014). However, for actual classic pulsars, the values of
$n_{\mathrm{br}}$
spread in
$\sim -10^5..10^5$
, being even negative in half of the cases (Hobbs et al., Reference Hobbs, Lyne, Kramer, Martin and Jordan2004; Biryukov et al., Reference Biryukov, Beskin, Karpov and Chmyreva2007). This challenges the adopted spin-down model. Assuming that the pulsar evolves with a huge, constant braking index, one gets that it have to either stop spinning down on timescales
$\sim \tau_{\mathrm{ch}}/n_{\mathrm{br}}$
for positive
$n_{\mathrm{br}}$
, or completely lose its rotational energy (
$P = 0$
) after
$\sim \tau_{ch}$
for
$n_{\mathrm{br}} \lt 0$
. This is the problem of anomalous braking indices. Many solutions have been proposed for it so far, but none have been widely accepted.
On the other hand, pulsar population studies show that smooth spin-down laws (9)-(10) really help to satisfactorily reproduce pulsar statistics. The population synthesis undertaken in this study provides further evidence for this. This contradiction can be resolved by assuming that braking indices oscillate with a large amplitude on timescales much shorter than the pulsar lifetime. Earlier, we have estimated this timescale as
$\sim 10^3..10^4$
years (Biryukov et al., Reference Biryukov, Beskin and Karpov2012). The effect of such oscillations is averaged out during long-term evolution, which the current work is focused on. Therefore, we believe that the model (9)-(10) is still reliable and produces physically reasonable results.
5.2 Orthogonal natal kick
In our analysis, we have focused on a spin-aligned natal kick of pulsars. An alternative scenario – when kick velocity is nearly orthogonal to the spin axis – was not considered in detail. Nevertheless, we have performed simple Monte Carlo calculations for this case to ensure completeness. Thus, we treated neutron stars as rotators with magnetic angles that are strictly equal to 0 or 90 degrees. In other words, two extreme cases were considered in terms of the magnetic alignment. Of course, magnetic angle
$\alpha = 0$
suppresses any pulsar effect. But, in our calculations, we are interested in statistics of the spin-velocity angles, and want to probe the strongest possible effect. Thus, adopting
$\alpha = 0$
is still makes sense in this particular case.
To generate the new synthetic subset of neutron stars, we modelled the velocity vector of each star in two ways. The first method involved setting the angle to the spin axis as normally distributed with a dispersion of 15 degrees. The second method is the same, but for the angle between the velocity vector and the star’s equatorial plane. The azimuthal orientation of the velocity was chosen uniformly in the 0.360 degrees range. The absolute values of the velocity were taken from a Maxwellian distribution with an average of 250 km s
$^{-1}$
. By modelling
$10^7$
such synthetic pulsars, we calculated the theoretical distributions of transverse velocities, i.e. the projections onto a common plane. These distributions are shown in Figure 5. The case of a spin-aligned kick is shown on the top plot of this figure. The top plot of this figure shows the case of a spin-aligned kick. The difference in the transverse velocities of strictly aligned and orthogonal pulsars is exactly the same as that observed (see Figure 2) and obtained in the population synthesis (see Figure 4).
These are the results of a simple Monte Carlo calculation of the distribution of transverse velocities of pulsars. The top plot shows the distribution of spin-aligned pulsar velocities for magnetically aligned and orthogonal pulsars. The difference between these distributions is consistent with that observed in population synthesis. In contrast, the bottom plot shows a similar distribution for a kick perpendicular to the spin axis. The split between the two distributions still occurs in this case, but the sign is different. Magnetically orthogonal pulsars, in particular, have qualitatively smaller and less dispersed transverse velocities. This supports the conclusion that Figure 2 reflects spin-kick alignment, rather than orthogonality. See the text for details.
Conversely, a spin-orthogonal kick produces the inverse difference in these distributions. In particular, magnetically orthogonal pulsars show smaller and less dispersed velocities than magnetically aligned ones. Therefore, we conclude that the observed diversity in the velocities of the two types of pulsar indicates alignment, rather than orthogonality of their spin and 3D velocity direction. This is the first time that such evidence has been obtained for a subset of pulsars, but not for individual ones, independently of their polarisation.
6. Conclusions
This work has shown the power of the pulsar velocities analysis method in revealing the spin-kick alignment phenomenon. This effect was initially established after analysing 54 pulsars (Noutsos et al., Reference Noutsos, Kramer, Carr and Johnston2012) more than a decade ago. However, recent massive estimations of pulsar polarisation parameters, including magnetic angles for more than 400 southern (Johnston et al., Reference Johnston, Kramer and Karastergiou2023) and almost 200 northern (Wang et al., Reference Wang, Han and Xu2023) pulsars, open a new possibility to refine underlying distributions and provide a more detailed understanding of this phenomenon. At the same time, it is evident that transverse velocities of relatively young orthogonal pulsars are the best estimations of their actual 3D velocities, which can be constrained for the first time. This fact could also be used for further analysis, which will be presented in future works.
Acknowledgement
The authors are thankful to Sergey Popov and Andrei Igoshev for their useful suggestions.
Funding statement
The work was carried out within the framework of the state assignment of the Special Astrophysical Observatory of the Russian Academy of Sciences, approved by the Ministry of Science and Higher Education of the Russian Federation
Competing interests
None.
Data availability statement
All the data related to the described analysis can be provided upon request.
Appendix A. Analyzed pulsar subset details and expanded results of the population synthesis
Pulsars under consideration. These classical isolated pulsars are assumed to be either weakly or strongly inclined in terms of their magnetic angle. At the same time, they all have distance estimates and full proper motions. The ages here are characteristic ages
$P/2\dot P$
in
$10^6$
years, while the surface fields are magnetodipolar estimates of
$3.2\times 10^{19} \sqrt{P \dot P}$
in
$10^{12}$
Gs. The last column lists the references where a particular pulsar has been found to be either aligned or orthogonal.
aDM is for dispersion measure-based distances according to Galactic free electrons distribution by [REF]; PX for the parallax-based distances.
bReferences that support the type of the obliquity: [1] Lyne & Manchester (Reference Lyne and Manchester1988); [2] Rankin (Reference Rankin1990); [3] Rankin (Reference Rankin1993b); [4] Maciesiak et al. (Reference Maciesiak, Gil and Ribeiro2011); [5] Keith et al. (Reference Keith, Johnston, Weltevrede and Kramer2010); [6] Malov & Nikitina (Reference Malov and Nikitina2013); [7] Johnston & Kramer (Reference Johnston and Kramer2019).
Pulsar population synthesis results show distributions of six observables used to find optimal model parameters: pulsar spin period P, magnetic field
$B_{\mathrm md} \propto \sqrt{P \dot{P}}$
, dispersion measure DM, galactic coordinates (l, b), and flux F at 1.4 GHz. The solid orange line shows synthetic distributions assuming a kick velocity aligned with the spin axis, while the dashed blue line represents an isotropic kick. Both models use the same initial parameters from Table 1. Observed distributions for 1057 single classical pulsars from Parkes and Swinburne surveys are shown as grey bars. Population synthesis qualitatively reproduces the observed statistics well, and the plotted quantities are insensitive to the possible spin-kick correlation.
Further results of population synthesis. This figure is organised similarly to the previous one: the grey-shaded bars represent the observed distributions (where applicable), while the orange solid and blue dashed lines are for synthetic distributions. (a) Magnetic angle
$\alpha$
between the spin and magnetic angle of synthetic pulsars. Surprisingly, this parameter is insensitive to the specific spin-kick orientation scenario. This is probably due to the nearly anisotropic distribution of pulsars relative to the observer. (b) Pulsar pulse widths at 10% of maximum in miliseconds. The observed distribution is shown for 381 pulsars with measured
$W_{10}$
from the control subset of 1057 objects. (c) Distances to pulsars relative to the Solar System barycenter. Observed values are mostly based on the dispersion measure. Here, all 1057 control pulsars are shown. (d) Proper motion along the galactic longitude. The observed distribution is shown for 106 pulsars with known
$\mu_l$
from the initial control subset. (e) Similar, but for galactic latitudinal direction. (f) Transverse velocities relative to the Solar system barycenter. Relative to this particular control subset of pulsars, the synthetic population shows an excess in nearby and faster pulsars. This could be a result of a combination of selection effects in the estimation of pulsar proper motions, systematic errors in estimations of distances, as well as some incompleteness of the kick distribution model.
Open access
