EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 687 (July 2024) A&A, 687 (2024) A272 Full HTML
Open Access
Issue
A&A
Volume 687, July 2024
Article Number A272
Number of page(s) 15
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202449996
Published online 25 July 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

Many astrophysical objects experience a velocity kick during their evolution. The most notable case is that of neutron stars, and in particular pulsars, which have space velocities significantly higher than their progenitors and thus must have received an additional velocity – a natal kick – at birth (Gunn & Ostriker 1970; Lyne & Lorimer 1994; Cordes & Chernoff 1998). Other likely recipients of kicks include stellar-mass black holes (Jonker & Nelemans 2004; Fragos et al. 2009; Repetto et al. 2012; Mandel 2016; Andrews & Kalogera 2022) and binary systems, since binaries receive a systemic kick due to mass loss after one of the two components goes supernova (Blaauw 1961; Hills 1983; Brandt & Podsiadlowski 1995; Van den Heuvel et al. 2000; Willems et al. 2004; Zhang et al. 2013; Atri et al. 2019; Gandhi et al. 2019; Fortin et al. 2022; Zhao et al. 2023).

There remain several open questions regarding these kicks, the most pressing of which for most systems is their magnitudes. For example, some studies found that pulsars were best described with a single-peaked distribution with relatively high velocities (e.g. Hobbs et al. 2005; Faucher-Giguère & Kaspi 2006), while more recent studies found a double-peaked distribution with a second significant peak at lower values (Arzoumanian et al. 2002; Verbunt et al. 2017; Igoshev 2020). The distribution of kick magnitudes, in turn, has implications not only for our understanding of the nature of the kicks, for example whether they are due to binary disruption or are intrinsic to an asymmetric supernova (Dewey & Cordes 1987; Van Paradijs & White 1995; Iben & Tutukov 1996; Fryer & Kalogera 1997; Van den Heuvel & Van Paradijs 1997), but also for the progenitors’ evolution (e.g. Kalogera 1996; Podsiadlowski et al. 2004; Van den Heuvel 2007; Beniamini & Piran 2016), the supernova mechanism (e.g. Janka 2012, 2013; Wongwathanarat et al. 2013; Coleman & Burrows 2022; Burrows et al. 2024), the merger rates of compact object binaries (e.g. Dominik et al. 2012; Giacobbo & Mapelli 2018; Kruckow et al. 2018; Vigna-Gómez et al. 2018; Iorio et al. 2023), and the locations of the transients they produce (e.g. Fryer et al. 1999; Perna & Belczynski 2002; Voss & Tauris 2003; Belczynski et al. 2006; Church et al. 2011; Jiang et al. 2020; Zevin et al. 2020; Mandhai et al. 2022).

Kicks kinematically decouple a population of objects from their progenitors, thus determining their present-day locations and velocities. In the context of Galactic pulsars, several authors have modelled their Galactic trajectories in order to constrain the phase-space distribution and compare it to observations. This method has been used to check whether pulsars could have been a viable gamma-ray burst progenitor (Hartmann et al. 1990; Paczynski 1990), to constrain the evolution of pulsars spins and magnetic fields (Bhattacharya et al. 1992; Hartman et al. 1997), and to study the population of isolated old neutron stars that might become observable due to accretion of the interstellar medium (Blaes & Rajagopal 1991; Blaes & Madau 1993; Zane et al. 1995; Popov et al. 2000; Treves et al. 2000). More recent works modelling the positions and velocities of pulsars across the Galaxy include those of Kiel & Hurley (2009), Ofek (2009), Sartore et al. (2010), and Sweeney et al. (2022).

The distribution of velocities in a population of kicked objects evolves with time, but more importantly, it also varies with the Galactic location where it is probed. That is to say, the velocity distribution we might be able to observe could be different from that of the total population. Lyne & Lorimer (1994), for instance, observed that young pulsars have higher velocities than older ones but do not believe this is caused by an actual dynamical evolution. Instead, they link it to the work of Cordes (1986), who argued that high-velocity pulsars rapidly move away from the solar neighbourhood, reducing their likelihood of detection and resulting in an observational bias towards lower velocities for older pulsars. Hansen & Phinney (1997) investigated this argument, using a Monte Carlo simulation of masses moving in a Galactic potential and analysing the galactocentric velocities of the masses that get close to the solar neighbourhood. They find that: (1) at first, these masses have a high velocity, (2) after approximately 10 Myr, the radial velocity of the set increases due to masses born close to the Galactic centre arriving at the Solar System, and (3) after approximately 100 Myr, only the low-velocity tail of the distribution is left.

In this paper we are interested in tracking the apparent velocity shifts first observed in pulsars by Lyne & Lorimer (1994), and understanding their origin. We aim to investigate how kicks affect the velocity of astrophysical objects, and how these effects can be detected in the solar neighbourhood at times that may be long after the kick was imparted to the object. To examine the velocity evolution of kicked objects, we used a method similar to that used by Hansen & Phinney (1997): a Monte Carlo simulation of point masses moving through the Galactic potential (also, cf. Sartore et al. 2010, and references therein).

This paper is structured as follows. In Sect. 2 we describe our model and simulation. Then, in Sect. 3, we give our main results, which show that the deceleration of kicked objects due to the Galactic potential is dominant in shaping the final velocity distribution. In Sect. 4 we explain these results by using grids of simulations that exemplify the dynamics. Moreover, in Sect. 5 we compare our results to 2D and (estimated) 3D velocities of pulsars. Finally, in Sect. 6 we summarise our findings and their implications.

2. Simulation

First, we describe our model and simulation. In particular, we describe our Milky Way (MW) model, which is used to seed the point masses in our Monte Carlo simulation (Sect. 2.1). Then, we elaborate on how we initialised and evaluated the trajectories of the objects (Sect. 2.2).

2.1. Milky Way model

To describe the initial positions of the point masses, we used the MW model of Gaspari et al. (2024), which is based on the model of Chrimes et al. (2021) and entails a synthetic MW image that we used to sample the initial positions of the simulated objects. The MW model consists of three parts: an exponential disc, a bar-like bulge, and spiral arms. Here, we give a brief overview of these components, while more detailed descriptions are given by Gaspari et al. (2024) and Chrimes et al. (2021).

The bulge is modelled with the triaxial boxy Gaussian distribution of Grady et al. (2020), where the bar makes a 27° angle with the line connecting the Solar System to the Galactic centre, in the xy plane (Wegg & Gerhard 2013). For the spiral arms, we used the data from Chrimes et al. (2021), who used the Reid et al. (2016, 2019) method to map the spiral arms, using water and methanol masers. In our 3D MW model, we set the 2D spiral arms of Chrimes et al. (2021) at z = 0. Also, since their arm-profile only covers the lower half of the MW, we rotated them for the upper half. For the disc, we used an exponential density profile:

ρ disc ( x , y , z ) ∝ exp ( − x 2 + y 2 R d ) exp ( − | z | z d ) , $$ \begin{aligned} \rho _{\mathrm{disc}}\left(x,{ y},z\right)\propto \exp \left(-\dfrac{\sqrt{x^2+{ y}^2}}{R_d}\right)\exp \left(-\dfrac{|z|}{z_d}\right), \end{aligned} $$(1)

where we chose Rd  =  2.6 kpc and zd  =  0.3 kpc (Bland-Hawthorn & Gerhard 2016). The density profiles of the arms, bulge, and disc were used to create a synthetic image of the MW, with a scale of 250 pc per pixel, and create the synthetic MW image by showing the densities on a logarithmic scale (effectively converting luminosity to magnitude).

The three components in this image are weighted based on their luminosity. We used the I-band luminosities from Flynn et al. (2006), who estimated a luminosity of 1 × 1010 L⊙ for the bulge and a luminosity of 3 × 1010 L⊙ for the combination of spiral arms and exponential disc. Also, we chose an arm strength of 0.15 (Gaspari et al. 2024), meaning the arms get a luminosity of 0.45 × 1010 L⊙ and the disc gets a luminosity of 2.55 × 1010 L⊙. The density profiles are normalised and multiplied by these luminosity. Moreover, we added a Gaussian blur with a full width at half maximum of 2 pixels (similar to Chrimes et al. 2021) in order to slightly broaden the arms and create a more realistic model.

To seed the initial positions of the point masses in our simulation, we used this MW model. We seeded our point masses in the thin disc since kicks often occur in the context of massive – and therefore young – stars, and this is where the population of young, massive stars is situated (e.g. Feltzing & Bensby 2008). We therefore set the initial positions of our point masses at z = 0 because the scale height of the thin disc is approximately the size of one pixel. Furthermore, we sampled the initial positions in the xy plane from a combination of the spiral arms and the exponential disc (Eq. (1), where z = 0). Figure 1 shows the synthetic MW image and the sampled initial positions. Our simulation consists of 104 point masses, of which 500 are shown in the figure. The image contains only the lower half of the xy plane, but the points are sampled in the upper half as well (which uses a rotated version of the lower half).

