Pulsar Kick: Status and Perspective

1. Introduction

Astrophysical compact objects like pulsars serve as remarkable cosmic markers, deepening our understanding of cosmic evolution and enabling tests of physics beyond the Standard Model (SM). Pulsars are rapidly rotating, magnetized Neutron Stars (NSs), first discovered by J. Bell-Burnell through observations of a radio source [1,2]. Pulsars provide unique advantages over laboratory-based experiments, functioning as natural, highly sensitive detectors in space. Their extreme environments, expansive space as a detection medium, prolonged exposure times, and their complementary role to terrestrial experiments make pulsars particularly powerful tools for studying the fundamental physics at large scales. These characteristics collectively enhance interaction probabilities and sensitivity by several orders of magnitude compared to traditional laboratory setups.

Binary pulsars, in particular, represent exceptional laboratories for testing Einstein’s General Relativity (GR) under intense gravitational fields. Observations of orbital period decay in binary compact star systems provided the first indirect evidence of Gravitational Waves (GWs), and their precise timing measurements allow for predictions with high accuracy, consistently aligned with GR [3,4].

NSs exhibit extreme densities and pressures at their cores, reaching around ${10}^{18}\mathrm{kg}/{\mathrm{m}}^3$ [5] and approximately ${10}^{34}\mathrm{Pa}$ [6], respectively. Under such conditions, quarks and gluons may exist in a deconfined state called quark-gluon plasma, a phase thought to have existed in the early universe just microseconds after the Big Bang [7]. This state, likely unique to NSs, enables pulsars to probe conditions otherwise unattainable in the observable universe.

A typical pulsar has a mass of approximately $1.4M_{\odot }$, a radius near $10\mathrm{km}$, and a magnetic field strength between ${10}^{14}-{10}^{16}\mathrm{G}$ [8]. Newly formed pulsars begin with a temperature around $\sim {10}^{12}\mathrm{K}$, which rapidly decreases to about ${10}^6\mathrm{K}$ over a few years due to substantial energy loss through neutrino emission [5]. Pulsars emit electromagnetic radiation at highly regular intervals, functioning as precise cosmic clocks.

Pulsars also play a crucial role in constraining the Equation of State (EoS) of nuclear matter. Measurements of pulsar mass, radius, and tidal deformability provide essential insights into the properties of ultra-dense matter. Pulsar Timing Arrays (PTAs), utilizing the extreme periodicity of pulsars, are used to detect the Stochastic Gravitational Wave Background (SGWB), potentially arising from SuperMassive Black Holes (SMBHs) [9] or exotic physics from the early universe [10].

Pulsars’ strong magnetic fields create a unique platform for studying magnetohydrodynamics and plasma physics under extreme conditions. The gravitational fields of pulsars may capture Dark Matter (DM), altering their properties, which could affect spin-down rates and spectral lines through DM-nucleon interactions [11,12,13,14,15,16,17]. Precise measurements thus serve as valuable tests for new physics scenarios, offering unique insights into the fundamental structure of the universe.

When a supergiant star reaches the end of its life, it undergoes a supernova explosion, collapsing inward and expelling its outer layers to form a gas and dust cloud known as a SuperNova Remnant (SNR). The core collapses into a NS or pulsar, which is why most NSs and pulsars are observed near the centers of SNRs. In cases where the supernova collapse is asymmetric, the resulting pulsar can receive a "kick" in a specific direction. In [18], a statistical analysis of the proper motion of 233 pulsars reveals that the highest $15\%$ of the velocity distribution reaches speeds exceeding 1000 km/s—an order of magnitude greater than that of their progenitors [19,20]. This phenomenon, known as the pulsar kick, provides a unique window into the underlying physics of supernova explosions [21], neutrino interactions [22], dense matter properties [23,24], and even the dynamics of galaxy evolution [25,26]. Pulsar velocities can be measured using timing residuals, Doppler shifts, and parallax methods. Additionally, the expansion of associated SNRs can provide valuable insights into pulsar velocity.

Monte Carlo simulations estimate the average birth velocity of pulsars to be around $\sim (250-300)\mathrm{km}/\mathrm{s}$, following a Maxwellian distribution with a dispersion of approximately $\sigma \sim 190\mathrm{km}/\mathrm{s}$ [27]. Pulsars with longer characteristic ages exhibit an asymmetric drift, which suggests they are dynamically older. This population may emerge from the lower-velocity segment of younger pulsars, though influenced by their binary origins and the effects of evolving within the Galactic potential. The findings align well with the characteristics of binaries containing NSs and closely match the velocities predicted by numerical supernova simulations. The observed velocity distribution of these pulsars is shaped by both their origins in binary systems and the influence of their later movement through the Galaxy.

In the study by [28], the average transverse velocity of thirteen MilliSecond Pulsars (MSPs) was determined to be around $85\pm 13\mathrm{km}/\mathrm{s}$, which is about four times slower than that of typical pulsars. Unlike younger ordinary pulsars, which tend to move away from the Galactic plane, half of these MSPs are moving toward it. Notably, half of them have ages comparable to or exceeding the age of the Galactic disk, indicating they are much older.

Johnston et al [29] investigated the alignment between pulsar spin and velocity vectors by studying 25 pulsars at 1369 MHz, with 21 also observed at 3100 MHz. The study used precise polarization measurements, accurate Rotation Measures (RMs), and refined velocity data to assess the relationship between the pulsar’s rotation and velocity directions. Although this alignment is challenging to confirm due to the presence of orthogonal polarization modes, additional insights from optical and X-ray observations for specific pulsars—such as PSR B0656+14, the Crab, and Vela—help indicate that the velocity vector aligns with the rotation axis. Furthermore, many pulsars show linear polarization with the position angle predominantly perpendicular to the magnetic field lines. However, if the supernova kicks impart both linear and angular momentum to the newborn NS, it may obscure any pre-existing stellar alignment before the explosion. Thus, observations of double NS systems cannot conclusively rule out aligned kicks.

Figure 1 shows an image of the pulsar kick, with false colors representing brightness temperatures ranging from $5.5\mathrm{K}$ to a maximum of $8.9\mathrm{K}$. The image has an angular resolution of approximately $1^{\prime }$, with a field of view of about $1^{\circ }.9\times 1^{\circ }.1$. A green cross marks the geometric center of the SNR [30], while circles indicate the location of PSR J0002+6216 [31]. A faint emission tail extends from the pulsar back toward the SNR’s geometric center.

