Constraints on Weak Supernova Kicks from Observed Pulsar Velocities

The observed single pulsar population is characterized by typical speeds of hundreds of $\rm{km}\,\rm{s}^{-1}$, generally attributed to large natal kicks due to asymmetric mass ejection during supernovae (Hobbs et al., 2005; Burrows, 2013; Verbunt et al., 2017; Deller et al., 2019). Such large kicks are consistent with simulations of core-collapse supernovae (Wongwathanarat et al., 2012; Müller, 2020). On the other hand, the existence of pulsars in globular clusters, where escape velocities can be ${\lesssim}50\,\mathrm{\rm{km}\,\rm{s}^{-1}}$, points to the need for a subpopulation of pulsars with low kicks (e.g., Sigurdsson, 2003). Meanwhile, natal kicks that disrupt the binary inhibit the formation of double neutron stars (e.g., Beniamini & Piran, 2016; Vigna-Gómez et al., 2018). Therefore, at least some fraction of neutron stars must receive low natal kicks.

Theoretically, low kicks (${\lesssim}30\,\mathrm{\rm{km}\,\rm{s}^{-1}}$) are often associated to electron-capture supernovae and ultra-stripped supernovae. USSNe only occur as the second SN in very tight binaries with minimal mass loss and almost never disrupt the binary; thus, we do not consider them further as they do not contribute to the velocity distribution of isolated pulsars (Tauris et al., 2015). ECSNe are thought to arise from a subset of the super-Asymptotic Giant Branch (sAGB) stars, which span a zero-age main sequence (ZAMS) mass range of ${\approx}6.5-12\,\mathrm{\rm{M}_{\odot}}$, though the precise ZAMS mass range for ECSN progenitors is highly uncertain (see Doherty et al. 2017 for a review).

In this paper, we use the COMPAS rapid binary population synthesis code (Stevenson et al., 2017; Vigna-Gómez et al., 2018) in order to reconcile the paucity of observed low-velocity pulsars with the need for a population of NSs with low natal kicks. We find that by restricting low natal kicks to occur only in stars which have previously transferred mass to a binary companion, our model reproduces the observed fraction of low-velocity isolated pulsars, without inhibiting the DNS yield. We interpret this as an upper limit on the ZAMS mass range for effectively single progenitors of ECSNe, while leaving the channel for stripped stars unaffected.

Our proposal builds on indications that ECSNe from single stars are less common than initially thought (Miyaji et al., 1980; Nomoto, 1984; Poelarends et al., 2008), but could be enhanced in binaries because envelope stripping by Roche-lobe overflow onto a companion suppresses second dredge-up (Podsiadlowski et al., 2004a; van den Heuvel, 2010; Ibeling & Heger, 2013; Dall’Osso et al., 2014; Poelarends et al., 2017).

In Sec. 2, we discuss the data set of pulsar velocities observed with very long baseline interferometry and the associated selection effects. In Sec. 3, we describe our population synthesis prescriptions. We present the results in Sec. 4. In Sec. 5, we discuss the caveats of our analysis and its implications for stellar evolution modellers.

The data used in this study are astrometric measurements of isolated pulsars obtained with very long baseline interferometry (VLBI). We use exclusively VLBI measurements rather than the larger data sets based on pulsar timing or dispersion measure because VLBI provides precise measurements and we aim to avoid concerns about systematic uncertainty in other pulsar velocity measurements (Deller et al., 2019).

We select a total of 81 pulsars from a variety of studies, primarily Deller et al. (2019) but also Bailes et al. (1990); Fomalont et al. (1999); Chatterjee et al. (2001); Brisken et al. (2002, 2003); Dodson et al. (2003); Chatterjee et al. (2004, 2009); Deller et al. (2009) and Kirsten et al. (2015). We remove millisecond pulsars, pulsars known to be in binaries, and pulsars in globular clusters in order to obtain a relatively homogeneous data set of pulsars whose velocities are primarily set by SN natal kicks (and any previous isolated binary evolution) rather than dynamics. Some of the apparently single pulsars may still be in binaries, as binaries with orbital periods much longer than the observational baseline are challenging to detect.

The data for each pulsar contain sets of bootstrapped fits for the parallax, position, and proper motion. We treat the fits as posterior samples. We re-weigh them by applying a prior on the parallax to be above 0.05 mas, i.e., the distance to be below 20 kpc, consistent with Galactic pulsars, as well as a prior on the transverse velocities to be less than 2000 $\rm{km}\,\rm{s}^{-1}$ , to exclude anomalously large values.