thumbnail Fig. 1.

Initial positions of point masses in our simulation (blue dots), plotted over the synthetic MW image, based on the works of Gaspari et al. (2024) and Chrimes et al. (2021), as described in the main text. The figure contains 500 of the 104 simulated point masses. The Sun symbol represents the Solar System, at R⊙ = 8.122 kpc and z⊙ = 20.8 pc.

2.2. Trajectories

After the initial positions have been sampled, the initial velocities of the points are determined. To model the kick, we used the kick velocity distribution from Hobbs et al. (2005), who estimated that the kick velocities of pulsars are best described by a Maxwellian with σ = 265 km s−1. We do not argue that their kick distribution accurately describes pulsar kicks, but used their work to show how such a distribution evolves over time. Indeed, as we show in Appendix B, the final velocity distribution is not strongly dependent on the shape of the initially assumed kick distribution. For each point, we sampled a kick velocity (vkick) from the Hobbs et al. (2005) distribution. Assuming the kick direction is isotropic, we oriented the kick vector by uniformly sampling an angle α within the xy plane, and sampling an angle β with a probability proportional to cos β, which is the declination relative to the xy plane. This way, the kick vectors are oriented isotropically.

We also assumed the points are in a circular orbit around the Galactic centre before they experience the kick, meaning they get an additional circular velocity (vcirc), which has a magnitude determined by the gravitational potential. We used the MW potential from McMillan (2017) to determine these circular velocities.

After combining the kick and the circular velocity, the initial velocity vector (v0), in the cylindrical coordinates of azimuth (ϕ), galactocentric radius (R), and vertical height (z), becomes

v 0 = ( v ϕ v R v z ) = ( v kick sin α cos β + v circ v kick cos α cos β v kick sin β ) . $$ \begin{aligned} \boldsymbol{v}_0=\begin{pmatrix}v_{\phi }\\ v_{R}\\ v_{z}\end{pmatrix}=\left(\begin{array}{@l@}v_{\mathrm{kick}}\sin \alpha \cos \beta +v_{\mathrm{circ}}\\ v_{\mathrm{kick}}\cos \alpha \cos \beta \\ v_{\mathrm{kick}}\sin \beta \end{array}\right). \end{aligned} $$(2)

Having determined the initial positions and velocities of the point masses, we evolved their trajectory in the McMillan (2017) potential using the GALPY1 v.1.9.0 package (Bovy 2015). In initialising the potential, we assumed that the circular velocity at the position of the Sun equals the azimuthal velocity of the Sun.

We evolved the trajectories of the point masses for 200 Myr, while evaluating the position and velocity of each point at several times during their trajectory. When we evaluated the trajectories at a certain time, we considered three populations in our simulation:

  • The solar neighbourhood, which we define as all points in the simulation that come within a distance of 2 kpc to a solar orbit. We determined this by initialising 12 orbits with the initial conditions of the Solar System, each one azimuthally rotated by 2π/12. Effectively, this creates 12 non-overlapping spheres with a radius of 2 kpc, moving in a Sun-like orbit (cf. Hansen & Phinney 1997). Azimuthal rotation does not affect the results significantly, due to the MW model being almost cylindrically symmetric. We used the solar velocity,

    v ⊙ = ( v ϕ v R v z ) = ( 245.6 km s − 1 ( GRAVITY Collab . et al . 2018 ) − 12.9 km s − 1 ( Drimmel & Poggio 2018 ) 7.78 km s − 1 ( Reid & Brunthaler 2004 ) ) , $$ \begin{aligned} \boldsymbol{v}_{\odot }=\begin{pmatrix}v_{\phi }\\ v_{R}\\ v_{z}\end{pmatrix}=\left(\begin{array}{@l@}245.6\,\mathrm{km\,s^{-1}}\,(\mathrm{GRAVITY\ Collab.\ et\ al.}\ 2018)\\ {-}12.9\,\mathrm{km\,s^{-1}}\,(\mathrm{Drimmel\ \& \ Poggio}\ 2018)\\ 7.78\,\mathrm{km\,s^{-1}}\,(\mathrm{Reid\ \& \ Brunthaler}\ 2004)\end{array}\right), \end{aligned} $$(3)

    and the solar position (rotated for each nth sphere, where 0 ≤ n < 12),

    r ⊙ = ( ϕ R z ) = ( n · 2 π / 12 8.122 kpc ( GRAVITY Collab . et al . 2018 ) 20.8 pc ( Bennett & Bovy 2019 ) ) , $$ \begin{aligned} \boldsymbol{r}_{\odot }=\begin{pmatrix}\phi \\ R\\ z\end{pmatrix}=\left(\begin{array}{@l@}n\cdot 2\pi /12\\ 8.122\,\mathrm{kpc}\,(\mathrm{GRAVITY\ Collab.\ et\ al.}\ 2018)\\ 20.8\,\mathrm{pc}\,(\mathrm{Bennett\ \& \ Bovy}\ 2019)\end{array}\right), \end{aligned} $$(4)

    as initial conditions for these solar spheres. These values correspond to the ones implemented in GALPY v.1.9.0. We co-evolved these solar orbits and at each point in time evaluated the galactocentric velocities of the points that are within one of these spheres.

  • The thick disc. As a rough approximation of the thick disc, based on Fig. 1, we defined this region as R ≤ 15 kpc and |z|≤2 kpc.

  • The total population, including all point masses, even the ones that get a high kick and escape the Galaxy.

Considering these three populations allows us to compare the velocities of the points that get close to a solar orbit – and therefore have a higher probability of being detected – to the velocities of the other points. This way, we can estimate the role of observational bias in context of, for example, older pulsars appearing to have lower velocities.

3. Deceleration

Analysing the results from the simulation described in Sect. 2, we first investigated how the velocity distribution changes over time (Sect. 3.1), and then examined how the components of the velocity vector (i.e. vϕ, vR, and vz) evolve as the points move through the Galaxy (Sect. 3.2).

3.1. Velocity magnitude

First, we note that at t = 0 Myr, the galactocentric velocity distribution of the point masses does not precisely match the Hobbs et al. (2005) kick distribution (in all three populations). This is caused by the circular velocity that is added to the kick velocity (i.e. the vcirc term in Eq. (2)). In other words, the velocity distribution is shifted to higher velocities at first, compared to the kick distribution, and then decreases to lower velocities. We find that if we consider all ages up to 20 Myr, the total distribution averages out to approximately the kick distribution. For this reason, we compared the velocities below 20 Myr to the ones above 20 Myr.

This comparison considers the distributions of the galactocentric velocities, meaning they are, for example, not corrected for the motion of the Solar System. In general, we find that the velocity distributions show a deceleration in all three populations. It is important to note here that with “deceleration”, we are referring to the fact that after about 20 Myr, the objects tend to have velocities lower than their initial velocity, this does not mean that the objects at these points in time are necessarily decelerating (i.e. have dv/dt < 0). Figure 2 displays the three distributions, for ages below and above 20 Myr, and indeed shows that for ages below 20 Myr, the velocities approximate the kick distribution, while above 20 Myr the velocities are significantly lower. In Appendix A, we zoom in on the decrease in the velocity distribution, and show how it remains stable beyond 200 Myr up to at least 1 Gyr.

thumbnail Fig. 2.

Distributions of galactocentric velocities (v = |v|, top row) and peculiar velocities in the local standard of rest of the point mass (vLSR, bottom row) for the simulated point masses. For the first 20 Myr, we evaluated the velocities every 1 Myr and show the total, normalised distribution in histograms with a bin width of 20 km s−1 (light blue curves). After 20 Myr, we evaluated the velocities every 10 Myr and also show the total, normalised distribution in similar histograms (dark blue curves). The dashed black curves show the kick distribution (i.e. a Maxwellian with σ = 265 km s−1, as used by Hobbs et al. 2005). The three panels show the evaluated velocities for the three populations: the solar neighbourhood (left), the thick disc (centre), and the total population (right). These populations have sizes of 24 817, 144 496, and 210 000 for ages below 20 Myr, and 7079, 54 320, and 180 000 for ages above 20 Myr, respectively.