Several hypotheses have been proposed to explain these high birth velocities of pulsars. The most widely accepted explanation for pulsar kicks is anisotropy in the supernova explosion at the time of their birth [33,34,35,36,37,38,39]. However, 3D hydrodynamic simulations indicate that this mechanism struggles to produce kicks exceeding 1000 km/s [23,33,35,36,40]. Additionally, these models tend to predict a strong correlation between mass and velocity that is not supported by observations [33,35]. This also makes it challenging to explain low-mass NSs with nearly symmetric ejecta experiencing significant kicks [41]. Furthermore, the models have difficulty reproducing the observed bimodal velocity distribution [20,23].

Asymmetrical ejection of matter or, potentially, neutrinos during Core-Collapse SuperNovae (CCSNe) not only imparts a natal kick to newborn pulsars but also contributes to their initial spin at birth, as explored in [42]. This study examines how these asymmetric forces create both linear and rotational motion in pulsars, shaping their velocity and spin characteristics from the moment they form.

Anisotropic neutrino emission or asymmetries in the collapsing star’s magnetic field could produce such high speeds. Hydrodynamic instabilities within the supernova may impart net momentum to the NS, resulting in a kick. For progenitor stars in binary systems, the supernova explosion can disrupt the system, transferring momentum to the pulsar. Additionally, certain beyond-SM physics scenarios may contribute to the natal kick of pulsars. In the following, we provide an overview of some leading hypotheses to explain the phenomenon of pulsar kicks.

The review article is organized as follows: Section 2 explores various mechanisms proposed to explain the phenomenon of pulsar kicks. Finally, Section 3 provides concluding remarks and discusses future directions for utilizing pulsar kicks as a tool to investigate fundamental physics.

We use natural units $c=1$ (speed of light in vacuum) and $\hslash =1$ (reduced Planck constant) in the paper unless stated otherwise.

2. Mechanisms of Explaining Pulsar Kick

Below, we outline the leading theories proposed to explain pulsar kicks, along with their challenges and limitations. It is also possible that pulsar kicks result from a combination of all these contributing factors.

2.1. Asymmetric Neutrino Emission

During core collapse, a vast number of neutrinos are emitted. If this emission is slightly anisotropic, it creates a momentum imbalance that can impart a significant kick to the NS. This mechanism is particularly intriguing since even a minor asymmetry in neutrino emission can produce high velocities due to the enormous energy neutrinos carry. As a proto-NS cools, neutrinos carry away about $99\%$ of its gravitational binding energy, approximately ${10}^{53}\mathrm{erg}$, with a momentum of around ${10}^{43}\mathrm{g}·\mathrm{cm}/\mathrm{s}$. Therefore, only a $1\%$ anisotropy in the neutrino emission distribution is sufficient to explain the pulsar kick [43,44].

In the dense environment of a NS, electron neutrinos (${\nu }_e$) have a shorter mean free path than other neutrino flavors. If resonant oscillations occur, causing muon and tau neutrinos (${\nu }_{\mu ,\tau }$) to oscillate into ${\nu }_e$ above the ${\nu }_{\mu ,\tau }$ neutrinosphere but below the ${\nu }_e$ neutrinosphere, then ${\nu }_e$ will be absorbed by the medium due to charge and neutral current interactions. The resonance point effectively determines the ${\nu }_{\mu ,\tau }$ neutrinosphere position. When neutrinos propagate in the presence of a magnetic field, their self-energy acquires a term proportional to $\mathbf{B}·\mathbf{k}$, where $\mathbf{B}$ is the magnetic field and $\mathbf{k}$ the neutrino momentum vector. This interaction distorts the effective ${\nu }_{\mu ,\tau }$ neutrinosphere, depending on the orientation of $\mathbf{B}$ and $\mathbf{k}$. Consequently, the neutrinosphere deviates from a spherical shape, and neutrinos escaping along the magnetic field direction have different temperatures than those emitted in the opposite direction. This directional asymmetry in neutrino momentum can create enough imbalance to account for the observed pulsar kick [44].

The condition for neutrino oscillation resonance in the presence of a magnetic field is given as [44]

$${\displaystyle \frac{\Delta m^2}{2k}}cos2\theta =\sqrt{2}{G}_F{N}_e+{\displaystyle \frac{e{G}_F}{\sqrt{2}}}{\left({\displaystyle \frac{3N_e}{{\pi }^4}}\right)}^{{\displaystyle \frac{1}{3}}}\widehat{\mathbf{k}}·\mathbf{B},$$

(1)

where, $\Delta m^2$ represents the squared mass difference between neutrinos, $\theta$ is the neutrino mixing angle in a vacuum, $N_e$ is the electron number density, $G_F$ denotes the Fermi constant, e is the electron charge, and $\widehat{\mathbf{k}}$ is the unit vector in the direction of the neutrino momentum. In this context, the neutrino momentum asymmetry is obtained as [44]

$${\displaystyle \frac{\delta k}{k}}=1\%{\left({\displaystyle \frac{3\mathrm{MeV}}{T}}\right)}^2\left({\displaystyle \frac{B}{3\times {10}^{14}\mathrm{G}}}\right),$$

(2)

where, T represents the background temperature. To achieve a $1\%$ anisotropy in the neutrino momentum, the neutrino mass must be approximately $100\mathrm{eV}$ with a small mixing angle which should decay in the cosmological time scale not to overclose the universe. This analysis assumes that the resonance point is unaffected by the temperature distribution.

The main limitation of this mechanism is that it requires a heavy neutrino mass, which conflicts with results from neutrino oscillation experiments and cosmological constraints [45,46]. Corrections to the neutrino transport and hydrostatic equations in a proto-NS indicate that a magnetic field of at least $({10}^{15}-{10}^{16})\mathrm{G}$ [47] is needed to produce a $1\%$ asymmetry in the neutrino momentum. Further refinement of the neutrino oscillation mechanism for generating pulsar kicks has considered variations in neutrino opacities and the neutrino absorption process, ${\nu }_e+n\to e^{-}+p^{+}$ [48]. An asymmetry in this neutrino absorption on opposite sides of the effective neutrinosphere leads to a pulsar kick, differing from the process described above, where the kick results from the distortion of the effective neutrinosphere by the magnetic field where neutrinos traveling along the magnetic field have different momenta compared to those moving in the opposite direction relative to the field. A self-consistent model of neutrino transport in the stellar atmosphere suggests that an even stronger magnetic field, approximately ${10}^{17}\mathrm{G}$ near the star’s surface, is necessary to account for pulsar kicks through the oscillation mechanism [49]. Later studies suggest that resonant neutrino conversion may account for pulsar kicks, with momentum asymmetry tied to the logarithmic derivative of the energy flux. The kick arises from the radial flux variation and deformation of the resonance surface, primarily driven by the geometric distortion of the neutrinosphere rather than any change in a conserved plane-parallel flux by the star’s overall structure [50]. The fractional anisotropy of momenta in this scenario is given as [50]