We extract an intrinsic transverse velocity relative to the pulsar’s local standard of rest (LSR) by correcting for Galactic rotation and the motion of the Sun. We assume a flat Galactic rotation curve, with constant speed 230 $\rm{km}\,\rm{s}^{-1}$ (Bhattacharjee et al., 2014), and take solar velocity (U, V, W)_⊙ = (11.1, 244, 7.25) $\rm{km}\,\rm{s}^{-1}$ (Schönrich et al., 2010). This correction assumes that the LSR velocity vector of the pulsar today is similar to its LSR at birth, which relies on the assumption of insignificant acceleration in the Galactic potential. This is valid for very young or slow-moving pulsars, but is not justified, e.g., for a 20 Myr pulsar moving at 500 $\rm{km}\,\rm{s}^{-1}$ , which would have travelled for ${\sim}10\,\mathrm{kpc}$ in the absence of external forces. However, only pulsars with high velocities (comparable to or larger than the local Galactic rotational velocity) will be impacted, so this is not expected to affect our analysis of low-velocity pulsars. A bigger concern is that the data were not collected in a manner intended to be a complete survey. Indeed, Deller et al. (2019) note that their sample focused specifically on systems at high Galactic latitudes, which could lead to a paucity of low-velocity pulsars. However, the distributions of characteristic ages and heights off the Galactic plane are indistinguishable between our sample and the broader sample of all isolated, non-millisecond pulsars in the ATNF pulsar catalogue¹¹1 https://www.atnf.csiro.au/research/pulsar/psrcat/.

It is also conceivable that pulsars with low velocities are naturally fainter and harder to detect because of some correlation between their formation physics and an observational characteristic relevant to radio detectability that is not yet understood (Sigurdsson, 2003). We test for this possible bias by examining the correlation between pulsar transverse velocities and distances. A Malmquist bias favours detecting only radio-bright pulsars at greater distances. Therefore, if low-velocity pulsars tend to be fainter in the radio band, we should see a positive correlation between distance and transverse velocity. While there is a mild correlation between distance and velocity, we confirmed that its magnitude is fully consistent with arising from the construction of transverse velocity as a product of proper motion and distance, with significant uncertainties in the distance (typically of order 20%). We therefore see no evidence for a selection effect against low-velocity pulsars once this correlation is accounted for. Further, we find that the correlation between the transverse velocities and heights of the pulsars out of the Galactic plane is no greater than the correlation between the transverse velocities and pulsar distances, which further supports the absence of a significant bias due to preferentially observing pulsars at high Galactic latitudes. Fig. 1 shows the transverse velocity distributions at the low end of our catalog.

We use the COMPAS rapid binary population synthesis code (Team COMPAS: Riley et al., 2021) to generate synthetic pulsar transverse velocity distributions for several different SN kick prescriptions. We simulate a total of $10^{6}$ binaries per prescription, with primary ZAMS mass drawn from the Kroupa (2001) initial mass function between 5 and 150 $\rm{M}_{\odot}$ , mass ratio drawn uniformly from 0.01 to 1, and semi-major axis distributed uniformly in the log between 0.01 and 1000 AU (Öpik, 1924). Here, very wide, non-interacting binaries represent the single star population. All binaries are initially circularized ($e_{0}=0$) and use solar metallicity $\rm{Z}_{\odot}=0.0142$ (Asplund et al., 2009). Unless otherwise specified, we follow the default COMPAS prescriptions (Stevenson et al., 2017; Vigna-Gómez et al., 2018; Vinciguerra et al., 2020; Team COMPAS: Riley et al., 2021).

We assume that a star with a helium core mass above 2.25 $\rm{M}_{\odot}$ at the base of the asymptotic giant branch will undergo CCSN once the mass of its carbon-oxygen core reaches the threshold set in Eq. (75) of Hurley et al. (2000), in which we replace the Chandrasekhar mass with 1.38 $\rm{M}_{\odot}$ (Belczynski et al., 2008). We set the CCSN remnant mass according to the Fryer et al. (2012) delayed prescription. If the helium core mass is between 1.6 and $2.25\,\mathrm{\rm{M}_{\odot}}$ at the base of the asymptotic giant branch and the carbon-oxygen core mass subsequently reaches $1.38\,\mathrm{\rm{M}_{\odot}}$, the star is assumed to form a NS with mass $1.26\,\mathrm{\rm{M}_{\odot}}$ in an ECSN.