Figure 2 shows that the total distributions below 20 Myr approximate the kick distribution. The distributions, however, start to shift to lower velocities, becoming relatively stable after 20–30 Myr. Interestingly, the decelerated velocity distributions look similar in all three populations: they show a peak around 200 km s−1. The main difference between the three populations can be found in the high-velocity tail, which is significantly smaller in the solar neighbourhood. This is due to the objects that get a large kick, escape the Galaxy, and are therefore not found in the solar neighbourhood but still contribute to the total population. The thick disc also shows more high-velocity objects than the solar neighbourhood, which is probably caused by objects that get a significant kick but cannot escape the Galaxy and fall back to the Galactic centre.

The bottom row of Fig. 2 shows the distributions of vLSR, which equals the velocity in the local standard of rest (LSR) frame. That is, for an object this velocity is calculated by subtracting the local vcirc from its velocity vector. The distributions of vLSR look similar to the galactocentric velocities, and show a similar deceleration. The deceleration of the kicked objects is therefore not limited to the galactocentric frame.

We argue that because the distributions of the point masses above 20 Myr decelerate and look similar in all three populations, it is difficult to explain the deceleration with an observational bias (as proposed by Cordes 1986; Lyne & Lorimer 1994). After all, the difference between the three populations is mainly the high-velocity tail, and if the deceleration was merely apparent, and the result of an observational bias, it should only be present in the solar neighbourhood and absent from the total population.

In the left panel of Fig. 3, we show how the velocities after 20 Myr relate to the initial velocities. The main deceleration, as can be seen in the figure, occurs for velocities below approximately 500 km s−1. Here, a higher v0 results in a lower velocity after 20 Myr. This can be explained by the fact that higher kicks disturb the initial circular orbits more significantly and cause the objects to follow new, eccentric, orbits with an apocentre at higher galactocentric radii. In these orbits, the objects spend most of their time near the apocentre, where they have velocities lower than their initial circular velocity, due to the gravitational potential. This is why objects in the high-velocity tail (i.e. v > 500 km s−1) retain velocities closer to their initial velocity: these objects get a sufficiently large kick to escape the Galaxy and are therefore less affected by the Galactic potential if they get a higher kick. The right panel in Fig. 3 shows the galactocentric radius after 20 Myr, compared to the initial galactocentric radius (R0), and it indeed shows that 83% of the objects are, after 20 Myr, situated at radii larger than their initial R0.

thumbnail Fig. 3.

Initial conditions versus decelerated conditions, for the total population of our simulation. The left panel shows the initial velocities of the objects (v0) versus the decelerated velocities (vdec), which we define as any velocity for t > 20 Myr. We evaluated the decreased velocities every 10 Myr and show the resulting distribution in a 2D histogram with bins of 20 km s−1. The right panel shows the initial galactocentric radial position of the objects (R0) versus the galactocentric radial position of the objects for t > 20 Myr (Rdec), evaluated every 10 Myr and shown in a 2D histogram with bins of 0.5 kpc. Here, 83% of the objects are located at Rdec > R0.

We therefore argue that interaction with the gravitational potential of the Galaxy can provide a better explanation, because (1) the deceleration is prevalent in all three populations, (2) the kicks disturb the initial circular orbits of the objects, after which they obtain eccentric orbits in which they spend most of their time at higher galactocentric radii, where they have velocities lower than their initial velocity (due to the Galactic potential), and (3) this effect does not depend on the Hobbs et al. (2005) distribution being exactly right, since Maxwellians with a different σ show a similar deceleration (in Appendix B, we show how varying the σ of these Maxwellians gives indistinguishable decelerated distributions in the solar neighbourhood).

The hypothesis that the MW potential determines this deceleration is also supported by the initial galactocentric radii of the objects that arrive in the solar neighbourhood. Figure 4 shows the R0 distributions of the objects in the solar neighbourhood, both below and above 20 Myr. At t = 0 Myr, all objects in the solar neighbourhood are born at an initial galactocentric radius where |R0 − R⊙|≤2, and the distribution for age ≤20 Myr still shows a vast majority of objects born in this region. For age > 20 Myr, however, enough time has passed for the kicked objects to have travelled relatively large distances, because of which the initial bias towards locally born objects disappears and the distribution starts to resemble the initial radii of the total population. After 20 Myr, therefore, the solar neighbourhood is a representative sample of the total population, both in velocity magnitude and birth locations. In other words: after a certain amount of time, the objects in the solar neighbourhood are mostly born significantly closer to the Galactic centre, meaning they migrated to larger radii and lost kinetic energy due to the Galactic potential.

thumbnail Fig. 4.

Initial galactocentric radii (R0) of the objects that are in the solar neighbourhood below 20 Myr (light blue) and above 20 Myr (dark blue), shown in normalised histograms with bins of 1 kpc. The dotted curve shows the initial radii of the total population, in a similar histogram. The distribution of the total population peaks slightly above the half-light radius, Rd = 2.6 kpc (Eq. (1)), because the spiral arms shift the initial positions to higher radii compared to the disc.

3.2. Velocity direction

Besides the deceleration in velocity magnitude, we are also interested in how the cylindrical components of the velocity vectors (i.e. vϕ, vR, and vz) change over time. In Fig. 5, we show the evolution of the velocity magnitude and the separate components. The top row in the figure shows how the velocity magnitude distributions evolve, and how they become relatively stable after 20–30 Myr. Also, in the total population the high-velocity tail is visible, as we also found in Fig. 2.

thumbnail Fig. 5.

Evolution of the galactocentric velocity (v) and the cylindrical components of the velocity vector (vϕ, vR, and vz). The figure shows normalised distributions, in 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1, 20 km s−1, and 10 km s−1 for the solar neighbourhood, the thick disc, and the total population, respectively. For clarity, we set the maximum values of the distributions at 4 × 10−5 s km−1 Myr−1 ≈ 13 × 10−19 km−1, even though the peaks in the last row exceed this maximum to a certain degree. The black lines show the mean velocities at each point in time.

The second row shows the azimuthal velocity (vϕ), which at t = 0 Myr is not symmetrically distributed around vϕ = 0 km s−1: the added vcirc term in Eq. (2) shifts the distribution to positive vϕ. After this initial shift, however, the azimuthal velocities decrease. In the total population, these velocities are lower than in the solar neighbourhood or thick disc, because this population includes the point masses that escape the Galaxy and travel to large distances, meaning the centripetal force caused by the gravitational attraction of the Galaxy becomes relatively small.

The third row shows the radial velocity (vR). Interestingly, these distributions start symmetrically distributed around vR = 0 km s−1, but start to increase to positive velocity, before decreasing again after approximately 20 Myr. Here, the total population differs in the fact that the deceleration does not cause a significant fraction of the sample to get negative vR, as opposed to the solar neighbourhood and the thick disc. Again, this is due to the subset of point masses that get a large kick and escape (and therefore retain a positive vR).

The fourth, and final, row shows the vertical velocities (vz). Similar to vR, the vertical velocities start symmetrically distributed around vz = 0 km s−1. However, in all three populations, both positive and negative vz decrease significantly, even within the first 10 Myr. This creates a narrow peak around vz = 0 km s−1.

In general, we note the large similarities between the velocity distributions in all three populations, and argue that this implies the population of objects that get close to the Solar System is a relatively representative sample of the total population. Therefore, we argue that the observational bias, as proposed by Cordes (1986) and Lyne & Lorimer (1994), is insufficient to explain the fact that older pulsars are observed to have lower velocities. Instead, the apparent deceleration seems to be caused by the actual deceleration of objects as their orbit is changed by the kick, after which they spend most of their time at larger galactocentric radii, effectively exchanging kinetic energy for potential energy.

4. Dynamics

Our simulation shows several interesting aspects of the velocity evolution of kicked objects, such as the general deceleration and the initial increase in vR. To provide explanations for these aspects, we built a grid of simulations for different R0 and vkick so that the effects of these parameters can be examined (Sect. 4.1). Also, we simulated objects that get a total initial velocity completely in one of the cylindrical directions, to show how the direction of the initial velocity affects the trajectory and velocity evolution (Sect. 4.2).

4.1. Grid

We are interested in how the different initial galactocentric radii (R0) and kick velocities influence the velocity evolution. To do this, we used a grid of smaller simulations (103 objects instead of 104), for several values of R0 (i.e. 1 kpc, 5 kpc, and 9 kpc, instead of the initial positions shown in Fig. 1) and vkick (i.e. 100 km s−1, 400 km s−1, and 700 km s−1, instead of sampling from the Hobbs et al. 2005 distribution). Figure 6 shows the velocity evolution of the total population for each value of R0 and vkick. The initial distributions (i.e. at t = 0 Myr) range from |vcirc(R0)−vkick| to vcirc(R0)+vkick, because of the angle between the initial circular velocity at R0 and the kick. For t > 0 Myr, the objects can have velocities outside of this range.