$${\displaystyle \frac{\delta k}{k}}=-{\displaystyle \frac{1}{18}}{\displaystyle \frac{\delta }{R_r}},$$

(3)

where $R_r$ represents the radius of the resonance surface, while the value of $\delta$ depends on the proto-NS model, such as the polytropic [51,52,53,54] or spherical Eddington models [55]. Although, in this scenario, a magnetic field on the order of ${10}^{16}\mathrm{G}$ and a neutrino mass around $100\mathrm{eV}$ are required.

2.2. Sterile Neutrino, Dark Matter and Other New Physics Scenarios

2.2.1. Sterile Neutrino and Dark Matter

The need for a heavy active neutrino mass to explain the pulsar kick, which conflicts with neutrino oscillation and cosmological predictions, can be mitigated by considering active-sterile neutrino oscillations. For active-active oscillations (${\nu }_e\leftrightarrow {\nu }_{\mu ,\tau }$), thermalization occurs when incoming ${\nu }_{\mu ,\tau }$ either enter their own neutrinosphere within the proto-NS or reconvert to ${\nu }_e$ after crossing the resonance region, placing them within their neutrinosphere. In contrast, for sterile neutrinos, only the latter process applies. Outgoing neutrinos alone contribute to the net momentum flux, with asymmetry arising from the differential areas of hemispherical emission surfaces for neutrinos moving in opposite directions with equal momentum, generating a non-zero fractional momentum asymmetry. Explaining pulsar kicks through active-sterile neutrino oscillations requires a sterile neutrino mass around the $\mathrm{keV}$ scale, making it a potential candidate for warm DM, along with a very small active-sterile mixing angle[56,57]. Additionally, a magnetic field of approximately ${10}^{17}\mathrm{G}$ is needed, which could plausibly exist in the interior of a proto-NS.

Introducing ultralight DM in the proto-star environment can reduce the need for extreme magnetic fields [58]. The alignment of the DM velocity relative to neutrino momentum in different directions creates an asymmetry that could account for the pulsar kick. The ultralight scalar DM couples with neutrinos through derivative coupling as

$$-\mathcal{L}\supset g_{\alpha \beta }^{\prime }{\partial }_{\mu }\varphi {\overline{\nu }}_{\alpha }{\gamma }^{\mu }{\nu }_{\beta },$$

(4)

where $g_{\alpha \beta }^{\prime }$ denotes the derivative coupling between DM and active neutrinos, and $\varphi$ denotes the oscillating DM scalar field given as $\varphi =(\sqrt{2{\rho }_{\mathrm{DM}}}/m_{\varphi })cos(m_{\varphi }t+\theta )$. Here, $m_{\varphi }$ denotes the mass of the scalar field, $\theta$ denotes a random phase and ${\rho }_{\mathrm{DM}}$ denotes the DM density. Therefore, the dispersion relation for neutrino becomes

$$E_{\nu }=\left|\mathbf{k}\right|-g_{\alpha \beta }^{\prime }\dot{\varphi }+g_{\alpha \beta }^{\prime }\widehat{\mathbf{k}}·\nabla \varphi +{\displaystyle \frac{m^2}{2\left|\mathbf{k}\right|}}+\mathcal{O}\left({g_{\alpha \beta }^{\prime }}^2{(\nabla \varphi )}^2\right).$$

(5)

The third term in the right hand side of Equation (5) is responsible for the pulsar kick. Note, the scalar coupling of $\varphi$ with neutrino does not change the neutrino momentum and the corresponding coupling will not give any effect on the pulsar kick.

The resonance condition for active-sterile neutrino oscillations within an ultralight scalar DM background is expressed as [58]

$${\displaystyle \frac{\Delta m^2}{2k}}cos2\theta =V_{\nu }+g_{\alpha \beta }^{\prime }\sqrt{2{\rho }_{\mathrm{DM}}}\widehat{\mathbf{k}}·{\mathbf{v}}_{\varphi },$$

(6)

where $V_{\nu }$ represents the standard matter potential, and ${\mathbf{v}}_{\varphi }$ is the DM velocity vector. The coupling of scalar with sterile neutrinos is assumed to be small. The resulting neutrino momentum asymmetry under a DM background becomes [58]

$${\displaystyle \frac{\delta k}{k}}={\displaystyle \frac{g_{\alpha \beta }^{\prime }{v}_{\varphi }\sqrt{2{\rho }_{\mathrm{DM}}}}{18V_{\nu }{R}_r(h_p^{-1}+h_{V_{\nu }}^{-1})}},$$

(7)

where $h_p^{-1}$ and $h_{V_{\nu }}^{-1}$ are quantities that depend on the proto-star model. Below, we discuss one of the proto-star models, specifically the polytropic EoS.