We consider a variety of different kick prescriptions, all based on the following “fiducial” prescription. Natal kicks for pulsars formed in CCSNe (hereafter, CC-pulsars, and similarly for EC-pulsars) are drawn directly from the observed pulsar velocity distribution. To obtain a 3D kick magnitude for the CC-pulsars, we randomly draw a velocity sample from the union of the transverse velocity posterior distributions of pulsars in our catalogue, selecting only among values of at least 50 $\rm{km}\,\rm{s}^{-1}$ , and divide by a random projection coefficient under the assumption of an isotropic viewing angle.

This is in contrast to the COMPAS default CCSN kick prescription (which uses the the Hobbs et al. 2005 Maxwellian model) and other evolution codes which draw from parametrized, analytical distributions, because these distributions differ from the observations at high velocities. However, this difference is due directly to the assumed shape of the natal kick distribution. The bulk of the isolated pulsar population have very high speeds ${\gtrsim}300\,\mathrm{\rm{km}\,\rm{s}^{-1}}$, much larger than typical pre-SN orbital velocities,

where $M$ is the total binary mass and $a$ is the separation of a circular binary. Consequently, the pulsar is nearly always unbound, with a final speed dominated by the magnitude of the natal kick, with some scatter on the order of $v_{\text{orb}}$. By drawing CCSN kicks directly from the observations, we are able to focus on the poorly-understood low-velocity regime.

In the Fiducial prescription (henceforth abbreviated as Fid), all EC-pulsars have natal kick magnitudes drawn from a Maxwellian with $\sigma_{\text{1D rms}}=30\,\mathrm{\rm{km}\,\rm{s}^{-1}}$.

All SN natal kicks are applied isotropically in the progenitor rest frame (cf. Bray & Eldridge 2018 and Giacobbo & Mapelli 2020, who investigate anisotropic kicks). If the SN disrupts the binary through a combination of symmetric mass loss (the Blaauw kick; Blaauw, 1961) and the natal kick, the speed of the NS is given by its asymptotic speed once it has escaped the gravitational potential of the companion. If the binary remains intact, the speed of the NS is given by the sum in quadrature of its speed in the center-of-mass frame and the binary velocity relative to the LSR, if any. We then project the NS speed onto the plane of the sky assuming an isotropically distributed viewing angle to obtain a transverse velocity prediction. NSs in wide binaries with orbital periods above 10 years are categorized together with isolated NSs as “apparently isolated”; varying the period threshold between 10 and 100 years has a negligible impact. Wide binary pulsars comprise ${\sim}20\%$ of the total apparently isolated NS population, immediately following the SN, in all our prescriptions. This fraction is a few times higher than the estimated fraction of pulsars in wide binaries obtained by Antoniadis (2021) ($\lesssim 10\%$), who considered observability through Gaia and radio pulsar surveys under a very simple model of binary properties (see also Igoshev & Perets 2019). This discrepancy is in part due to the fact that our 20% estimate does not account for the amount of time that a pulsar is observable in a wide binary, which could be disrupted by a second supernova. It would be interesting to compare the properties of wide binaries predicted by our models against observations through a similar modelling of selection effects.

ECSNe are the primary source of low-velocity pulsars in our models, so we consider variations on the Fid prescription described above that impact ECSNe. In the first variation, Kick: 1 $\rm{km}\,\rm{s}^{-1}$ (K1), we reduce the natal kicks of ECSNe down to a Maxwellian with $\sigma_{\text{1D rms}}=1\,\mathrm{\rm{km}\,\rm{s}^{-1}}$, motivated in part by hydrodynamical simulations (Gessner & Janka, 2018). The second variation, No Wide ECSNe (NW), flags and ignores ECSNe in non-interacting binaries, allowing for ECSNe only in stars which previously lost mass through Roche-lobe overflow. The final variation, No Wide ECSNe, Kick: 1 $\rm{km}\,\rm{s}^{-1}$ (NW-K1), combines the first and second variations. A list of variations is given in Table 4.