thumbnail Fig. 6.

Grid of simulations of 103 kicked objects, with a fixed initial galactocentric radius, R0 (rows: 1 kpc, 5 kpc, and 9 kpc), and a single value for vkick (columns: 100 km s−1, 400 km s−1, and 700 km s−1). The figure shows the evolution of the galactocentric velocity (v) of the total population, in 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1. Even though these simulations use a single value for vkick, the velocity distributions do not start as delta functions, because of the vcirc term in Eq. (2).

The figure shows that low kicks (i.e. 100 km s−1) perturb the initial circular orbits of the objects only slightly, because of which the objects start to oscillate around their initial orbit. However, larger kicks (i.e. 400 km s−1 and 700 km s−1) perturb the initial orbits significantly, meaning the objects start to migrate through the Galaxy. That is, they will spend most of their time at galactocentric radii larger than the initial R0, causing the deceleration. The figure also shows that this deceleration is more rapid at lower R0. After all, the velocities of objects born near the Galactic centre decrease more rapidly as they move outwards, due to the shape of the MW potential.

Figure 7 shows the galactocentric velocities of the objects that are in the solar neighbourhood, for the grid of simulations. For vkick = 100 km s−1, objects can only enter the solar neighbourhood if their perturbed orbit intersects with the region. When vkick = 100 km s−1, this is not the case for initial radii, R0 = 1 kpc, but for R0 = 5 kpc the figure shows that part of the orbits do intersect with the solar neighbourhood. This part is even larger if R0 = 9 kpc, since these objects are born within the solar neighbourhood. When vkick = 400 km s−1, the objects that arrive in the solar neighbourhood will have more eccentric orbits, and the figure shows the deceleration (similar to the deceleration shown in Figs. 5 and 6). For vkick = 700 km s−1, most objects get a velocity greater than the escape velocity of the MW, which means they temporarily cross the solar neighbourhood, but most will not return to this region.

thumbnail Fig. 7.

Grid of simulations of 103 kicked objects, similar to Fig. 6, but showing only the objects that are in the solar neighbourhood. Again, these simulations have a fixed initial galactocentric radius, R0 (rows: 1 kpc, 5 kpc, and 9 kpc), and a single value for vkick (columns: 100 km s−1, 400 km s−1, and 700 km s−1). Similarly to Fig. 6, the evolution of the galactocentric velocity (v) is shown through 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1. However, since the sizes of the populations in the panels differ significantly, they are shown on different colour scales.

These objects, the ones that receive a high kick, influence the high-velocity tail, and cause the main difference between the solar neighbourhood population and the total population. Since the high-kicks only influence the high-velocity tail, and this tail does not contribute significantly to the population of objects observed in the solar neighbourhood (as shown in Figs. 2, 5, and 7), the initial kick distribution can be changed greatly without altering the distribution of objects older than 20 Myr, observed in the solar neighbourhood (as we show in Appendix B).

Moreover, Fig. 7 shows how different initial conditions can bring objects to the solar neighbourhood with identical galactocentric velocity. For instance, both when R0 = 1 kpc and vkick = 400 km s−1, and when R0 = 9 kpc and vkick = 100 km s−1, objects with a velocity of approximately 200 km s−1 can be found in the solar neighbourhood. This shows that an observed velocity in the solar neighbourhood of about 200 km s−1 could both be the result of a 100 km s−1 kick as well as a 400 km s−1 kick. This means that, besides obvious statements such as the fact that an observed object of any velocity is extremely unlikely to be the result of a 100 km s−1 kick at R0 = 1 kpc, we find that in order to determine the initial conditions of an individual observed object (of unknown age) based on its velocity magnitude, one needs to assume (1) the distribution of initial positions, and (2) the distribution of kick velocities. In other words, one needs to weigh the simulations in our grid according to assumed distributions of R0 and vkick, before being able to determine the probable origin of an observed object with a certain velocity. The speed of an object alone contains little information about the kick at birth.

4.2. Extrema

We are interested in explaining why the direction of the velocity vectors changes over time, as shown in Fig. 5. In particular, we are interested in explaining why initially vR increases, while vϕ and vz decrease immediately. To do this, we investigated the special cases where the initial velocity is pointed entirely in the ϕ or (negative) R directions. In other words, in these two dimensions we examined the extrema of the orientation of the initial velocity vector (i.e. the initial velocities, being a combination of vkick and vcirc as defined in Eq. (2), pointed entirely in one of the cylindrical directions).

We hypothesise that the reason vR initially increases while vϕ and vz decrease is the fact that the orientation of the velocity vector of a kicked object changes within the cylindrical coordinate system, because of their motion through the Galactic potential. The cylindrical coordinate system can, then, favour reorientation in the positive R direction, causing the acceleration. In Fig. 8, we sketch how initial velocity vectors oriented in the ϕ and negative R direction can be reoriented by their motion through the Galactic potential in such a way that they obtain a positive vR component.

thumbnail Fig. 8.

Schematic diagrams of hypothetical paths (light orange curves) of objects (orange dots) resulting from an initial velocity (orange arrows) completely in the azimuthal (left panel) or negative radial (central panel) direction. The shown initial velocity vectors and paths are merely qualitative estimates of possible paths that potentially contribute to the initial acceleration of vR (as found in Fig. 5). In Fig. 9, we show examples of orbits that are a result of initial velocities in these directions.

The figure shows the following two effects, which can reorient the velocity vector in the positive R direction (as sketched in the panels of Fig. 8):

  • If the initial velocity, v0, is oriented in the azimuthal direction, the probability of the velocity being exactly the one needed for a circular orbit is extremely small. This means that the resulting orbit will not be circular: if v0 > vcirc, the object will travel to larger R, and if v0 < vcirc, the object will fall inwards to smaller R. The initially azimuthal velocity will therefore probably reorient in the cylindrical coordinate system, obtain a radial component, and, because a Hobbs et al. (2005) kick is typically larger than vcirc, this radial component will most likely be positive.

  • If v0 is oriented in the negative radial direction, it will travel towards the Galactic centre and overshoot to the positive R direction (neglecting possible collisions with other objects). Because of the cylindrical coordinate system, a negative radial velocity will turn into a positive radial velocity, while the opposite is not true.

Because of these effects, the vectors of azimuthal and negative radial velocities are able to obtain a positive radial component. However, positive radial motion does not reorient itself in a similar way. Positive radial motion can either escape the Galaxy or decelerate enough to turn into negative radial motion.

To reproduce these effects in our simulation, we examined the orbits of objects that get an initial velocity entirely in one of these cylindrical dimensions. We put objects at R0 = 2.6 kpc – since this is the half-light radius of the disc, approximating the peak in Fig. 4 – and give them initial velocities of 100 km s−1, 400 km s−1, and 700 km s−1, in the azimuthal or radial direction. Figure 9 shows the orbits resulting from these initial conditions.

thumbnail Fig. 9.

Trajectories of objects starting at R0 = Rd = 2.6 kpc that get an initial velocity (v0) either completely in the ϕ direction (top row) or in the R direction (bottom row). For each direction, we show the trajectories of v0 = 100 km s−1 (light orange), 400 km s−1 (orange), and 700 km s−1 (dark orange), and add –400 km s−1 (dotted orange) and –700 km s−1 (dotted dark orange) for the radial direction. The first column shows the galactocentric velocity (v), and the second and third column show the R and z positions of the objects, respectively. In the fourth column, we show the azimuthal velocity (vϕ) and in the fifth column we show the radial velocity (vR). The vertical velocity (vz) does not obtain non-zero values in these orbits.

The figure shows that azimuthal velocity can indeed reorient into radial velocity. Within the first 20 Myr, the azimuthal velocity can decrease rapidly while the radial velocity increases, and this effect is greatest for the largest kick. The smallest v0 (i.e. 100 km s−1) is smaller than the local vcirc, causing the object to fall inwards (as sketched in Fig. 8). The radial velocity, in contrast, cannot reorient into other dimensions. However, as the simulated orbits show, a negative radial velocity can turn positive when crossing the Galactic centre, and a positive radial velocity can turn negative when it is decelerated completely and falls back to the Galactic centre. As Fig. 9 shows, the former tends to happen on a shorter timescale, causing an increase in vR before 20 Myr. Then, these objects are decelerated, where the object with the highest v0 (i.e. 700 km s−1) escape the Galaxy.

Kicked objects can of course also obtain a v0 in the z direction. Because of the gravitational attraction to the Galactic centre, these orbits are also able to obtain a positive vR component. However, for most kicks this will mostly occur at positions with a relatively large |z|, that is to say, outside of the solar neighbourhood. This effect, therefore, probably does not contribute significantly to the initial increase in vR, since this increase looks similar in the solar neighbourhood and the total population.