• The Polytropic model: The isotropic neutrinosphere in the polytropic model is described by the following equations [51,52,53,54]

  $$\begin{array}{ccc}\hfill {\displaystyle \frac{dP\left(r\right)}{dr}}& =& -{\displaystyle \frac{GM\left(r\right)\rho \left(r\right)}{r^2}},\left(\mathrm{hydrodynamical}\mathrm{equilibrium}\right)\hfill \end{array}$$

  (8)

  $$\begin{array}{ccc}\hfill F\left(r\right)& =& -{\displaystyle \frac{1}{36}}{\displaystyle \frac{1}{\kappa \rho \left(r\right)}}{\displaystyle \frac{d{T}^2}{dr}},\left(\mathrm{energy}\mathrm{transport}\right)\hfill \end{array}$$

  (9)

  $$\begin{array}{ccc}\hfill F\left(r\right)& =& {\displaystyle \frac{L_c}{4\pi r^2}},\left(\mathrm{flux}\mathrm{conservation}\right)\hfill \end{array}$$

  (10)

  where P denotes the pressure, F denotes the flux, G denotes Newton’s gravitational constant, $M\left(r\right)=4\pi {\int }_0^r{r^{\prime }}^{2}d{r}^{\prime }{\rho }_T\left(r^{\prime }\right)$, ${\rho }_T$ is the total density of matter, and $L_c$ denotes the luminosity of the proto-star. Assuming the proto-star is filled with relativistic nucleons of polytropic gas, the EoS can be written as

  $$P\left(r\right)=K{\rho }^{\Gamma },K={\displaystyle \frac{T_c}{m_n{\rho }_c^{1/3}}}=5.6\times {10}^{-5}{\mathrm{MeV}}^{-4/3},$$

  (11)

  where $\Gamma$, $m_n$, $T_c$, and ${\rho }_c$ denote the adiabatic index, nucleon mass, temperature and density of the core respectively and $T_c=40\mathrm{MeV}$ and ${\rho }_c={10}^{14}\mathrm{gm}/{\mathrm{cm}}^3$. Using Equations (8) and (11), we can write

  $${\displaystyle \frac{d{\rho }^{\Gamma -1}}{dr}}=-{\displaystyle \frac{{\lambda }_{\Gamma }{r}_c{\rho }_c^{\Gamma -1}}{M_c}}{\displaystyle \frac{M\left(r\right)}{r^2}},$$

  (12)

  where $r_c$, ${\rho }_c$ and $M_c$ denote the radius, density and mass of the core of the proto-star. The quantity ${\lambda }_{\Gamma }$ is given as

  $${\lambda }_{\Gamma }=(G{M}_c/r_c{\rho }_c^{\Gamma -1})(\Gamma -1)/K\Gamma \simeq 0.29{\displaystyle \frac{2G{M}_c}{\mathrm{km}}}{\displaystyle \frac{10\mathrm{km}}{r_c}}{\displaystyle \frac{40\mathrm{MeV}}{T_c}},$$

  (13)

  where $\Gamma =4/3$ and $M_c=1M_{\odot }$. The approximate solution of Equation (12) is obtained as

  $${\rho }^{\Gamma -1}\left(r\right)={\rho }_c^{\Gamma -1}\left[{\lambda }_{\Gamma }\left({\displaystyle \frac{r_c}{r}}-1\right)m\left(r\right)+1\right],$$

  (14)

  where $m\left(r\right)=\mu +(1-\mu )r_c/r$ with $m\left(r_c\right)=1$. Also, by defining $x=r_c/r$, and $\alpha =(1-\mu ){\lambda }_{\Gamma }$, $\beta =(2\mu -1){\lambda }_{\Gamma }$, $\gamma =1-\mu {\lambda }_{\Gamma }$, we write Equation (14) as

  $${\rho }^{\Gamma -1}\left(x\right)={\rho }_c^{\Gamma -1}(\alpha x^2+\beta x+\gamma ).$$

  (15)
  The parameter $\mu$ is obtained by setting the condition $\rho \left(R_s\right)=0$, where $R_s$ is the radius of the star and we obtain

  $$\mu =\left[{\displaystyle \frac{R_s}{{\lambda }_{\Gamma }(R_s-r_c)}}-{\displaystyle \frac{r_c}{R_s}}\right]{\displaystyle \frac{R_s}{R_s-r_c}}.$$

  (16)
  Also, combining Equations (9) and (10) we can write the expression of the temperature profile in terms of density distribution as

  $${\displaystyle \frac{d{T}^2}{dr}}=-{\displaystyle \frac{9\kappa L_c}{\pi r^2}}\rho \left(r\right),$$

  (17)

  where the core luminosity is $L_c\sim 9.5\times {10}^{51}\mathrm{erg}/\mathrm{s}$ and $\kappa \sim 5.6\times {10}^{-9}{\mathrm{cm}}^4/{\mathrm{erg}}^3{\mathrm{s}}^2$. The solution of Equation (17) is obtained as

  $$T\left(r\right)=T_c\sqrt{2{\lambda }_c\left[\chi (r_c/r)-\chi \left(1\right)+1\right]},$$

  (18)

  where

  $${\lambda }_c={\displaystyle \frac{9}{2\pi }}{\displaystyle \frac{\kappa L_c{\rho }_c}{T_c^2{r}_c}}\sim 1.95{\displaystyle \frac{{\rho }_c}{{10}^{14}\mathrm{g}/{\mathrm{cm}}^3}}{\displaystyle \frac{10\mathrm{km}}{r_c}}{\left({\displaystyle \frac{40\mathrm{MeV}}{T_c}}\right)}^2,$$

  (19)

  and $\chi$ is a polynomial function x given as

$$\begin{array}{c}\hfill \chi \left(x\right)={\gamma }^{3}x+{\displaystyle \frac{3}{2}}\beta {\gamma }^2{x}^2+\gamma (\alpha \gamma +{\beta }^2)x^3+{\displaystyle \frac{\beta }{4}}(6\alpha \gamma +{\beta }^2)x^4+\\ \hfill +{\displaystyle \frac{3\alpha }{5}}(\alpha \gamma +{\beta }^2)x^5+{\displaystyle \frac{\beta {\alpha }^2}{2}}{x}^6+{\displaystyle \frac{{\alpha }^3}{7}}{x}^7.\end{array}$$

(20)

To achieve a 1% momentum asymmetry required for the pulsar kick, the coupling and mass parameters are estimated as $g^{\prime }\sim 1.23\times {10}^{-11}{\mathrm{eV}}^{-1}$ and $m_{\varphi }\sim 1.02\times {10}^{-14}\mathrm{eV}$. The asymmetric neutrino emission during core collapse induces a lasting strain in the surrounding spacetime, detectable by second- and third-generation GW detectors [59]. The maximum characteristic GW amplitude from this neutrino emission is given by [58,59]

$$\left|h\left(t\right)\right|\simeq {\displaystyle \frac{2G}{r}}\alpha L_{\nu }\tau \lesssim 1.06\times {10}^{-19}\left({\displaystyle \frac{\alpha }{0.01}}\right)\left({\displaystyle \frac{E_{\mathrm{tot}}}{2\times {10}^{53}\mathrm{erg}}}\right)\left({\displaystyle \frac{1\mathrm{kpc}}{r}}\right),$$

(21)

where $\alpha$ represents the anisotropy parameter, $L_{\nu }$ the neutrino luminosity, and r the distance from the source. The characteristic neutrino emission timescale $\tau$ is defined as $E_{\mathrm{tot}}=L_{\nu }\tau$.

The anisotropy parameter depends on the DM mass and coupling, and alternative DM models, such as vector ultralight DM, can also be constrained by pulsar kick observations. Consequently, GW detectors offer a unique avenue for probing DM interactions.