In Fig. 2, we compare the transverse velocity cumulative distribution functions of apparently isolated NSs from the Fid (a) and NW (b) prescriptions to those of the observed pulsars. Each observed transverse velocity CDF (black) is constructed by randomly sampling one posterior transverse velocity sample per pulsar. The spread in these CDFs indicates the uncertainty in the velocity measurements. Each synthetic CDF (colored) is constructed by randomly drawing as many velocities as the total number of observed pulsars $N_{P}{=}81$ from the modelled population. Their spread indicates the impact of small-number statistics. The CDFs match at velocities ${\gtrsim}400\,\mathrm{\rm{km}\,\rm{s}^{-1}}$ (both panels), validating the use of the CCSN prescription used in this study.

Meanwhile, the mismatch between the modelled and observed velocity distributions in (a) at low velocities indicates that the Fid prescription overpredicts the number of low-velocity pulsars. The preferred NW prescription (b) reduces the disagreement in the low-velocity regime by removing NSs from ECSNe in non-interacting binaries.

We devise a statistic in order to quantify the goodness of fit between model and data. For a specified cutoff velocity $v_{\textrm{cut}}$, let $f_{\textrm{low}}$ represent the fraction of low-speed pulsars with transverse velocity $\leq$ $v_{\textrm{cut}}$, as predicted by a given model variation. We choose a value of $v_{\textrm{cut}}=50\,\mathrm{\rm{km}\,\rm{s}^{-1}}$. While this value is somewhat ad hoc, it serves to differentiate between high kicks, which disrupt binaries and eject pulsars from globular clusters, and low kicks, which do not. Let $N_{\textrm{slow}}$ be the number of observed pulsars with velocity $\leq$ $v_{\textrm{cut}}$ out of $N_{P}$. The probability of observing $N_{\textrm{slow}}$ low-velocity pulsars out of $N_{P}$ is described by the binomial distribution. The probability of observing $N_{\textrm{slow}}$ pulsars or fewer is thus given by

In practice, $N_{\textrm{slow}}$ is not known precisely due to the measurement uncertainty in observed pulsar transverse velocities, so it is marginalized out. Under the null hypothesis that the data are drawn from the model variation, ${P({\leq}N_{\textrm{slow}}|f_{\textrm{low}})}$ is uniformly distributed, allowing us to rule out any variations which yield very small or very large values for ${P({\leq}N_{\textrm{slow}}|f_{\textrm{low}})}$.

${P({\leq}N_{\textrm{slow}}|f_{\textrm{low}})}$ is plotted in Fig. 3 in black. The simulated $f_{\textrm{low}}$ values are shown as colored vertical lines; the uncertainty due to the finite simulated population is within the width of the lines. The null hypothesis is in tension with the data when a model, specified by $f_{\textrm{low}}$, intersects the black curve outside the indicated $\pm 2\sigma$ interval (see Table 4 for exact values of ${P({\leq}N_{\textrm{slow}}|f_{\textrm{low}})}$).

Fig. 3(a) shows $f_{\textrm{low}}$ values calculated after including low-velocity pulsars from both CCSNe and ECSNe. Model variations with $f_{\textrm{low}}$ values which cross the ${P({\leq}N_{\textrm{slow}}|f_{\textrm{low}})}$ curve below the $-2\sigma$ threshold over-produce low-velocity pulsars from a combination of ECSNe and CCSNe. This is the case for the Fid and K1 variations. Meanwhile, both variations which mask out ECSNe in non-interacting binaries cross well above the threshold and cannot be ruled out.

In order to distinguish the relative importance of ECSNe kicks and low-velocity CCSNe kicks, in Fig. 3(b), we consider an alternative population in which we include only EC-pulsars as low-velocity pulsars. Explicitly, we assume that all CC-pulsars are given very large kicks. Now, only the K1 model still crosses below the $-2\sigma$ cutoff, so it is the only model that can be confidently ruled out on the basis that it over-produces low-velocity pulsars from ECSNe alone. The low-velocity EC-pulsars in this model predominately come from binaries disrupted by Blaauw kicks, though ${\sim}30\%$ are in intact, very wide binaries.

The Fid model now falls just within the $2\sigma$ confidence interval and cannot be confidently ruled out. However, the high production of low-velocity EC-pulsars alone indicates a likely tension with observations. The two models which mask non-interacting ECSNe, NW and NW-K1, now nearly cross above the $+2\sigma$ level, suggesting that they significantly under-produce low-velocity pulsars. Since we have intentionally removed CC-pulsars here, this merely suggests that some fraction of low-velocity pulsars may be contributed by CC-pulsars.