Moreover, if the kick vectors are oriented isotropically, the probability of the initial velocity being pointed entirely in one of the directions is negligible. We therefore do not argue that the discussed effects are a comprehensive description of the simulated orbits, but instead conclude they do show that the velocity vector initially tends to reorient in the positive radial direction, which can help explain the initial increase in vR in our simulation.

5. Pulsars

Since in our simulation kicked objects decelerate due to the Galactic potential, we can investigate whether any such signal is visible in observational data. To do this, we looked at pulsar data since they experience kicks and their age can be estimated. We converted the 3D velocities from our simulation into 2D proper motions, which we compared to observed proper motions of pulsars (Sect. 5.1). Then, we used a Monte Carlo estimation to determine the 3D galactocentric velocities of the pulsars, and compared them to our simulation (Sect. 5.2).

5.1. Proper motions

We used the pulsar data from the ATNF pulsar catalogue (Manchester et al. 2005), and selected all the pulsars with estimates for position, proper motion, distance, and age, except those with a binary companion. In total, the ATNF pulsar catalogue contains 280 pulsars that meet these requirements. We compared 2D proper motions of the solar neighbourhood population in our simulation with the pulsars that have estimated distances below 2 kpc, and compared the total population in our simulation with estimated 3D velocities of all 280 pulsars.

We used ASTROPY2 v.6.0.0 (Astropy Collaboration 2013, 2018, 2022) to convert the simulated objects in the solar neighbourhood from a galactocentric frame to a heliocentric frame and converted the 3D velocity vectors into 2D proper motion vectors. Since we needed to account for all 12 spheres that we used to determine the objects in the solar neighbourhood, we rotated the coordinates in our simulation by n ⋅ 2π/12 for each nth sphere and add the proper motions of the objects within 2 kpc of the Solar System to the distribution. We also corrected for the motion of the Sun by rotating the objects at time t by −t ⋅ vϕ, ⊙/R⊙ (using the values from Eqs. (3) and (4)).

For each pulsar, we sampled 103 values from Gaussians describing the observed proper motion and its uncertainty. After doing this for the proper motion in both right ascension (μα) and declination (μδ), we combined them to get the total 2D proper motion vector ( μ 2 = μ α 2 + μ δ 2 ) $ \left(\mu^2=\mu_{\alpha}^2+\mu_{\delta}^2\right) $. Then, we summed all posterior distributions of pulsars with a certain age to create the proper motion distributions of young and old pulsars.

There are many reasons why the pulsar data should not match our simulation. For example:

  • In our simulation, we assumed the kick velocities are distributed according to the Hobbs et al. (2005) distribution. However, whether this distribution accurately describes neutron star kicks is still a matter of debate, with more recent studies finding alternative distributions, some of which are bimodal (e.g. Verbunt et al. 2017; Igoshev et al. 2021; O’Doherty et al. 2023).

  • While our simulation has cylindrical symmetry, with the Galactic centre in the origin, the pulsar data have a bias towards nearby pulsars. After all, the closer a pulsar is to the Solar System, the lower the luminosity threshold for detectability (e.g. Chrimes et al. 2021).

  • The distance estimates used in the ATNF catalogue may be inaccurate (as argued by e.g. Verbunt et al. 2017), which would influence our estimation of the 3D velocity based on the 2D proper motion.

  • The pulsar ages stated in the ATNF catalogue (Manchester et al. 2005) are the characteristic spin-down ages, which equal P/(2Ṗ) for a period P, and assume that the braking index equals 3, and the initial period is significantly smaller than the current period (e.g. Ostriker & Gunn 1969; Shapiro & Teukolsky 1983). However, these assumptions may not be completely accurate for all pulsars (e.g. H.E.S.S. Collaboration 2011), meaning the true age of a pulsar probably differs from the characteristic age (Jiang et al. 2013).

  • The Monte Carlo estimates of the 3D pulsar velocities (as described in Sect. 5.2) assume an isotropically distributed velocity vector. However, our simulation shows that even after a few million years, the velocity vectors are far from isotropic (see e.g. Fig. 5).

  • There may be a spin-kick alignment (Johnston et al. 2005), which means that there is an observational bias towards pulsars moving radially with respect to us. Furthermore, if kicked objects preferentially end up on more radial orbits this contradicts the assumption of isotropy, causing an underestimation of the kick velocity (Mandel & Igoshev 2023).

However, despite all of these arguments for differences, we find that the pulsar data and our simulation match reasonably well.

Figure 10 shows the cumulative distributions of the proper motions, and comparing the observed pulsars younger than 20 Myr to the ones older than 20 Myr shows a deceleration similar to the one in our simulation. The proper motion distribution of points in our simulation below 20 Myr approximates the distribution of pulsars younger than 20 Myr. For older pulsars, there is only a small difference between the distribution of pulsars with ages in between 20 Myr and 200 Myr, and the distribution of all pulsars older than 20 Myr (which also includes pulsars older than 200 Myr); both of these distributions are well described by our simulation. This means that, despite all the observational biases and other reasons for differences between simulation and observation, they match relatively well.

thumbnail Fig. 10.

Comparison between pulsar data and our simulation. The pulsar data are taken from the ATNF catalogue (Manchester et al. 2005), as described in the text. The top-left panel shows the observed proper motions (μ) of pulsars within 2 kpc of the Solar System, for ages below 20 Myr (light red, consists of 72 pulsars), between 20 Myr and 200 Myr (red, consists of 22 pulsars), and above 20 Myr (dark red, consists of 49 pulsars), in cumulative histograms with a bin width of 10 mas yr−1. In the top-right panel, we show the proper motions of the simulated objects in the solar neighbourhood, for ages below 20 Myr (light blue) and above 20 Myr (dark blue), in a similar histogram. The bottom-left panel shows the sum of the Monte Carlo estimated posteriors of the pulsar galactocentric velocities (v), in histograms with a bin width of 20 km s−1, together with the Hobbs et al. (2005) kick distribution. These distributions also include pulsars at distances greater than 2 kpc from the Solar System, resulting in populations with ages below 20 Myr, between 20 Myr and 200 Myr, and above 200 Myr with sizes of 179, 45, and 101, respectively. The bottom-right panel shows similar histograms, for the total population in our simulation, as is also shown in Fig. 2. Similarly to Fig. 2, we evaluate our simulation every 1 Myr for the first 20 Myr, after which we evaluate the simulation every 10 Myr.

5.2. Galactocentric velocities

Besides the 2D proper motion comparison, we also investigated whether there is an agreement between the 3D galactocentric velocities of the pulsars and our simulation. However, only the 2D proper motion can be observed, meaning we only know the projection of the 3D velocity vector on the sky. To estimate the galactocentric velocity based on the proper motion, we used a Monte Carlo estimation. For each pulsar, we constructed three Gaussians that correspond to the mean and uncertainty of the proper motion (μα and μδ), and distance to the pulsar (D), as stated in the ATNF catalogue. Then, we sampled 103 values from these Gaussian distributions (where we require that D > 0), based on which we determine the transverse part of the 3D velocity vector. To determine the total 3D velocity, we assumed galactocentric isotropy and uniformly sampled a value u between −1 and 1, for each set of values of μα, μδ, and D. After correcting for the motion of the Solar System, we set the angle between the transverse velocity and the total velocity to be cos−1(u) (see Gaspari et al. 2024, for a more detailed description). We determined the magnitude of the 3D galactocentric velocity vector for all 103 samples, and thus a posterior velocity distribution for each pulsar. To compare the observations with our simulation, we summed the posteriors of pulsars with a certain age.

Figure 10 shows the results of our Monte Carlo estimation, including the posteriors for pulsars with ages below 20 Myr, in between 20 Myr and 200 Myr, and above 20 Myr (including pulsars older than 200 Myr). The young pulsars start with a velocity distribution relatively close to the Hobbs et al. (2005) kick distribution, though not as close as our simulation. Moreover, the older pulsars show a decelerated velocity distribution, which looks similar to the decelerated distribution of our simulation (for the total population). Both the posteriors of the observations and the distribution in our simulation show a peak around 200 km s−1.