Figure 2 illustrates the variation of GW intensity with frequency, providing insights into pulsar kicks driven by asymmetric neutrino emission induced by ultralight DM. For, ultralight vector DM, $\alpha =1.8\times {10}^{-4}{\sigma }_1$, where ${\sigma }_1=g_{\alpha \beta }^{\prime }/m_A^{\prime }{x}_{\mathrm{res}}\eta$ and for ultralight scalar DM, $\alpha =1.8\times {10}^{-7}{\sigma }_2$ where ${\sigma }_2=g_{\alpha \beta }^{\prime }/x_{\mathrm{res}}\eta$. Here, $m_{A^{\prime }}$ denotes the mass of ultralight vector DM, $x_{\mathrm{res}}=r_c/r_{\mathrm{res}}$ and $\eta <1$ is a parameter depends on the proto-star model. The GW intensity decreases with increasing radiation frequency, while higher ultralight DM-neutrino coupling enhances the wave intensity. Signals from nearby supernovae with appropriate rotation frequencies are expected to be detectable by future GW detectors. The neutrino-induced gravitational memory effect in the presence of ultralight DM, resulting from a supernova, lies below the sensitivity of current detectors like LIGO and Virgo but can be observed by second- and third-generation detectors. Rapidly rotating NSs with larger anisotropy parameters may be detected by advanced detectors such as aLIGO and aVirgo, whereas slowly rotating NSs are within the sensitivity range of LISA. GW detectors operating in the deci-Hz range [60,61,62,63], such as DECIGO (DECi-hertz Interferometer Gravitational Wave Observatory) [64] and BBO (Big Bang Observer) [64], which achieve a sensitivity of GW strain around ${10}^{-24}$ could detect memory signals from pulsar kicks. Furthermore, the Einstein Telescope (ET) [65], with its unprecedented sensitivity (${10}^{-24}-{10}^{-25}$), will explore GW signals from cosmological distances, enabling studies of supernova neutrinos and the associated gravitational memory effects. These advanced detectors will also be capable of probing smaller anisotropy parameters for asymmetric neutrino emissions from sources as distant as a few megaparsecs. Detecting such gravitational memory signals will offer valuable constraints on sterile neutrino masses, ultralight DM interactions, and other beyond SM scenarios.

2.2.2. Majoron Emission

If, rather than ultralight DM or sterile neutrinos, active neutrinos couple to a massless pseudoscalar majoron, an asymmetric emission of majorons could account for the pulsar kick, provided that majorons carry away at least $10\%$ of the star’s total binding energy [66]. A magnetic field around ${10}^{16}\mathrm{G}$ is necessary for this explanation. In such a strong magnetic field, electrons in the medium become polarized, introducing an additional term to the neutrino potential given by [66]

$$\delta V=-\sqrt{2}{G}_F{n}_B\langle {\lambda }_e\rangle cos{\theta }_{k}diag\left({\displaystyle \frac{3}{2}},{\displaystyle \frac{3}{2}},{\displaystyle \frac{1}{2}}\right),$$

(22)

where $n_B$ represents the baryon density, $Y_e=(n_e-{\overline{n}}_e)/n_B$ is the electron fraction, ${\theta }_k$ is the angle between the neutrino momentum and polarization direction, and $\langle {\lambda }_e\rangle$ denotes the fraction of polarized electrons, defined as [66]

$$\langle {\overrightarrow{\lambda }}_e\rangle =-{\displaystyle \frac{e\overrightarrow{B}}{2}}{\left({\displaystyle \frac{3}{{\pi }^4}}\right)}^{1/3}{n}_e^{-2/3}.$$

(23)

Since $\delta V$ depends on the direction of the neutrino momentum, the majoron emission processes, $\nu \nu \to \varphi$ and $\overline{\nu }\to \nu \varphi$, will vary based on the initial neutrino directions, leading to anisotropic emission. The magnetic field strength necessary to account for the pulsar kick via the process ${\nu }_e{\nu }_e\to \varphi$ is given by [66]

$$|\overrightarrow{B}|=3\times {10}^{16}\mathrm{G}{\displaystyle \frac{M}{1.4M_{\odot }}}{\displaystyle \frac{v}{500\mathrm{km}/\mathrm{s}}}{\displaystyle \frac{3\times {10}^{53}\mathrm{erg}}{E_{\mathrm{tot}}}}{\displaystyle \frac{V_{\nu }}{0.5\mathrm{eV}}}{\left({\displaystyle \frac{0.05{\mathrm{fm}}^{-3}}{n_e}}\right)}^{1/3}{\displaystyle \frac{0.5}{x}}.$$

(24)

2.2.3. Neutrino Spin-Flavour Oscillation

The pulsar kick mechanism can also be attributed to the spin polarization of matter in a supernova’s magnetic field, incorporating neutrino spin-flavor precession. This magnetic field distorts the resonance surface for neutrino spin-flavor transitions, leading to anisotropic neutrino emissions. In cases where resonant spin-flavor precession occurs into sterile neutrino states, the magnetic field strength required to produce the pulsar kick is approximately $B\gtrsim 4\times {10}^{15}\mathrm{G}$, assuming a Dirac transition magnetic moment of ${\mu }_{\nu }\gtrsim {10}^{-14}{\mu }_B$ [67]. Likewise, for resonant spin-flavor precession between active neutrinos and anti-neutrinos, a magnetic field of $B\gtrsim 2\times {10}^{16}\mathrm{G}$ is needed, with a Majorana transition magnetic moment of ${\mu }_{\nu }\gtrsim {10}^{-15}{\mu }_B$.

The chirality flip of neutrinos, driven by the presence of a neutrino magnetic moment in strange quark matter NSs, offers an explanation for the pulsar kick. For this mechanism to account for the average pulsar velocity, the required neutrino magnetic moment is approximately ${\mu }_{\nu }\sim 3.6\times {10}^{-18}{\mu }_B$ [68].

The spin-gravity coupling of neutrinos in the gravitational field of a rotating NS can lead to a pulsar kick. Taking into account the magnetic and gravitational interactions, the spin flavour conversion for neutrinos is governed by the evolution equation [69]

$$\iota {\displaystyle \frac{d}{d\lambda }}\left(\begin{array}{c}{\nu }_{fL}\\ {\nu }_{f^{\prime }L}\\ {\nu }_{fR}\\ {\nu }_{f^{\prime }R}\end{array}\right)=H\left(\begin{array}{c}{\nu }_{fL}\\ {\nu }_{f^{\prime }L}\\ {\nu }_{fR}\\ {\nu }_{f^{\prime }R}\end{array}\right),$$

(25)

where the effective Hamiltonian in the chiral basis can be written as [69]

$$H=\left(\begin{array}{cc}{H}_L& H_{f{f}^{\prime }}^{*}\\ H_{f{f}^{\prime }}& H_R\end{array}\right),$$

(26)

and

$$H_L=\left(\begin{array}{cc}{V}_{{\nu }_f}+{\Omega }_G-\delta c_2& \delta s_2\\ \delta s_2& V_{{\nu }_{f^{\prime }}}+{\Omega }_G+\delta c_2\end{array}\right),$$

(27)

$$H_R=\left(\begin{array}{cc}-\delta c_2& \delta s_2\\ \delta s_2& \delta c_2\end{array}\right),H_{f{f}^{\prime }}=B_{\perp }\left(\begin{array}{cc}{\mu }_{ff}& {\mu }_{f{f}^{\prime }}\\ {\mu }_{f{f}^{\prime }}& {\mu }_{f^{\prime }{f}^{\prime }}\end{array}\right),$$

(28)

where $c_2=cos2\theta$, $s_2=sin2\theta$, $\delta =(m_2^2-m_1^2)/4E$, $m_{1,2}$ are the two neutrino mass eigenstates, $B_{\perp }=Bsin\alpha$ is the component of the magnetic field orthogonal to neutrino momentum, and [69]

$${\Omega }_G\left(r\right)=8\times {10}^{-13}\left({\displaystyle \frac{M}{M_{\odot }}}\right){\left({\displaystyle \frac{R}{10\mathrm{km}}}\right)}^2{\left({\displaystyle \frac{10\mathrm{km}}{r}}\right)}^3{\displaystyle \frac{{\omega }^{\prime }cos\beta }{{10}^4\mathrm{Hz}}}\mathrm{eV},$$

(29)

where ${\omega }^{\prime }$ is related to the angular velocity of NS and $\beta$ is the angle between the neutrino momentum and the angular velocity of NS. The anisotropy in neutrino emission arises from the alignment of neutrino momenta relative to the star’s angular velocity. The fraction of pulsar momentum asymmetry in this scenario is given as [69]

$${\displaystyle \frac{\delta k}{k}}={\displaystyle \frac{4GM{R}^2{\omega }^{\prime }}{5r^3}}{\displaystyle \frac{1}{9V_{{\nu }_f}{R}_r(h_p^{-1}+h_{V_{{\nu }_f}}^{-1})}}.$$

(30)

In this scenario, spin-flavor oscillation is considered for calculating the momentum asymmetry of the pulsar, with a required neutrino magnetic moment of approximately ${\mu }_{\nu }\lesssim {10}^{-11}{\mu }_B$ to account for the kick.

2.2.4. Lorentz and CPT Violation

Pulsar kick can also be explained from the active-sterile neutrino oscillation in the Standard Model Extension (SME) framework, where the Lorentz and/or CPT invariance can be violated. The neutrino dispersion relation in such framework can be written as [70]

$$E\simeq k+{\displaystyle \frac{m^2}{2k}}+\Omega ,$$

(31)

where

$$\Omega =-{\displaystyle \frac{c_L^{\mu \nu }{k}_{\mu }{k}_{\nu }}{k}}+{\displaystyle \frac{a_{L\mu }{k}^{\mu }}{k}},$$

(32)

where c violates Lorentz invariance and a violates both CPT and Lorentz invariance. The resonance condition for the active-sterile neutrino oscillation in this SME framework can be written as [70]

$$2\delta c_2=V_{{\nu }_f}+\delta \Pi ·\widehat{\mathbf{k}},$$

(33)

where $\delta \Pi ·\widehat{\mathbf{k}}=(\delta \mathbf{a}-2k\delta \mathbf{c})·\widehat{\mathbf{k}}$ and $\delta \mathbf{c}={\mathbf{c}}^{\left({\nu }_f\right)}-{\mathbf{c}}^{\left({\nu }_s\right)}$, $\delta \mathbf{a}={\mathbf{a}}^{\left({\nu }_f\right)}-{\mathbf{a}}^{\left({\nu }_s\right)}$. Therefore, the fractional asymmetry in the neutrino momenta can be written as [70]

$${\displaystyle \frac{\delta k}{k}}=\left(2k\right|\delta {\mathbf{c}}_{\mathbf{L}}|+|\delta {\mathbf{a}}_{\mathbf{L}}\left|\right){\displaystyle \frac{1}{9R_r{V}_{{\nu }_f}(h_p^{-1}+h_{V_{{\nu }_f}}^{-1})}}.$$

(34)

For $1\%$ of momentum asymmetry of pulsar, the SME parameters can be constrained as $\delta {\mathbf{c}}_L\lesssim 4\times {10}^{-17}$ and $\delta {\mathbf{a}}_L\lesssim 1.6\times {10}^{-9}\mathrm{eV}$. The SME parameters and the ultralight scalar DM parameters are related in the pulsar kick scenario as $a_L=g^{\prime }|\nabla \varphi |$ and $c_L=g^{\prime }|\nabla \varphi |/2k$.

2.2.5. Massless Neutrino Framework

A pulsar kick can be explained even with massless neutrinos if a strong magnetic field is present, assuming a slight violation of weak universality. In this scenario, one does not need any neutrino magnetic moment.

The observed pulsar velocity can also be attributed to tensor neutrino-gravity coupling, which causes anisotropic neutrino emission under strong magnetic fields if neutrinos violate the equivalence principle at a level of ${10}^{-9}-{10}^{-10}$ [71]. In scenarios where neutrinos interact with the gradient of the gravitational potential (as can occur in string theory), the dimensional parameter characterizing this Violation of Equivalence Principle (VEP) needed to account for the pulsar kick is around $({10}^{-2}-{10}^{-3})\mathrm{cm}$.

Neutrino oscillations induced by VEP, based on the Parametrized Post-Newtonian (PPN) framework, can account for the pulsar kick [72]. In this scenario, VEP effects generate anisotropies in both the linear and angular momenta of emitted neutrinos, potentially explaining both the translational and rotational motions observed in pulsars. The massless neutrino dispersion relation under these conditions is given by [73]

$$E=k\left[1+h_{oi}{\widehat{\mathbf{k}}}_i-{\displaystyle \frac{1}{2}}{h}_{ij}{\widehat{\mathbf{k}}}_i{\widehat{\mathbf{k}}}_j-{\displaystyle \frac{1}{2}}{h}_{00}U\right],$$

(35)

where

$$h_{00}=2{\gamma }^{\prime }U+\mathcal{O}\left({\omega }^4\right),$$

(36)