The fact that our simulation resembles the pulsar data well (especially for older pulsars), despite all the reasons why it should not, is either coincidental or implies that the deceleration due to the Galactic potential is a major factor in shaping the velocity distributions of pulsars older than 20 Myr. Some aspects of the match between simulation and observation are accidental. For example, the high-velocity tail in the total population of our simulation consists of objects that get a high kick and escape the Galaxy, while the high-velocity tail in the pulsar data is more likely due to uncertainty in our Monte Carlo estimation. Nevertheless, we argue that the general match between the deceleration seen in the pulsars and the deceleration in our simulation, both in proper motion and galactocentric velocity, is unlikely to be coincidental, and shows how deceleration by the Galactic potential has a major influence on the pulsar velocities. After all, we find that the pulsar velocities are well described by the velocity distribution of the total population, meaning observational biases or selection effects do not influence the results significantly.

It is not surprising that the distributions of older pulsars match our simulation better than the distribution of young pulsars, since, as we mentioned in Sects. 3.1 and 4.1, we find that the velocity distribution of older objects in the solar neighbourhood remains constant when varying the kick distribution (Appendix B). In other words, the velocities of older pulsar may not depend significantly on the kick distribution of Hobbs et al. (2005) being accurate.

6. Conclusion

In this work we simulated the trajectories of objects that migrated through the Galaxy after receiving a velocity kick (Sect. 2) and find that they are decelerated by the Galactic potential (Sect. 3). This deceleration is a result of the kicks that increase the trajectories’ apocentres (Sect. 4), and is consistent with the deceleration observed in the velocities of Galactic pulsars at later ages (Sect. 5). Based on this, we conclude the following:

  • Objects that receive a kick tend to have their trajectories’ apocentres increased when they do not get unbound. As a result, they spend most of their time at larger radii, where they have lower velocities and are hence decelerated with respect to their initial velocity (Fig. 3). This causes the velocity distribution to shift towards lower values within 20–30 Myr after the kick, after which the median velocity stabilises at around 200 km s−1 (Fig. 2).

  • In the solar neighbourhood, objects older than about 20 Myr are, approximately, a representative sample of the total population, both in terms of velocity (Fig. 2) and initial galactocentric radius (Fig. 4).

  • Objects that receive a small kick (e.g. 100 km s−1) are only slightly disturbed in their initial circular orbit, while higher kicks (e.g. 400 km s−1 or higher) significantly disturb the initial orbit and result in the deceleration described above. Therefore, objects that received small kicks keep travelling at around the same initial speed, while the slowest objects result from the upper tail of the kick distribution (Fig. 7).

  • Within the first 20 Myr, vϕ and vz decrease, while vR increases (Fig. 5). This is due to the reorientation of the initial velocity vector into the positive radial direction (Fig. 8), and suggests that nearby objects have non-isotropic velocities within this period.

  • Objects in the solar neighbourhood with identical velocities and ages can have vastly different origins in terms of vkick and R0 (Fig. 7).

  • Using the Galactic pulsars from the ATNF catalogue (Manchester et al. 2005), we find that both the observed 2D proper motions and the estimated 3D galactocentric velocities of the pulsars resemble our simulation relatively well (Fig. 10). This indicates that the deceleration due to migration through the Galactic potential plays a major role in determining the pulsar velocities.

As a more general conclusion, we argue that the observational bias posed by Cordes (1986) and Lyne & Lorimer (1994) – the hypothesis that high-velocity objects leaving the solar neighbourhood sooner causes a selection effect – is insufficient in explaining the apparent deceleration of pulsars as they age. Instead, we find that the difference between the solar neighbourhood and the total population of kicked objects is merely the high-velocity tail, and the deceleration can be found in all populations: it produces a peak at approximately 200 km s−1. Therefore, we argue that the deceleration of pulsars as they age is not a selection effect but a genuine deceleration, caused by the Galactic potential.

It would be interesting to further investigate the problem of how to determine the kick velocity of an object older than 20 Myr based on its observed proper motion and position. After all, we find that the (scalar) velocity distribution of objects in the solar neighbourhood is mostly insensitive to the kick distribution (Appendix B). Since our simulation is a relatively simplistic one, it can be compared to any population of kicked objects, such as binaries and – as we discuss in this work – pulsars. Such a comparison could potentially help explain the observed velocities of these binaries (e.g. Zhao et al. 2023) and pulsars (e.g. Cordes & Chernoff 1998). Moreover, it could be interesting to repeat the Monte Carlo estimation of the pulsar velocities but using the distributions of vϕ, vR, and vz we show in Fig. 5 instead of assuming isotropy.

Acknowledgments