$$h_{0i}=-{\displaystyle \frac{7}{2}}{\Delta }_1{V}_i-{\displaystyle \frac{1}{2}}{\Delta }_2{W}_i+\left({\alpha }_2-{\displaystyle \frac{1}{2}}{\alpha }_1\right)v_{i}U-{\alpha }_2{v}_j{U}_{ji}+\mathcal{O}\left({\omega }^4\right),$$

(37)

$$h_{ij}=2\gamma U{\delta }_{ij}+\Gamma U_{ij}+\mathcal{O}\left({\omega }^4\right).$$

(38)

Here, $\gamma$, ${\gamma }^{\prime }$, $\Gamma$, ${\Delta }_{1,2}$, ${\alpha }_{1,2}$ and $\mathbf{v}$ denote the PPN expansion parameters. Also, $\omega$ denotes the velocity of the source. The parameters ${\alpha }_{1,2}$ vanish in the Lorentz covariant theories however, if there exists a preferred frame characterize by velocity $\mathbf{v}$, then ${\alpha }_{1,2}$ are non-zero. Also, U, $V_i$, $W_i$, $U_{ji}$ are the potential terms. In this scenario, the neutrino evolution equation can be written as [72]

$$\iota {\displaystyle \frac{d}{dr}}\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\end{array}\right)={\displaystyle \frac{{\Delta }_0}{2}}\left(\begin{array}{cc}-cos2{\theta }_g& sin2{\theta }_g\\ sin2{\theta }_g& cos2{\theta }_g\end{array}\right)\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\end{array}\right),$$

(39)

where ${\theta }_g$ denotes the mixing angle between the flavour and gravitational eigenstate and ${\Delta }_0$ for a rotating proto-NS can be written as [72]

$$\begin{array}{c}\hfill {\Delta }_0=[-(\delta {\gamma }^{\prime }+\delta \gamma )U-\delta \Gamma J-\delta \Gamma I{(\widehat{\mathbf{r}}·\widehat{\mathbf{k}})}^2+\left\{(\delta {\alpha }_2-{\displaystyle \frac{1}{2}}\delta {\alpha }_1)U-\delta {\alpha }_{2}J\right\}\mathbf{v}·\widehat{\mathbf{k}}-\\ \hfill \delta {\alpha }_{2}I(\widehat{\mathbf{r}}·\mathbf{v})\widehat{\mathbf{r}}·\widehat{\mathbf{k}}-{\displaystyle \frac{1}{2}}(7\delta {\Delta }_1+\delta {\Delta }_2)J\Omega \times \mathbf{r}·\widehat{\mathbf{k}}]E,\end{array}$$

(40)

where $\delta \gamma ={\gamma }^2-{\gamma }^1$, I and J are related to the potentials and ${\Delta }_0$ plays the same role as $(m_2^2-m_1^2)/4E$.

2.3. Hydrodynamic Instabilities

The globally anisotropic ejection of mass and subsequent NS acceleration in a supernova core can result from convective overturn and low-order oscillations in the neutrino-heated layer, which is called the Standing Accretion Shock instability (SASI). Specifically, low-mode $(l=1,2)$ convection arises from initial random perturbations, particularly when the explosion onset is gradual, allowing convective structures to merge [74]. This kick mechanism does not strictly align the pulsar’s spin with its spatial velocity, although the rotation axis serves as a natural orientation, potentially limiting extreme misalignment. Both gravitational and hydrodynamic forces from the asymmetric supernova ejecta propel the NS over a timescale exceeding one second, achieving speeds above $500\mathrm{km}/\mathrm{s}$.

Figure 3 reveals that NSs in the high-velocity, high-acceleration group (darker shading) continue to accelerate to much higher velocities beyond one second, whereas only a few low-velocity stars exceed 200 km/s [74]. After one second of post-bounce evolution, simulated for most models, NS velocities reach up to 800 km/s. The models split into two distinct populations: one with velocities below 200 km/s and low acceleration, and another with velocities above 200 km/s, exhibiting higher average acceleration.

Hydrodynamic recoil resulting from neutrino-driven CCSNe provides a straightforward mechanism for accelerating NSs and pulsars, eliminating the need for assumptions of anisotropic neutrino emission or other exotic models. In a two-dimensional radiation-hydrodynamic simulation, the collapse leads to substantial acceleration of a newly formed NS as a result of a neutrino-driven explosion. During this explosion, a $10\%$ anisotropy in the ejected low-mass, high-velocity material causes a recoil in the high-mass NS, allowing the pulsar to reach recoil speeds of $150\mathrm{km}/\mathrm{s}$. This mechanism offers a consistent explanation for pulsar velocities based on the dynamics of supernova ejecta [39].

The asymmetry of ejecta in Supernova Remnants (SNRs) provides a possible explanation for the observed pulsar kicks [75]. Researchers have studied the connection between NS kick velocities and the X-ray morphologies of 18 SNRs using data from the Chandra X-ray Observatory and ROSAT. They analyzed SNR asymmetries through a power-ratio method, which is a multipole expansion technique that focuses on dipole, quadrupole, and octupole power ratios to quantify asymmetries. The results showed no significant correlation between the power ratio magnitudes and NS kick velocities overall. However, for cases like Cas A and G292.0+1.8, where the X-ray emissions primarily represent ejecta distribution, the associated NSs tended to move in the opposite direction to the bulk of X-ray emissions. Similarly, PKS 1209–51, CTB 109, and Puppis A also exhibited this trend, but since their emissions are dominated by circumstellar or interstellar material, these asymmetries may not directly reflect the distribution of the original supernova ejecta.

Figure 4 presents the quadrupole power ratio ($P_2/P_0$, left) and octupole power ratio ($P_3/P_0$, right) plotted against NS transverse velocities, using the center of emission as the origin for the multipole expansion [75]. No significant correlation is observed between $P_2/P_0$ or $P_3/P_0$ and NS velocities. This finding remains consistent for both proper motion estimates and geometric offset velocity estimates, even when accounting for the considerable uncertainties in the geometric measurements. Additionally, neither non-interacting nor interacting SNRs, represented as circles and squares, respectively, in Figure 4, exhibit a clear trend in power ratios relative to NS velocity. However, it is noteworthy that most SNRs with high $P_2/P_0$ and $P_3/P_0$ values are associated with interacting remnants.

The findings support the hypothesis that NS kicks are more likely due to asymmetries in ejecta rather than anisotropic neutrino emissions. In summary, no relationship was found between NS velocities and power ratio magnitudes within the 0.5–2.1 keV X-ray band, regardless of whether the center of emission or explosion site was chosen as the reference point for the multipole expansions. This indicates that the small-scale X-ray asymmetries and the NS kicks may result from distinct processes or be influenced by different factors.

The simulation has also explored the possibility of explaining high recoil velocities through a neutrino-driven supernova explosion mechanism, where non-radial hydrodynamic instabilities break symmetry and induce strong kicks. However, these studies are based on two-dimensional simulations that assume axial symmetry, with the polar axis treated as a coordinate singularity, restricting fluid flow across it.

Hydrodynamical simulations offer valuable insights into pulsar kicks, yet have limitations. High-resolution, three-dimensional simulations of hydrodynamics and neutrino transport are computationally intensive, which restricts accuracy. While simulations can account for pulsar speeds of a few hundred km/s, achieving velocities above $1000\mathrm{km}/\mathrm{s}$ remains challenging and may require magnetic or neutrino-driven anisotropies. Additionally, the results are highly sensitive to initial conditions, such as the progenitor’s rotation and magnetic field configuration, which influence large-scale asymmetries.

2.4. Other Mechanisms

• Chiral anisotropy conversion: Pulsar kicks can be accounted for by the anisotropic emission of neutrinos, which arises from their scattering with the background axial electron current—a result of the chiral separation effect [76]. Achieving a pulsar recoil requires anisotropy in either the magnetic field or density in momentum space. In this framework, a magnetic field strength of approximately ${10}^{16}\mathrm{G}$ can drive pulsar velocities exceeding $1000\mathrm{km}/\mathrm{s}$.
• Evanescent proto-neutron star: If a CCSN results in a rapidly rotating proto-NS, it may subsequently cool and undergo fragmentation, forming a binary proto-NS system in a very close orbit [77]. In this scenario, the lighter companion could eventually be tidally disrupted, imparting a significant kick to the remaining proto-NS. This mechanism has the potential to generate kick velocities exceeding $1000\mathrm{km}/\mathrm{s}$.
• Rocket effect: An asymmetry in the magnetic field configuration of a pulsar’s strong magnetic field could generate a small, continuous electromagnetic force. If the magnetic moment is misaligned with the pulsar’s rotation axis, this could lead to a "rocket effect," exerting a gradual push on the pulsar [78]. However, a significant drawback of this mechanism is that it cannot achieve kick velocities as high as those produced by neutrino-driven or hydrodynamic processes.

3. Conclusions and Discussions

The velocity distribution of pulsars is strongly influenced by the internal structure of the proto-NS. Various EoS have been proposed to model the proto-star’s evolution, including the polytropic model, the spherical Eddington model, and the Murnaghan EoS [79]. These equations describe the relationship between pressure, density, and temperature within the star. However, each of these models has its own inherent limitations and assumptions.

The theory that pulsar kicks arise from asymmetric neutrino emission does not predict a direct correlation between velocity and magnetic field ($B-V$ correlation) because the magnetic field relevant for the kick is located inside the hot NS during the initial seconds after the supernova collapse. Astronomical observations, however, can only infer the surface magnetic field of NSs millions of years after this event. The magnetic fields at early and late stages of NS evolution are not straightforwardly connected due to the complex changes that occur during the star’s cooling process. Future X-ray observatories like XRISM will have the capability to detect weak emission lines from trace elements, potentially providing new insights into linking pulsar kick theories with observational evidence. For example, the magnetic field at the central region has been considered as ${10}^{18}\mathrm{G}$. The magnetic field of the proto-NS depends on the baryon density and its parametrized as [80,81,82]

$$B_{\mathrm{PNS}}=B_s+B_c[1-e^{-\beta {(\rho /{\rho }_s)}^{\gamma }}],$$

(41)

where $B_s={10}^{14}\mathrm{G}$, $B_c={10}^{18}\mathrm{G}$ are the magnetic fields at the surface and the core respectively, $\beta ={10}^{-5}$, $\gamma =3$, and ${\rho }_s={10}^{11}\mathrm{g}/{\mathrm{cm}}^3$. Thus, the requirement of a strong magnetic field to drive neutrino-induced kicks is not completely ruled out.

Since no strong correlation exists between pulsar velocity and other pulsar properties, statistical analyses of pulsar populations neither confirm nor disprove any existing models in the literature. Most of the proposed mechanism for pulsar kicks does not depend on a correlation between pulsar velocity and magnetic field. These result aligns with recent studies suggesting that a B-V correlation is unnecessary to explain pulsar kicks [27,28,83].

One benefit of explaining pulsar kicks through the hypothesis of massless neutrinos and equivalence principle violation is that it avoids conflicts with cosmological neutrino constraints, unlike the model involving active neutrino oscillations proposed in [44].

While establishing a direct relationship between the magnetic field and pulsar velocity is challenging, the duration of neutrino or majoron emission significantly exceeds the initial rotation period of a newborn pulsar. This extended emission period suggests that there may be an alignment between the rotation axis and the direction of the pulsar kick [84,85,86]. Although only a few pulsars have known rotational axes that appear to align with their kick directions, this alignment trend hints at a potential correlation. However, further data is essential to reach a definitive conclusion on this connection.

It remains uncertain to what extent pronounced $l=1$ modes in the ejecta distribution and sustained downflows of matter onto the NS can develop in a fully 3D environment, or how frequently these features might occur, though initial 3D simulations using the current setup and input physics show promise. Additionally, the expected distribution of NS recoil velocities from 3D models is still unclear. The hydrodynamic kick mechanism is inherently stochastic and requires extended simulation durations, but the vast computational resources needed make such large-scale, long-term simulations currently infeasible.

Studies [77,83,84,87,88,89] suggest that low-velocity NSs tend to have space velocities aligned with their spin axes, while high-velocity NSs move primarily perpendicular to their spins, resulting in two distinct (low- and high-velocity) distributions.

This review addresses the longstanding challenge of explaining pulsar kicks in astrophysics, along with potential solutions to this phenomenon. Pulsars receive a high-velocity kick shortly after their formation in supernova explosions, with birth speeds ranging from several hundred to thousands of kilometers per second, often placing them far from their origin. The three leading explanations for pulsar kicks are asymmetric supernova explosions, asymmetric neutrino emissions, and hydrodynamic instabilities. However, no single mechanism has fully accounted for the observed pulsar velocity distribution. It is likely that a combination of these mechanisms contributes in varying degrees, depending on the specific conditions of each supernova. Advancements in supernova observations, more refined simulations, and deeper insights into neutrino physics in supernova environment may improve our understanding of these powerful kicks. Additionally, studying pulsar kicks provides a valuable opportunity to constrain potential new physics scenarios.