We thank Frank Verbunt, Gijs Nelemans, and Ilya Mandel for useful discussions regarding this Master’s project, and the referee, Andrei Igoshev, for valuable comments that helped improve this paper. A.J.L. was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 725246), and N.G. acknowledges studentship support from the Dutch Research Council (NWO) under the project number 680.92.18.02. In this work, we used pulsar data from the ATNF catalogue (Manchester et al. 2005, available at https://www.atnf.csiro.au/research/pulsar/psrcat). Moreover, we made use of NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), GALPY (Bovy 2015), and ASTROPY, a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

References

  1. Andrews, J. J., & Kalogera, V. 2022, ApJ, 930, 159 [NASA ADS] [CrossRef] [Google Scholar]
  2. Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 2002, ApJ, 568, 289 [NASA ADS] [CrossRef] [Google Scholar]
  3. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  6. Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116 [NASA ADS] [CrossRef] [Google Scholar]
  7. Belczynski, K., Perna, R., Bulik, T., et al. 2006, ApJ, 648, 1110 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beniamini, P., & Piran, T. 2016, MNRAS, 456, 4089 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [Google Scholar]
  10. Bhattacharya, D., Wijers, R. A. M. J., Hartman, J. W., & Verbunt, F. 1992, A&A, 254, 198 [NASA ADS] [Google Scholar]
  11. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, N.J.: Princeton University Press) [Google Scholar]
  12. Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265 [NASA ADS] [Google Scholar]
  13. Blaes, O., & Madau, P. 1993, ApJ, 403, 690 [NASA ADS] [CrossRef] [Google Scholar]
  14. Blaes, O., & Rajagopal, M. 1991, ApJ, 381, 210 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  16. Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
  17. Brandt, N., & Podsiadlowski, P. 1995, MNRAS, 274, 461 [NASA ADS] [CrossRef] [Google Scholar]
  18. Burrows, A., Wang, T., Vartanyan, D., & Coleman, M. S. B. 2024, ApJ, 963, 63 [NASA ADS] [CrossRef] [Google Scholar]
  19. Chrimes, A. A., Levan, A. J., Groot, P. J., Lyman, J. D., & Nelemans, G. 2021, MNRAS, 508, 1929 [NASA ADS] [CrossRef] [Google Scholar]
  20. Church, R. P., Levan, A. J., Davies, M. B., & Tanvir, N. 2011, MNRAS, 413, 2004 [NASA ADS] [CrossRef] [Google Scholar]
  21. Coleman, M. S. B., & Burrows, A. 2022, MNRAS, 517, 3938 [NASA ADS] [CrossRef] [Google Scholar]
  22. Cordes, J. M. 1986, ApJ, 311, 183 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cordes, J. M., & Chernoff, D. F. 1998, ApJ, 505, 315 [NASA ADS] [CrossRef] [Google Scholar]
  24. Dewey, R. J., & Cordes, J. M. 1987, ApJ, 321, 780 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [Google Scholar]
  26. Drimmel, R., & Poggio, E. 2018, Res. Notes Am. Astron. Soc., 2, 210 [Google Scholar]
  27. Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332 [CrossRef] [Google Scholar]
  28. Feltzing, S., & Bensby, T. 2008, Phys. Scr. Vol. T, 133, 014031 [CrossRef] [Google Scholar]
  29. Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [Google Scholar]
  30. Fortin, F., García, F., Chaty, S., Chassande-Mottin, E., & Simaz Bunzel, A. 2022, A&A, 665, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Fragos, T., Willems, B., Kalogera, V., et al. 2009, ApJ, 697, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fryer, C., & Kalogera, V. 1997, ApJ, 489, 244 [NASA ADS] [CrossRef] [Google Scholar]
  33. Fryer, C. L., Woosley, S. E., & Hartmann, D. H. 1999, ApJ, 526, 152 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gandhi, P., Rao, A., Johnson, M. A. C., Paice, J. A., & Maccarone, T. J. 2019, MNRAS, 485, 2642 [Google Scholar]
  35. Gaspari, N., Levan, A. J., Chrimes, A. A., & Nelemans, G. 2024, MNRAS, 527, 1101 [Google Scholar]
  36. Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011 [Google Scholar]
  37. Grady, J., Belokurov, V., & Evans, N. W. 2020, MNRAS, 492, 3128 [NASA ADS] [CrossRef] [Google Scholar]
  38. GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Gunn, J. E., & Ostriker, J. P. 1970, ApJ, 160, 979 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569 [NASA ADS] [Google Scholar]
  41. Harris, C. R., Millman, K. J., Van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hartman, J. W., Bhattacharya, D., Wijers, R., & Verbunt, F. 1997, A&A, 322, 477 [NASA ADS] [Google Scholar]
  43. Hartmann, D., Epstein, R. I., & Woosley, S. E. 1990, ApJ, 348, 625 [NASA ADS] [CrossRef] [Google Scholar]
  44. H.E.S.S. Collaboration (Abramowski, A., et al.) 2011, A&A, 533, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Hills, J. G. 1983, ApJ, 267, 322 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  48. Iben, I., Jr., & Tutukov, A. V. 1996, ApJ, 456, 738 [NASA ADS] [CrossRef] [Google Scholar]
  49. Igoshev, A. P. 2020, MNRAS, 494, 3663 [NASA ADS] [CrossRef] [Google Scholar]
  50. Igoshev, A. P., Chruslinska, M., Dorozsmai, A., & Toonen, S. 2021, MNRAS, 508, 3345 [NASA ADS] [CrossRef] [Google Scholar]
  51. Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426 [NASA ADS] [CrossRef] [Google Scholar]
  52. Janka, H.-T. 2012, Annu. Rev. Nucl. Part. Sci., 62, 407 [NASA ADS] [CrossRef] [Google Scholar]
  53. Janka, H.-T. 2013, MNRAS, 434, 1355 [Google Scholar]
  54. Jiang, L., Zhang, C.-M., Tanni, A., & Zhao, H.-H. 2013, Int. J. Mod. Phys. Conf. Ser., 23, 95 [NASA ADS] [CrossRef] [Google Scholar]
  55. Jiang, Z., Wang, J., Zhang, F., et al. 2020, MNRAS, 498, 926 [CrossRef] [Google Scholar]
  56. Johnston, S., Hobbs, G., Vigeland, S., et al. 2005, MNRAS, 364, 1397 [NASA ADS] [CrossRef] [Google Scholar]
  57. Jonker, P. G., & Nelemans, G. 2004, MNRAS, 354, 355 [Google Scholar]
  58. Kalogera, V. 1996, ApJ, 471, 352 [Google Scholar]
  59. Kiel, P. D., & Hurley, J. R. 2009, MNRAS, 395, 2326 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [CrossRef] [Google Scholar]
  61. Lyne, A. G., & Lorimer, D. R. 1994, Nature, 369, 127 [NASA ADS] [CrossRef] [Google Scholar]
  62. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 [Google Scholar]
  63. Mandel, I. 2016, MNRAS, 456, 578 [Google Scholar]
  64. Mandel, I., & Igoshev, A. P. 2023, ApJ, 944, 153 [NASA ADS] [CrossRef] [Google Scholar]
  65. Mandhai, S., Lamb, G. P., Tanvir, N. R., et al. 2022, MNRAS, 514, 2716 [NASA ADS] [CrossRef] [Google Scholar]
  66. McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
  67. O’Doherty, T. N., Bahramian, A., Miller-Jones, J. C. A., et al. 2023, MNRAS, 521, 2504 [CrossRef] [Google Scholar]
  68. Ofek, E. O. 2009, PASP, 121, 814 [NASA ADS] [CrossRef] [Google Scholar]
  69. Ostriker, J. P., & Gunn, J. E. 1969, ApJ, 157, 1395 [Google Scholar]
  70. Paczynski, B. 1990, ApJ, 348, 485 [NASA ADS] [CrossRef] [Google Scholar]
  71. Perna, R., & Belczynski, K. 2002, ApJ, 570, 252 [NASA ADS] [CrossRef] [Google Scholar]
  72. Podsiadlowski, P., Langer, N., Poelarends, A. J. T., et al. 2004, ApJ, 612, 1044 [NASA ADS] [CrossRef] [Google Scholar]
  73. Popov, S. B., Colpi, M., Treves, A., et al. 2000, ApJ, 530, 896 [NASA ADS] [CrossRef] [Google Scholar]
  74. Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872 [Google Scholar]
  75. Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77 [Google Scholar]
  76. Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131 [Google Scholar]
  77. Repetto, S., Davies, M. B., & Sigurdsson, S. 2012, MNRAS, 425, 2799 [Google Scholar]
  78. Sartore, N., Ripamonti, E., Treves, A., & Turolla, R. 2010, A&A, 510, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Shapiro, S. L., & Teukolsky, S. A. 1983, Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects (New York: Wiley) [CrossRef] [Google Scholar]
  80. Sweeney, D., Tuthill, P., Sharma, S., & Hirai, R. 2022, MNRAS, 516, 4971 [NASA ADS] [CrossRef] [Google Scholar]
  81. Treves, A., Turolla, R., Zane, S., & Colpi, M. 2000, PASP, 112, 297 [NASA ADS] [CrossRef] [Google Scholar]
  82. Van den Heuvel, E. P. J., 2007, in The Multicolored Landscape of Compact Objects and Their Explosive Origins, eds. T. di Salvo, G. L. Israel, L. Piersant, et al., AIP Conf. Ser., 924, 598 [Google Scholar]
  83. Van den Heuvel, E. P. J., & Van Paradijs, J. 1997, ApJ, 483, 399 [CrossRef] [Google Scholar]
  84. Van den Heuvel, E. P. J., Portegies Zwart, S. F., Bhattacharya, D., & Kaper, L. 2000, A&A, 364, 563 [NASA ADS] [Google Scholar]
  85. Van Paradijs, J., & White, N. 1995, ApJ, 447, L33 [NASA ADS] [CrossRef] [Google Scholar]
  86. Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MNRAS, 481, 4009 [Google Scholar]
  88. Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169 [Google Scholar]
  89. Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874 [Google Scholar]
  90. Willems, B., Kalogera, V., & Henninger, M. 2004, ApJ, 616, 414 [NASA ADS] [CrossRef] [Google Scholar]
  91. Wongwathanarat, A., Janka, H. T., & Müller, E. 2013, A&A, 552, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Zane, S., Turolla, R., Zampieri, L., Colpi, M., & Treves, A. 1995, ApJ, 451, 739 [NASA ADS] [CrossRef] [Google Scholar]
  93. Zevin, M., Kelley, L. Z., Nugent, A., et al. 2020, ApJ, 904, 190 [NASA ADS] [CrossRef] [Google Scholar]
  94. Zhang, Z., Gilfanov, M., & Bogdán, Á. 2013, A&A, 556, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Zhao, Y., Gandhi, P., Dashwood Brown, C., et al. 2023, MNRAS, 525, 1498 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Velocity evolution

To examine the deceleration in the velocity distribution more closely, we show the velocities of objects in the solar neighbourhood in Fig. A.1, similarly to Fig. 5 but on a logarithmic axis that is extended to 1 Gyr. The figure shows how the mean of the scalar velocity starts to decrease rapidly after about 10 Myr, and stays at approximately 200 km/s. This shows how, despite the fact that the system may not have reached full dynamical equilibrium yet (e.g. Binney & Tremaine 1987, ch. 4), the deceleration to mean velocities of ∼200 km/s in the solar neighbourhood appears to be stable beyond 200 Myr. This is true for the evolution in the cylindrical velocity components as well.

thumbnail Fig. A.1.

Velocity evolution of the objects in the solar neighbourhood, similar to Fig. 5, but on a logarithmic scale extended to 1 Gyr. The panels show the magnitude of the velocity vector (v), as well as the cylindrical velocity components (vϕ, vR, and vz). The first 25 Myr, the velocities are evaluated every 1 Myr, after which they are evaluated every 10 Myr, as shown in the 2D histogram with velocity bins of 30 km/s. The black line traces the mean velocity ( v ¯ $ \overline{v} $) at each point in time.

Appendix B: Sensitivity to the initial kick

In our simulation, we used the Hobbs et al. (2005) kick distribution, which is a Maxwellian with a σ of 265 km/s. However, we find that the peak at approximately 200 km/s (as shown in Fig. 2) also forms if we change the kick distribution. To show this, we ran our simulation for different kick distributions, all described by Maxwellians:

f ( v ) = 2 π v 2 σ 3 exp ( − v 2 2 σ 2 ) , $$ \begin{aligned} f(v)=\sqrt{\dfrac{2}{\pi }}\dfrac{v^2}{\sigma ^3}\exp \left(-\dfrac{v^2}{2\sigma ^2}\right)\quad , \end{aligned} $$(B.1)

with values for σ varying from 100 km/s to 500 km/s. Figure B.1 shows the results of varying the kick distribution in this way for the observed velocity distributions in the solar neighbourhood. As can be seen in the figure, varying the kick distribution does not affect the decelerated distribution meaningfully, since the distributions above 20 Myr are all nearly indistinguishable. This is not surprising, since a larger σ increases the fraction of the total population that get a high enough vkick to escape the Galaxy, and after all, these are not observed in the solar neighbourhood. We find, therefore, that it is extremely difficult to estimate the kick distribution based on the magnitude of the velocities of observed objects in the solar neighbourhood that are older than 20 Myr, although it might be possible to constrain kicks based on the complete velocity vector (e.g. Atri et al. 2019).

thumbnail Fig. B.1.

Our simulation, but for 103 point masses instead of 104, and for varying kick distributions. Similarly to the results shown in Fig. 2, we evaluate the points below 20 Myr every 1 Myr (light blue curve) and after 20 Myr every 10 Myr, which provides the velocity distributions shown in normalised histograms with a bin width of 30 km/s. For the kick distributions, we use Maxwellians (Eq. B.1) with a σ varying from 100 km/s to 500 km/s, with intervals of 50 km/s.

All Figures

thumbnail Fig. 1.

Initial positions of point masses in our simulation (blue dots), plotted over the synthetic MW image, based on the works of Gaspari et al. (2024) and Chrimes et al. (2021), as described in the main text. The figure contains 500 of the 104 simulated point masses. The Sun symbol represents the Solar System, at R⊙ = 8.122 kpc and z⊙ = 20.8 pc.
In the text
thumbnail Fig. 2.

Distributions of galactocentric velocities (v = |v|, top row) and peculiar velocities in the local standard of rest of the point mass (vLSR, bottom row) for the simulated point masses. For the first 20 Myr, we evaluated the velocities every 1 Myr and show the total, normalised distribution in histograms with a bin width of 20 km s−1 (light blue curves). After 20 Myr, we evaluated the velocities every 10 Myr and also show the total, normalised distribution in similar histograms (dark blue curves). The dashed black curves show the kick distribution (i.e. a Maxwellian with σ = 265 km s−1, as used by Hobbs et al. 2005). The three panels show the evaluated velocities for the three populations: the solar neighbourhood (left), the thick disc (centre), and the total population (right). These populations have sizes of 24 817, 144 496, and 210 000 for ages below 20 Myr, and 7079, 54 320, and 180 000 for ages above 20 Myr, respectively.
In the text
thumbnail Fig. 3.

Initial conditions versus decelerated conditions, for the total population of our simulation. The left panel shows the initial velocities of the objects (v0) versus the decelerated velocities (vdec), which we define as any velocity for t > 20 Myr. We evaluated the decreased velocities every 10 Myr and show the resulting distribution in a 2D histogram with bins of 20 km s−1. The right panel shows the initial galactocentric radial position of the objects (R0) versus the galactocentric radial position of the objects for t > 20 Myr (Rdec), evaluated every 10 Myr and shown in a 2D histogram with bins of 0.5 kpc. Here, 83% of the objects are located at Rdec > R0.
In the text
thumbnail Fig. 4.

Initial galactocentric radii (R0) of the objects that are in the solar neighbourhood below 20 Myr (light blue) and above 20 Myr (dark blue), shown in normalised histograms with bins of 1 kpc. The dotted curve shows the initial radii of the total population, in a similar histogram. The distribution of the total population peaks slightly above the half-light radius, Rd = 2.6 kpc (Eq. (1)), because the spiral arms shift the initial positions to higher radii compared to the disc.
In the text
thumbnail Fig. 5.

Evolution of the galactocentric velocity (v) and the cylindrical components of the velocity vector (vϕ, vR, and vz). The figure shows normalised distributions, in 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1, 20 km s−1, and 10 km s−1 for the solar neighbourhood, the thick disc, and the total population, respectively. For clarity, we set the maximum values of the distributions at 4 × 10−5 s km−1 Myr−1 ≈ 13 × 10−19 km−1, even though the peaks in the last row exceed this maximum to a certain degree. The black lines show the mean velocities at each point in time.
In the text
thumbnail Fig. 6.

Grid of simulations of 103 kicked objects, with a fixed initial galactocentric radius, R0 (rows: 1 kpc, 5 kpc, and 9 kpc), and a single value for vkick (columns: 100 km s−1, 400 km s−1, and 700 km s−1). The figure shows the evolution of the galactocentric velocity (v) of the total population, in 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1. Even though these simulations use a single value for vkick, the velocity distributions do not start as delta functions, because of the vcirc term in Eq. (2).
In the text
thumbnail Fig. 7.

Grid of simulations of 103 kicked objects, similar to Fig. 6, but showing only the objects that are in the solar neighbourhood. Again, these simulations have a fixed initial galactocentric radius, R0 (rows: 1 kpc, 5 kpc, and 9 kpc), and a single value for vkick (columns: 100 km s−1, 400 km s−1, and 700 km s−1). Similarly to Fig. 6, the evolution of the galactocentric velocity (v) is shown through 2D histograms with time bins of 10 Myr and velocity bins of 30 km s−1. However, since the sizes of the populations in the panels differ significantly, they are shown on different colour scales.
In the text
thumbnail Fig. 8.

Schematic diagrams of hypothetical paths (light orange curves) of objects (orange dots) resulting from an initial velocity (orange arrows) completely in the azimuthal (left panel) or negative radial (central panel) direction. The shown initial velocity vectors and paths are merely qualitative estimates of possible paths that potentially contribute to the initial acceleration of vR (as found in Fig. 5). In Fig. 9, we show examples of orbits that are a result of initial velocities in these directions.
In the text
thumbnail Fig. 9.

Trajectories of objects starting at R0 = Rd = 2.6 kpc that get an initial velocity (v0) either completely in the ϕ direction (top row) or in the R direction (bottom row). For each direction, we show the trajectories of v0 = 100 km s−1 (light orange), 400 km s−1 (orange), and 700 km s−1 (dark orange), and add –400 km s−1 (dotted orange) and –700 km s−1 (dotted dark orange) for the radial direction. The first column shows the galactocentric velocity (v), and the second and third column show the R and z positions of the objects, respectively. In the fourth column, we show the azimuthal velocity (vϕ) and in the fifth column we show the radial velocity (vR). The vertical velocity (vz) does not obtain non-zero values in these orbits.
In the text
thumbnail Fig. 10.

Comparison between pulsar data and our simulation. The pulsar data are taken from the ATNF catalogue (Manchester et al. 2005), as described in the text. The top-left panel shows the observed proper motions (μ) of pulsars within 2 kpc of the Solar System, for ages below 20 Myr (light red, consists of 72 pulsars), between 20 Myr and 200 Myr (red, consists of 22 pulsars), and above 20 Myr (dark red, consists of 49 pulsars), in cumulative histograms with a bin width of 10 mas yr−1. In the top-right panel, we show the proper motions of the simulated objects in the solar neighbourhood, for ages below 20 Myr (light blue) and above 20 Myr (dark blue), in a similar histogram. The bottom-left panel shows the sum of the Monte Carlo estimated posteriors of the pulsar galactocentric velocities (v), in histograms with a bin width of 20 km s−1, together with the Hobbs et al. (2005) kick distribution. These distributions also include pulsars at distances greater than 2 kpc from the Solar System, resulting in populations with ages below 20 Myr, between 20 Myr and 200 Myr, and above 200 Myr with sizes of 179, 45, and 101, respectively. The bottom-right panel shows similar histograms, for the total population in our simulation, as is also shown in Fig. 2. Similarly to Fig. 2, we evaluate our simulation every 1 Myr for the first 20 Myr, after which we evaluate the simulation every 10 Myr.
In the text
thumbnail Fig. A.1.

Velocity evolution of the objects in the solar neighbourhood, similar to Fig. 5, but on a logarithmic scale extended to 1 Gyr. The panels show the magnitude of the velocity vector (v), as well as the cylindrical velocity components (vϕ, vR, and vz). The first 25 Myr, the velocities are evaluated every 1 Myr, after which they are evaluated every 10 Myr, as shown in the 2D histogram with velocity bins of 30 km/s. The black line traces the mean velocity ( v ¯ $ \overline{v} $) at each point in time.
In the text
thumbnail Fig. B.1.

Our simulation, but for 103 point masses instead of 104, and for varying kick distributions. Similarly to the results shown in Fig. 2, we evaluate the points below 20 Myr every 1 Myr (light blue curve) and after 20 Myr every 10 Myr, which provides the velocity distributions shown in normalised histograms with a bin width of 30 km/s. For the kick distributions, we use Maxwellians (Eq. B.1) with a σ varying from 100 km/s to 500 km/s, with intervals of 50 km/s.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.