© The Authors 2024

1 Introduction

Neutron star binaries – and in particular their mergers – can provide insight into a variety of astrophysical phenomena. For example, mergers of these binary neutron stars (BNSs), consisting of two neutron stars orbiting each other, have produced gravitational waves that have been detected (e.g. Abbott et al. 2017b) and are used to gain insight into the structure of neutron stars (e.g. Radice et al. 2018). Moreover, BNS mergers are thought to produce short-duration gamma-ray bursts (SGRBs; Eichler et al. 1989; Narayan et al. 1992), which is why their rates and locations have been topics of research (for a review see Berger 2014). Furthermore, the mergers of such binaries are a prime source of heavy element (r–process) enrichment throughout the Universe and so far are the only site for which direct evidence of heavy element production has been found (Berger et al. 2013; Tanvir et al. 2013; Arcavi et al. 2017; Kasen et al. 2017; Metzger 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Villar et al. 2017; Rastinejad et al. 2022; Troja et al. 2022; Levan et al. 2023, for a review see Nakar 2020). Understanding both the details of the mergers (e.g. component masses, ejecta masses, composition) and the locations of the mergers relative to host galaxies is then critical to determining their role in cosmic chemical evolution, which is why we are interested in estimating the systemic kicks of BNS systems.

After all, the merger locations of BNS with respect to their host galaxies are determined by the velocity kicks these systems receive on formation (e.g. Fryer et al. 1999; Bloom et al. 1999; Bulik et al. 1999; Belczynski et al. 2002, 2006; Perna & Belczynski 2002; Voss & Tauris 2003; Church et al. 2011; Abbott et al. 2017a; Zevin et al. 2020). That is, both neutron stars are formed in a supernova, and because of this, they receive a natal kick, whose magnitude has been estimated by analysing the velocities of young pulsars (Lyne & Lorimer 1994; Hobbs et al. 2005; Faucher-Giguère & Kaspi 2006; Arzoumanian et al. 2002; Verbunt et al. 2017; Igoshev 2020; Igoshev et al. 2021). The total systemic kick the binary receives is then a combination of these two natal kicks and a Blaauw kick due to mass loss at each supernova (Blaauw 1961; Hills 1983; Van den Heuvel et al. 2000; Zhang et al. 2013; Tauris et al. 2017; Andrews & Zezas 2019; Zhao et al. 2023). The magnitude of the BNSs’ systemic kicks has been debated. For instance, high kicks can explain SGRBs that appear to be host-less (Berger 2010; Fong & Berger 2013; Tunnicliffe et al. 2014; O’Connor et al. 2022; Fong et al. 2022). In particular, while generically large kicks like those seen in pulsars can re-create the observed offsets to short GRBs reasonably well (Church et al. 2011), Behroozi et al. (2014) find that the galaxy-offset distribution of SGRBs, as found by Fong & Berger (2013), can be explained by a SGRB progenitor kick distribution where only 19% of the objects are kicked with velocities >150 km/s. Indeed, large kicks may not be necessary to explain SGRBs (Perets & Beniamini 2021), and Galactic BNSs also appear to have experienced small kicks, as suggested by their orbital parameters (Beniamini & Piran 2016) and their small peculiar velocities (Gaspari et al. 2024).

We are interested in investigating whether we can constrain the systemic kicks of the Galactic BNSs, for instance, by analysing their current speeds. Disberg et al. (2024) have shown, however, that kicked objects near the Solar System that are older than a few tens of millions of years will obtain galac-tocentric speeds that are not representative of their kicks, due to their motion within the Galactic potential. For this reason we are interested in estimating the ages of the Galactic BNSs, and we aim to give an overview of their spin-down ages based on braking in rotational speed of the neutron star in the binary that is observed as pulsar and their kinematic ages based on their distance from the Galactic plane. For some of the Galactic BNSs, there are already estimates of the spin-down ages (e.g. Lorimer et al. 2005; Kargaltsev et al. 2006; Andrews & Mandel 2019) and the kinematic ages (e.g. Arzoumanian et al. 1999; Wex et al. 2000; Willems et al. 2004, 2006), but we aim to give a complete overview for the Galactic BNS systems that have a known location, proper motion, and no association with a globular cluster (as e.g. listed by Ding et al. 2024).

In order to investigate the kicks of the BNSs, we analysed the magnitude of their current galactocentric velocities (Monte Carlo estimated following Gaspari et al. 2024), taking into account the estimates for their ages and the deceleration found by Disberg et al. (2024). However, our estimates contain the complete velocity vectors, not just their magnitude, which provides additional information that can potentially be used to constrain their kicks. Atri et al. (2019), for instance, used estimated velocity vectors to trace back trajectories of black hole X-ray binaries to the Galactic plane, where they are thought to be formed, and determined the peculiar velocity at the disc crossing as a potential kick velocity (cf. O’Doherty et al. 2023, who applied this method to binaries consisting of a neutron star with a low-mass companion).

Likewise, we used the estimated velocity vectors in order to trace back the trajectories of the Galactic BNSs, but we are interested in using the shape of the complete trajectory in order to estimate their kicks. That is, Disberg et al. (2024) find that the initially circular Galactic orbits of kicked objects are disturbed by the kick, resulting in more eccentric Galactic orbits (as also noted by e.g. Hoang et al. 2022). We are interested in investigating whether we can use this fact to constrain kicks based on the eccentricity of the Galactic orbits of the BNS systems.

This paper is structured as follows. In Sect. 2, we investigate the ages of the Galactic BNSs. Then, in Sect. 3, we analyse the current speeds of the BNSs and use the work of Disberg et al. (2024) to show that these harbour little information about their kicks. As an alternative method to constrain kicks, we show in Sect. 4 that the eccentricity of a Galactic orbit is a better indication of kick velocity, and we estimate these eccentricities for the Galactic BNSs. In Sect. 5, we determine the relationship between kick velocity and eccentricity of the Galactic orbit, and we use this relationship to kinematically constrain the kicks of the BNSs. Finally, in Sect. 6, we summarise our findings and their implications.

2 Ages

We are interested in investigating the ages of the Galactic BNSs, since a young population of kicked objects might look different from an old population of kicked objects, due to their motion through the Galaxy (Disberg et al. 2024). In order to do this, we used two quantities that can tell us something about their ages: spin-down ages (Sect. 2.1) and kinematic ages (Sect. 2.2).

2.1 Spin-down ages

The first quantity related to the ages of the BNS systems uses the spin properties of the observed pulsars in the binaries. During their formation pulsars get spun-up to high frequencies, corresponding to an initial period P₀. The spin-down time it takes for the pulsar to slow down to a period P, given a certain braking-index n, is then given by $${\tau }_{\text{sd}}=\frac{P}{(n-1)\dot{P}}\left[1-{\left(\frac{P_0}{P}\right)}^{n-1}\right].$$(1)

The characteristic age of the pulsar (τ_c) is determined by assuming the braking is caused by pure dipole radiation, meaning n = 3 – while multipole radiation results in values of n > 3, which could reach values op to n = 5 if the braking is dominated by gravitational radiation (e.g. Camilo et al. 1994; Ferrario & Wickramasinghe 2007) - and the initial period is small compared to the current period. This gives, as a function of P and P or as a function of the frequency v = P⁻¹ and its derivative $\dot{v}=-\dot{P}{P}^{-2}$ (e.g. Shapiro & Teukolsky 1983; Arzoumanian et al. 1999): $${\tau }_{\text{c}}=\frac{P}{2\dot{P}}=-\frac{v}{2\dot{v}}.$$(2)

For some pulsars, τ_c seems to approximate their true age, for example because it aligns with an “expansion age” for a supernova remnant based on ejecta velocities (e.g. Wyckoff & Murray 1977), though it is not always a good estimate of the true age of a pulsar (Jiang et al. 2013; Zhang et al. 2022). For a millisecond pulsar (MSP) in a BNS, for instance, the assumptions that P₀ << P and n = 3 may be inaccurate. These MSPs are thought to have been spun-up by accreting mass from their companion, through which they were “recycled” to millisecond periods (e.g. Bhattacharya & Van den Heuvel 1991).

In order to estimate τ_c, we constructed Gaussian distributions corresponding to the observed values of P and Ṗ or v and $\dot{v}$ and their uncertainties (as given in Table 1). Then, we sampled these distributions 10³ times for each BNS and for each instance determined τ_c through Eq. (2). The resulting distributions of the characteristic ages are shown in Fig. 1, and their median values and uncertainties are given in Table 2. The distributions in Fig. 1 show that τ_c is well constrained, due to the precision of P and Ṗ measurements, and does not differ from values that can be found in literature (e.g. Kargaltsev et al. 2006; Andrews & Mandel 2019). However, it remains difficult to determine the true ages of these pulsars based on a spin-down age, in particular due to the (binary) evolutionary history of the MSPs. This is, for instance, reflected by the fact that the median τ_c of J1518+4904 exceeds a Hubble time. Nevertheless, Maoz & Nakar (2024) argue that τ_c is a “reliable indicator of age,” and that corrections for binary evolution are expected to be relatively small (based on Kiziltan & Thorsett 2010).

2.2 Kinematic ages

As an additional way to investigate the ages of the Galactic BNSs, we estimated their kinematic ages. That is, a BNS should be formed through supernovae in the thin disc, at ɀ ≈ 0 kpc, after which it is kicked into an orbit containing heights |ɀ| > 0. We can use the motion of the binaries away from or towards the disc to determine the last time they crossed ɀ = 0. This could potentially give a lower bound for their true age, since the binaries may have crossed the disc multiple times since their birth. However, even if a BNS has not crossed the disc multiple times in its kinematic history, there are many reasons why its kinematic age might differ (significantly) from its actual age. After all, a BNS receives two systemic kicks, after the SNe of each component, which means that it could already start to migrate through the Galaxy before receiving its final systemic kick. In addition, the kinematic ages are made uncertain by the unknown radial velocity with respect to the Solar System (v_r) and the unknown birth location of the binaries (which could be at small but non-zero |ɀ|). For these reasons we do not consider the kinematic ages of the BNSs as accurate estimates of their actual ages, but since there are several studies concerning the kinematic ages of individual BNSs, we determined them here to give a complete overview.

We estimated the galactocentric speeds for the Galactic BNSs in our sample. For these binaries, the sky locations, proper notions, and distances are given in Table 3. In order to determine the magnitude of their galactocentric velocity vector, we estimated the unknown velocity in the radial direction (v_r) through two different isotropy assumptions (following the method of Gaspari et al. 2024). Firstly, we assumed that the full galacto-centric velocity is oriented isotropically. We refer to this kind of isotropy as “GC isotropy.” For a BNS where the proper motion and distance translates into a transverse velocity v_t, the GC-isotropic estimate of v_r equals $$v_r=v_t\cot \theta ,$$(3)

where θ = arccos u and u is uniformly sampled between 0 and 1. We sampled u 10³ times, with which we created a distribution of v_r for each BNS. Secondly, we assumed that the peculiar velocity of the BNSs is oriented isotropically in the local standard of rest (LSR), which we refer to as “LSR isotropy.” If a BNS is located at a position where the velocity vector of the LSR equals u_LSR, then the unknown v_r is estimated by subtracting the 2D transverse part of the LSR velocity vector (v_LSR,t) from the BNS’ v_t and adding the radial part of the LSR motion: $$v_r=\Vert v_t-v_{\text{LSR},t}\Vert \cot \theta +v_{\text{LSR},r},$$(4)

where θ is defined similar to Eq. (3), meaning this also resulted in a distribution of v_r containing 10³ values. The LSR isotropy assumption is appropriate if the BNS receives low kicks, and its velocity is therefore dominated by the initial circular velocity of the LSR (as determined by the Galactic potential).

We obtained a v_r distribution from assuming either GC isotropy or LSR isotropy, and then used the BNS sky-locations, distances, and 2D proper motions of the BNS systems (which are given in Table 3) to construct Gaussian distributions based on their mean and uncertainties, for each binary. We sampled these distributions 10³ times and determined the 3D velocity vectors for each instance. Having obtained 10³ sets of positions and velocities for each BNS, we used the GALPY¹ v .1.9.0 package (Bovy 2015) to integrate the orbits backwards in the Milky Way (MW) potential of McMillan (2017), assuming that the circular velocity at the position of the Sun equals its azimuthal velocity (i.e. v_LSR(8.122 kpc) = 245.6 km/s, GRAVITY Collaboration 2018). We evaluated these orbits every Myr, up to 1 Gyr, in order to determine the points at which they cross ɀ = 0. If ɀ(t_i) • ɀ(t_i+1) < 0, meaning one has a positive ɀ and the other a negative one, we set the time of crossing as the average between these consecutive timestamps. Using this method, we determined the distributions of τ_kin for the first (χ₁) and second (χ₂) time each BNS crosses the disc (when tracing back their kinematic history), where we neglect orbits that do not cross within 1 Gyr, or do so at galactocentric radii R > 30 kpc.

In Fig. 1, we show the distributions of τ_kin for χ₁ and χ₂, assuming either GC isotropy or LSR isotropy. For most BNSs, except perhaps J0737−3039A/B and J1518+4904, the distributions of τ_kin for GC isotropy and LSR isotropy agree relatively well. For J0737−3039A/B, the characteristic age overlaps with τ_kin for χ₁ assuming GC isotropy, and χ₂ assuming LSR isotropy. In general, the characteristic ages of the binaries exceed their kinematic ages significantly. While it is difficult to estimate the reliability of the spin-down ages, the difference of one or more orders of magnitude between the spin-down and the kinematic ages does imply they were likely formed at χ_≥3. In Table 2 we list the median kinematic ages and their uncertainties. Also, the boxplots in Fig. 1 hide structures such as potential bimodality in the τ_kin distributions, which is why we show the disc crossings in Appendix A. For most binaries, and particularly for LSR isotropy, the τkin estimates describe the disc crossings relatively well.

Willems et al. (2004) discussed the kinematic age of J0737− 3039A/B, depending on assumed values of v_r and the angle of the 1D proper motion on the sky (Ω). They conclude that for almost all values of v_r and Ω, τ_kin for χ₁ is less than 100 Myr and τ_kin for χ₂ exceeds 20 Myr, both of which match our findings. In addition, Lorimer et al. (2005) employed spin-down models to constrain the age of J0747−3039A/B to 30–70 Myr, which also coincides with the region of overlap between τ_c and τ_kin we find. Willems et al. (2006) found that certain values of v_r give rise to τ_kin ~ 5 Myr for χ₁. These solutions can also be seen in Appendix A. For B1534+12, Willems et al. (2004) found that for τ_kin spans a range between 1 and 210 Myr, and for χ₂ the kinematic age is at least 90 Myr, both of which are compatible with our τkin distributions. Arzoumanian et al. (1999), in turn, found that B1534+12 has likely crossed ɀ = 0 more than once, which is supported by the fact that τ_c likely exceeds τ_kin for χ₂. Arzoumanian et al. (1999) also noted that the BNS B1913+16 can be constrained well kinematically, since it has a low altitude (ɀ) and is moving away from the disc. They estimated τ_kin for χ₁ to be ~5 Myr, which is only slightly below our estimates, and τ_kin for χ₂ to be 60–80 Myr, which exceeds our estimates for χ₂ but aligns better with the characteristic age. The estimates of Willems et al. (2004) are similar, with τ_kin = 2–4 Myr for χ₁ and τ_kin ≳ 55 Myr for χ₂ (and are also in agreement with the analysis of Wex et al. 2000).

Moreover, we considered the merger times (τ_gw) of the binaries. These are of course not age-estimates themselves, since they give the amount of time needed for these binaries to merge in the future. However, the rate at which the orbital separation of the binaries shrink - due to the emission of gravitational waves – strongly depends on the current orbital separation, meaning the binaries spend most of their time close to their initial orbital separation. This means that, on a population level (but not necessarily on an individual level), we expect binaries with longer merger times to be older than systems with younger merger times (for further discussion see Maoz & Nakar 2024). We show the τ_gw values in Fig. 1 and list them in Table 2 (which agree with values given by e.g. Faulkner et al. 2005; Tauris et al. 2017). The merger times of several BNSs exceed a Hubble time and are significantly longer than the merger times of the other binaries. With the exception of J1930–1852, these binaries indeed have the largest values of τ_c.

Fig. 1 Age-related quantities for the Galactic BNSs as discussed in Sect. 2: the characteristic age τc (brown, through Eq. (2) and the kinematic ages τkin assumming GC isotropy (dark green for χ1 and light green for χ2) and kinematic ages assuming LSR isotropy (dark purple for χ1 and light purple for χ2). The boxes extend from the first to third quartile (25% to 75% percentile) of the distributions, while the whiskers show the 5% and 95% percentiles and the different coloured lines show the medians. The figure also contains lines corresponding to τgw (dash-dotted lines), the Hubble time τH = 13.8 Gyr (dotted line, Komatsu et al. 2011), as well as a background line at t = 40 Myr (dashed line, relevant for Sect. 3). The median values of these distributions and their uncertainties are listed in Table 2. For J0737−3039A/B we show τc of J0737−3039B.

Fig. 2 Total distributions of the magnitudes of the Monte Carlo estimated present-day galactocentric velocity vectors of the BNSs (v), shown in normalised histograms with bins of 20 km/s, for GC isotropy (green) and LSR isotropy (purple). The dotted lines show the corresponding median speeds.

3 Speeds

We are interested in the magnitude of kicks the Galactic BNSs experienced at their formation. As a first attempt to determine the kick velocities of the binaries, we investigated their current speeds, as estimated through the method discussed in Sect. 2.2. For each binary, we obtained 10³ positions and velocities, resulting in a distribution of the present-day galactocentric speeds for each binary. In Fig. 2, we show the combined speed distributions of the BNS systems in our sample (for the velocity distributions for some of the individual BNSs, see Gaspari et al. 2024). For both GC isotropy and LSR isotropy, the speed distributions peak around 200–250 km/s. Since this is comparable to the circular velocity of the Solar System, one could argue that this means the binaries are likely to have experienced low kicks. However, Disberg et al. (2024) show that this is not necessarily the case. Here, we briefly re-create and expand their argument for why it is difficult to infer kicks based on current (galactocentric) speeds. In order to do this, we simulated kicked objects in the Galactic potential (Sect. 3.1) and compared the results to the Galactic BNSs (Sect. 3.2).

3.1 Simulation

The argument of Disberg et al. (2024) is based on a Monte Carlo simulation of objects moving through the Galactic potential after receiving a kick. First, they seeded point-masses in an assumed prior distribution which is described by an exponential disc (convolved with a prescription for the spiral arms, using the work of Chrimes et al. 2021). Here, we created a similar simulation, using the prior distribution from Disberg et al. (2024), and compared it with a simulation using a different prior. That is, we compared it with a simulation using a Gaussian annulus as prior distribution for the positions of the objects. Such a distribution was first proposed by Faucher-Giguère & Kaspi (2006), and adopted by Sartore et al. (2010) to fit the pulsar distribution found by Yusifov & Küçük (2004). They adopted the following Gaussian annulus distribution: $$\rho (R)\propto \exp \left(-\frac{{(R-7.04\text{\hspace{0.17em}kpc})}^2}{2{(1.83\text{\hspace{0.17em}kpc})}^2}\right),$$(5)

where R is the galactocentric radius. In Fig. 3, we seeded 10³ points in both the exponential disc prior and the Gaussian annulus prior, and show their positions (together with the positions of the BNSs). The figure shows the significant difference between these priors: in the exponential disc most points are located near the Galactic centre, while in the Gaussian annulus most points are located at R ≈ R_⊙.

Similarly to Fig. 3, we seeded 10³ point-masses in each prior (determining their initial positions), gave them an initial velocity corresponding to the circular velocity at their initial R0 – as dictated by the MW potential of McMillan (2017) – and added a kick velocity that is isotropically oriented (for a more detailed description, see Disberg et al. 2024). We sampled the magnitude of the kick velocity from a Maxwellian distribution, and repeated this simulation for several Maxwellians, choosing values for σ of 10, 20, 50, 100, 200, and 500 km/s (cf. Hobbs et al. 2005, who estimated the σ for isolated pulsars to be 265 km/s). Having established the initial position and initial velocity vector of the 10³ point masses, we computed their trajectories using GALPY v.1.9.0 (Bovy 2015) and the MW potential of McMillan (2017), evaluating the positions and speeds of the point-masses every Myr for 200 Myr. Since we want to compare the speeds of the kicked objects to the observed speeds of the BNSs, we selected the objects that are, at a certain time, in the solar neighbourhood. We adopted the solar neighbourhood definition of Disberg et al. (2024), who – motivated by the cylindrical symmetry of the system – initialised the orbit of the Sun 12 times, each one azimuthally rotated by 2π/12, and evolved these together with the point-masses. Then, at each point in time, they evaluated the speeds of objects that are within 2 kpc of one of the solar obits.

3.2 Deceleration

We repeated the simulation described above for each Maxwellian (σ = 10, 20, 50, 100, 200, and 500 km/s), and obtained the galac- tocentric speeds shown in Fig. 4. In accordance with the results of Disberg et al. (2024), we find that after a certain amount of time the objects in the solar neighbourhood have decreased speeds: they have been decelerated by the Galactic potential. That is, for higher kicks the objects that we see in the solar neighbourhood have more eccentric orbits, while we are more likely to observe them closer to their apocentre where they have a lower speed (corresponding to the asymmetric drift found by Hansen & Phinney 1997). Because of this, the speeds of older objects in the solar neighbourhood are reduced, and have a median value of ∼150–250km/s, independent of the kick distribution. The speed distributions of these older objects do differ in their spread, paradoxically resulting in the lowest speeds being observed for the highest kicks. For young objects, the speeds do depend on the kicks, but these can include objects that become unbound and are therefore not seen in the solar neighbourhood at later times. Moreover, we find no significant difference between the exponential disc and the Gaussian annulus prior, when it comes to these speeds. After all, orbits with identical eccentricity can be formed at different values of R₀ (albeit with a different kick magnitude and orientation).

We do find that in general there tend to be more objects in the solar neighbourhood for the Gaussian annulus prior, because in this prior more objects are formed near or in the solar neighbourhood. However, the relative amount of objects that cross the solar neighbourhood in our simulation may not be representative of the observable BNSs in the Galaxy, because these are only observable during the limited time they are visible as a radio pulsar. Also, after a certain amount of time the BNS merges and is no longer observable as a BNS, and the merger times may – to some degree – depend on the kick magnitude. Nevertheless, while these effects may be important in estimating the formation rates of BNSs from the observed Milky Way population, it has little impact on our work, which seeks to infer the kicks.

Figure 4 shows how the median speeds decrease rapidly for the different kick distributions, obtaining similar values independent of their initial kicks. In Fig. 5, we plot these medians together, for the exponential disc and the Gaussian annulus prior, to show the timescales of this effect. The figure shows that for the exponential disc prior, the medians for all σ obtain similar values sometime between 20 and 30 Myr (as also found by Disberg et al. 2024), and that for the Gaussian annulus prior the medians become similar around 40 Myr. In this work, we conservatively adopted 40 Myr as the timescale above which we expect the objects in the solar neighbourhood to be decelerated and have similar median speeds (∼150–250 km/s) and therefore provide little information about the kick distribution.

In Fig. 1, we show a line at 40 Myr. We argue that for the Galactic BNSs, it is reasonable to estimate them all to be older than ∼40 Myr, since Fig. 1 and Table 2 show that (1) the characteristic ages of the BNSs are well above 40 Myr, meaning that their true ages retain values above this limit even if τ_c is a considerable overestimation (with the only exception being τ_c for J0737–3039B, but this is not a MSP), and (2) the kinematic ages of the binaries indicate that there is only one likely disc crossing below 40 Myr (the exception being B1913+16, which could have two). Because the Galactic BNSs are all likely to be older than 40 Myr, we expect them to have median speeds around ∼150 – 250 km/s, regardless of the kick distribution. This corresponds to the galactocentric speeds we find, as shown in Fig. 2, with median speeds ∼200–250 km/s. If we could estimate v_r more precisely and get more accurate estimates for the BNSs velocities, we could perhaps differentiate between the speed distributions shown in Fig. 4, but currently the uncertainty in the BNS velocity estimates prevents us from doing so. We therefore conclude that, based on the galactocentric speeds of the BNSs alone, we are not able to constrain their kicks.

Moreover, the simulation whose speeds we show in Fig. 4 allowed us to investigate the accuracy of the GC isotropy and LSR isotropy assumptions. We evaluated the simulated trajectories at t = 0, 50, 100, 150 and 200 Myr, and considered the velocity vectors of the objects that are at these times in the solar neighbourhood (with speeds equal to v). These vectors are decomposed in order to obtain a u_t vector, which we combined with Eqs. (3) and (4) to determine the GC-isotropic and LSR- isotropic estimated distributions of v_r . We obtained the speeds estimated through assumption of isotropy (v_iso) by combining the v_r estimates with the ut vectors. The v_iso /v distributions are then a proxy for the accuracy of the isotropy assumptions. In Fig. 6, we show how the v_iso/v distributions evolve over time, for the different Maxwellian kick distributions and for both GC isotropy and LSR isotropy. For σ ≥ 100 km/s the isotropies do not differ significantly in their accuracy, since high velocities are less affected by subtracting the LSR velocity. For low velocities (σ ≤ 50 km/s), however, the velocity vectors of the simulated objects are more or less aligned with uLSR, due to their initial circular velocity, because of which the GC-isotropic estimates underestimate the true speeds while the LSR-isotropic estimates remain remarkably accurate. This is why we are more confident in the LSR-isotropic velocity estimates shown in Fig. 2.

Fig. 3 Priors for the initial positions in our simulation, shown by seeding 103 points and displaying their density in normalised 2D histograms with bins of 0.5 kpc. The left panel shows the exponential disc prior from Disberg et al. (2024), and the right panel shows the Gaussian annulus prior as defined in Eq. (5). The black stars show the positions of the Galactic BNSs, determined as the mean of the Monte Carlo distributions described in Sect. 2.2 (cf. Gaspari et al. 2024). The dotted lines show R = R⊙ ± 2 kpc, which traces the edges of the solar neighbourhood (as defined by Disberg et al. 2024).

Fig. 4 Evolution of the galactocentric speeds (v) of the objects in the solar neighbourhood, for the different Maxwellian kick distributions (σ = 10, 20, 50, 100, 200, and 500 km/s, for each column), shown in 2D histograms with velocity bins of 30 km/s and time bins of 1 Myr. The top row shows the results using the exponential disc prior (cf. Appendix B of Disberg et al. 2024), while the bottom row shows the results using the Gaussian annulus prior. The black lines correspond to the median speed at each time bin. We note that the colour-scale differs for each panel, because the number of objects in the solar neighbourhood differs for each kick distribution.

Fig. 5 Medians of the galactocentric speeds (vmed) as shown in Fig. 4 for kick distributions described by σ = 10 km/s (light dotted), 20 km/s, (dotted), 50 km/s (dark dotted), 100 km/s (light), 200 km/s, and 500 km/s (dark).

Fig. 6 Accuracy of the velocity distributions determined through isotropy assumptions (viso) for the objects in our simulation using the Gaussian annulus prior (of which the speeds are shown in Fig. 4) that are in the solar neighbourhood at a certain point in time. This accuracy is estimated by dividing these distributions by the actual speeds (v) and evaluating viso/v every 50 Myr, for GC isotropy (top row, through Eq. (3)) and LSR isotropy (bottom row, through Eq. (4)), and for Maxwellian kick distributions with different values of σ (columns), shown in normalised 2D histogram with viso/v bins of 0.1 and time bins of 50 Myr, in units of (500 Myr)–1. The black lines show the medians of the distributions.

Fig. 7 Evolution of the eccentricities (ẽ, Eq. (6)) of the objects in the solar neighbourhood, for the different Maxwellian kick distributions (σ = 10, 20, 50, 100, 200, and 500 km/s, for each column), shown in 2D histograms with eccentricity bins of 0.05 and time bins of 1 Myr. The top row shows the results using the exponential disc prior, while the bottom row shows the results using the Gaussian annulus prior. The black dashed lines show t = 40 Myr, approximating the timescale after which the median velocities have been decelerated (as shown in Fig. 5). Similar to Fig. 4, we note that the colour-scale differs for each panel, because the number of objects in the solar neighbourhood differs for each kick distribution.

4 Orbits

Despite the galactocentric speeds of the BNSs harbouring little information about their kicks, we explored whether we can use the fact that the Galactic trajectories of the BNSs appear to be more or less circular (Gaspari et al. 2024) in order to constrain their systemic kicks. In order to do this, we determined the eccentricity of the Galactic orbits we simulated in Fig. 4 and discuss the relationship between this eccentricity and kick distribution (Sect. 4.1). Moreover, we analysed the Monte Carlo estimated BNS trajectories and estimated their eccentricity (Sect. 4.2).

4.1 Eccentricity

In Sect. 3, we analysed the galactocentric velocities of the BNSs. However, when it comes to inferring kicks it may be more useful to look at their peculiar velocities (which are obtained by subtracting the local uLSR, as determined by the MW potential). After all, for small kicks the orbits of the objects remain more or less circular, which means they have a lower peculiar velocity (v_pec) than objects that receive high kicks, and retain these low velocities after ~40 Myr. Nevertheless, for higher kicks v_pec decreases similarly to v (Disberg et al. 2024), meaning it may still be difficult to differentiate between large kicks using the magnitude of v_pec. Thus, for objects older than ~40 Myr the current magnitude of v_pec is only useful for inferring kicks insofar as it can reveal how circular the Galactic trajectory of an object is.

We therefore aim to investigate directly how the kick distribution affects the resulting Galactic orbits and their eccentricity, and whether we can use this to infer kicks. In particular, we are interested in the eccentricity of the orbits of the objects in Fig. 4, which received kicks from vastly different kick distributions. However, these Galactic orbits are not Keplerian, so they do not have an eccentricity as defined in Keplerian dynamics. Nevertheless, we define an analogue of eccentricity for Galactic orbits, by analysing the orbits of the kicked objects in our simulation and defining the minimum galactocentric radius an object obtains in their orbit as R_min , and the maximum value of R as R_max. Then, we computed our eccentricity analogue: $$\tilde{e}=\frac{R_{\max }-R_{\min }}{R_{\max }+R_{\min }},$$(6)

which we refer to as the eccentricity of the Galactic orbit, or just the BNS’ eccentricity (we note that this is not related to the eccentricity within the binary). Similarly to a Keplerian eccentricity: if R_min = R_max the eccentricity equals 0 and the orbit is perfectly circular, and if R_min = 0 or R_max ≫ R_min the eccentricity goes to 1.

If we compute the Galactic orbits of kicked objects long enough, the values we find for R_min and R_max will correspond to the peri- and apocentres of the Galactic orbits, respectively. However, it is possible that within our simulations some objects do not have enough time to cross their peri- and apocentres, meaning R_min and R_max can depend – to some degree – on the computation time of the simulation (t_max). Nevertheless, even if t_max does not allow for the simulated objects to cross their peri- and apocentres, our estimates of ẽ will still be an approximation of the Galactic eccentricity and tell us something about the circularity of their orbits. For example, if an object received a large kick, it might start moving radially outwards in such a way that R_min = R(t = 0) and R_max = R(t = t_max), for a certain t_max. Nevertheless, this will result in high values of ẽ, even if this value increases slightly if we increase t_max.

For the simulations shown in Fig. 4 (where t_max = 200 Myr), we determined ẽ for each of the 10³ objects. In Fig. 7, we show for each simulation the eccentricities of the objects that are at a certain time in the solar neighbourhood. For each object, their value of e does not change, but at each point in time there are different objects in the solar neighbourhood, resulting in different values of e being observed. The figure shows that for Maxwellian kick distributions with increasing σ, we find objects in the solar neighbourhood with more eccentric orbits. That is, larger kicks perturb the initially circular orbits of the objects more, resulting in higher values of ẽ. Moreover, the difference between t≤40 Myr and t > 40 Myr we observe in the galactocentric velocities of these objects (Fig. 4) is less prominent in their eccentricities. For t ≤ 40 Myr, the objects in the solar neighbourhood are slightly more eccentric, possibly due to objects receiving a significant kick because of which they do not return to the solar neighbourhood within t_max. However, the median eccentricities do depend on σ, contrary to the velocities in Fig. 4, meaning ẽ can potentially allow us to distinguish between different kick distributions. Lastly, we note that in Fig. 7, similarly to Fig. 4, we find no significant differences between the results that use exponential disc and the ones that use the Gaussian annulus prior.

Fig. 8 Kinematic histories of the Galactic BNSs (columns), assuming either GC isotropy (green) or LSR isotropy (purple) in the Monte Carlo estimation, integrated backwards in time for 200 Myr. These trajectories are extended to 1 Gyr when used in Sect. 2.2 to determine the kinematic ages (and obtain z values as shown in Appendix A). The dashed circles show R = R⊙ = 8.122 kpc (GRAVITY Collaboration 2018).

4.2 Trajectories

Since eccentricity can potentially provide insight in kicks, we investigated the eccentricities of the Galactic orbits of the BNSs.

Similarly to Sect. 2.2, we traced the BNS trajectories back by estimating their present-day velocity vector and assuming either GC isotropy or LSR isotropy. In Fig. 8, we show these trajectories traced back for 200 Myr (corresponding to t_max in Fig. 7). The figure shows that, despite the assumption of isotropy, most orbits are well-constrained, and some of these orbits (such as the LSR-isotropic orbits of J0737–3039A/B, J1829+2456, J1411+2551, J0453+1559, and J1518+4804) appear to have a remarkably low eccentricity, while others – and in particular the GC-isotropic orbits – seem to be more eccentric (such as the orbits of B1913+16, J0509+3801, B1534+12, and J1930–1852). However, not all orbits shown in Fig. 8 show plausible kinematic histories of the BNSs. In particular, the GC-isotropic orbits of J0737–3039A/B mostly trace back the binary to large galac- tocentric radii, and do not return closer to the Galactic centre within 200 Myr. Since J0737-3039B has a characteristic age of ~50, and it is unlikely that this binary was formed this far away from the Galactic centre, we do not deem these particular trajectories an accurate representation of the actual kinematic history of this binary. This may suggest that, at least in this case, LSR- isotropic orbits represent a more physical scenario. Nevertheless, the other orbits appear to be plausible, allowing us to estimate the eccentricities of these Galactic orbits.

We used the BNS trajectories to determine R_min and R_max, and then calculated ẽ (through Eq. (6)). This resulted in an ẽ-distribution consisting of 10³ eccentricities, for each binary. We determined these distributions for t_max = 200 Myr and for t_max = 1 Gyr, and repeated this for both isotropy assumptions. Figure 9 shows the resulting ẽ-distributions. For the LSR-isotropic orbits, the eccentricity (1) is well-constrained, (2) obtains values where ẽ ≲ 0.5 (with the exception of J0509+3801), and (3) does not change significantly between t_max = 200 Myr and t_max = 1 Gyr. This means that (1) the eccentricity distributions appear suitable for constraining the BNS kicks, (2) these kicks will probably be relatively small (considering the results from Fig. 7), and (3) within 200 Myr most of these orbits cross their peri- and apocentres (making ẽ mostly independent of t_max if t_max ≥ 200 Myr). For the GC-isotropic eccentricities in Fig. 9, the distributions are broader and obtain higher values of ẽ, meaning that if we deduce kicks from these distributions the GC-isotropic kicks will be larger. In fact, for some binaries the LSR-isotropic values of ẽ are below 0.5, while the GC-isotropic counterparts approach 1 (e.g. B1913+16, J1913+1102, J1829+2456, and J1930−1852). Moreover, some of the GC-isotropic distributions shift - to a certain degree - to higher values between t_max = 200 Myr and t_max = 1 Gyr. After all, for low kicks the assumption of GC isotropy underestimates v and therefore overestimates ẽ.

In particular, the GC-isotropic eccentricities of J0737−3039A/B are dependent on t_max. After all, these orbits are moving away from the Galactic centre (as shown in Fig. 8), meaning that increasing t_max will also increase R_max and therefore ẽ as well. The other GC-isotropic orbits show a similar effect, albeit less strong: there is a slight pile-up at ẽ ≳ 0.8 for t_max = 1 Gyr. Nevertheless, the distributions are mostly stable, due to most trajectories crossing their peri- and apocentres within 200 Myr. The eccentricity of B1534+12 is particularly well-constrained, showing a narrow-peak at ẽ ≈ 0.4, for both isotropies and both values of t_max.

Fig. 9 Eccentricities of the Galactic orbits of the BNSs determined by applying Eq. (6) to the Monte Carlo estimated orbits of each binary, assuming either GC isotropy (left panel) or LSR isotropy (right panel), and shown in 2D histograms with ẽ bins of 0.025. For each BNS, we computed ẽ for the trajectories integrated back up to tmax = 200 Myr (as shown in Fig. 8), and also for the trajectories with tmax = 1 Gyr (as used in Sect. 2.2), and compare both ẽ-distributions (tmax = 200 Myr first, and then tmax = 1 Gyr below).

5 Kicks

Since we are interested in estimating the BNS kicks based on the eccentricities of their Galactic orbits, we investigated the relationship between kick and eccentricity (Sect. 5.1). Then, we applied this relationship to the eccentricity distributions we found for the Galactic BNSs in order to kinematically constrain their kicks (Sect. 5.2).

5.1 Kicks versus eccentricity

In Fig. 7, we already found a relationship between kicks and eccentricities: higher kicks cause more eccentric orbits. In order to show the exact relationship between kick magnitude and eccentricity, we expanded this simulation. That is, we repeated our simulation (as described in Sect. 3.1), but instead of a Maxwellian kick distribution we used a delta function to describe the kicks, and repeated the simulation for values of v_kick between 0 km/s and 650 km/s, with steps of 10 km/s (and for both the exponential disc and Gaussian annulus prior). For these simulations we chose t_max = 200 Myr, since we found in Fig. 9 that the ẽ-values of the BNSs do not change significantly beyond 200 Myr. Increasing t_max beyond this value would therefore not affect our results meaningfully. Then, we evaluated the trajectories of the simulated objects every 1 Myr, and determined the eccentricities of the orbits that are in the solar neighbourhood at that point in time. Because young objects show slightly higher values of ẽ (as found in Fig. 7), we left out these objects and summed all eccentricities we find in between 40 Myr and t_max, effectively determining the relationship between v_kick and ẽ.

However, if we only consider the eccentricity of an orbit without taking its direction into account, this relationship will be bimodal. After all, if we observe an object with ẽ = 0 moving along with the Solar System, we estimate that v_kick = 0 km/s: its initial circular orbit has not been disturbed. If, on the other hand, we observe an object with ẽ = 0 orbiting the Galactic centre in the opposite direction compared to the Solar System, this would mean that it received a large kick that reversed its initial circular orbit that coincidentally resulted in a circular orbit in the opposite direction. This shows that in addition to ẽ, we need to take into account whether the Galactic orbit of the objects are in the same or opposite direction of the Solar System. We call these two cases “prograde” and “retrograde” orbits, respectively. For the BNS trajectories in Fig. 8 and the trajectories of the objects in our simulation, we know whether they are prograde or retrograde, since we can compare their initial velocity vector to the solar velocity, and if v_ϕ/v_ϕ,⊙ ≥ 0 the entire orbit is prograde (due to conservation of angular momentum). For each BNS, ≳90% of their LSR-isotropic trajectories are prograde. The GC-isotropic trajectories are >90% prograde, except for the GC-isotropic orbits of B1913+16 (54%), J1913+1102 (61%), J1829+2456 (77%), and J1930−1852 (73%).

In Fig. 10, we show the results of these simulations, separated into prograde and retrograde orbits, for the exponential disc and the Gaussian annulus prior. If v_kick = 0 km/s, all orbits are prograde with ẽ = 0, and the objects retain their initial circular orbit. Then, when we increase v_kick, the Galactic orbits become more eccentric. If v_kick equals the circular velocity (v_LSR) at a solar radius, which is approximately vLSR of the solar neighbourhood, it becomes possible for the kick to be of equal magnitude but in the opposite direction compared to the initial velocity vector. This results in a net velocity of 0 km/s, causing the object to fall to the Galactic centre, obtain R_min = 0 kpc, and therefore gives ẽ = 1. If v_kick > v_LSR (R_⊙), it is possible for the kick to be pointed in the opposite direction of the initial circular velocity and cause retrograde orbits in the solar neighbourhood. At this point, a larger kick means that the retrograde orbits can obtain velocities closer to v_LSR (R_⊙), resulting in a decrease of the retrograde values of ẽ. This continues up until v_kick = 2v_LSR (R_⊙), where the kick can cause perfectly circular, retrograde orbits. Increasing v_kick beyond this point would result again in more eccentric orbits, approaching ${\lim }_{v_{\text{kick\hspace{0.17em}}}\to \infty }\tilde{e}=1$.

We did not find significant differences between the results that used the exponential disc and the ones that used the Gaussian annulus prior. In particular, for ẽ ≲ 0.6, the distributions are remarkably similar. For ẽ ≳ 0.6, there are a few differences between the two rows in Fig. 9. The prograde orbits that use the Gaussian annulus prior and are extremely eccentric tend to be the results of slightly lower kicks than their exponential disc counterparts. This can be explained by the fact that the exponential disc prior seeds the objects closer the Galactic centre and thus deeper in the Galactic potential, meaning that it takes a larger kick for them to obtain orbits that travel to extremely large galactocentric radii and do not cross the solar neighbourhood between 40 Myr and 200 Myr. In other words, we expect that for highly eccentric orbits (i.e. ẽ ≳ 0.8), the results start to depend – to some de gree – on t_max. This can also be found in the retrograde orbits for the Gaussian annulus: v_kick ≳ 500 km/s results in the most eccentric orbits to not be found in the solar neighbourhood (these kicks are also large enough to cause objects to become unbound, depending on the alignment between kick and initial velocity vector). However, for our purposes the relationship between v_kick and ẽ found in Fig. 10 suffices, since the trajectories of the BNSs mostly contain prograde orbits with ẽ < 0.8, except for a few GC-isotropic orbits (which we are less confident in).

Fig. 10 Number of eccentricities found in the solar neighbourhood for each vkick, shown in a 2D histogram with ẽ bins of 0.05 and bins of 10 km/s on a logarithmic colour scale. For each vkick = n ⋅ 10 km/s, we simulate 103 objects using either the exponential disc (top row) or Gaussian annulus (bottom row) prior, and evaluate their orbits in between 40 Myr and tmax = 200 Myr. The resulting eccentricities are divided into prograde (left column) and retrograde (right column) orbits. The dashed lines show 1 and 2 times the circular velocity (vLSR) at a solar radius.

Fig. 11 Posterior kick distributions resulting from our eccentricity-constrained estimates for GC isotropy (green lines) and LSR isotropy (purple lines) and for the exponential disc (dark lines) and Gaussian annulus (light lines) prior. The panels show the normalised results from integrating the posteriors in Appendix B over ẽ, meaning the lines trace the values of a histogram with bins of 10 km/s. Each panel shows the distributions for an individual BNS except the lower-right panel, which combines all distributions.

5.2 Kick velocities

Now that we have estimated the distributions of ẽ for the BNSs (Fig. 9) and have determined the relationship between v_kick and ẽ (Fig. 10), we combined the two in order to estimate the BNS kicks. That is, for the BNSs we took each value of ẽ in their eccentricity distribution, determined whether the corresponding orbit is prograde or retrograde, and matched the value to the corresponding column in Fig. 10. Then, we took out that column, normalised it, and put it in a posterior distribution. We repeated this (1) for each ẽ in the distribution of Fig. 9, (2) for the exponential disc and the Gaussian annulus prior, and (3) for the GC-istropy and LSR isotropy assumption. Effectively, this means we weighted each column in Fig. 10 based on the ẽ distributions we found in Fig. 9 (and combined the prograde and retrograde orbits). In Appendix B we show the posterior distributions for each binary. We integrated these posteriors over ẽ and normalised the result, which gave the posterior probability distributions of v_kick.

In Fig. 11 we show our estimated distributions of v_kick for each individual BNS, as well as the combined kick distribution. Due to the LSR-isotropic orbits being more circular, they result in smaller estimated kicks compared to the GC-isotropic ones. GC isotropy, after all, tends to underestimate low velocities (Fig. 6), resulting in an overestimated eccentricity (Fig. 9), and therefore also overestimates the corresponding kick magnitude. Nevertheless, for some BNSs (i.e. J0509+3801, B1534+12 and J1411+2551) the LSR-istropic and GC-isotropic estimates agree well. Moreover, the exponential disc and Gaussian annulus prior result in similar kick estimates, with the GC-isotropic estimates for J0737−3039A/B, B1913+16, J1913+1102, and J1930−1852 showing the largest differences. This is caused by their highly eccentric orbits (Fig. 9) and the fact that the major differences between the two priors occur at high eccentricities (Fig. 10). Because of this, together with the fact that (1) some GC-isotropic orbits do not appear to be plausible kinematic histories of the BNSs, (2) the highly eccentric orbits in Fig. 10 start to be dependent on t_max, and (3) LSR-isotropic velocity estimates are more accurate for low velocities, we are more confident in the LSR-isotropic kick estimates.

The lower right panel in Fig. 11 contains the total kick distributions for the Galactic BNSs, showing that the LSR-isotropic estimates peak at ~40–50 km/s and the GC-isotropic estimates peak at ~150–200 km/s. In Appendix C we fit a log-normal and a Maxwellian distribution to these estimates, and show that the LSR-isotropic estimates – which we are more confident in than their GC-isotropic counterparts – peak only slightly below the kick distribution found by O’Doherty et al. (2023), who estimate the kicks of neutron stars with low-mass (non-NS) companions using the method of Atri et al. (2019). We find that our LSR-isotropic estimates (i.e. the average between the results from the exponential disc and Gaussian annulus prior) are well-described by a log-normal distribution that peaks at 43 ± 2 km/s, and We note that v_kick ≤ 50 km/s is usually seen as the region of low or no kicks, since it equals the typical escape velocity of globular clusters (e.g. Atri et al. 2019; Mandel & Müller 2020; Willcox et al. 2021; O’Doherty et al. 2023). In our log-normal fit to the LSR-isotropic estimates, 23% of the kick velocities are ≤ 50 km/s.

The method we used to constrain the BNS kicks is to some degree similar to the method Atri et al. (2019) applied to black hole X-ray binaries. They Monte Carlo estimated the 3D velocity vector of these binaries and then traced their trajectories back in time. Then, instead of analysing the eccentricity of the entire orbit, they looked at the peculiar velocities at the instances when the binaries crossed the disc and considered these as potential kick velocities. Using this method, they obtained a potential kick distribution for each binary (O’Doherty et al. 2023, we note that this method is slightly biased towards low kicks, since these produce more disc crossings, which can be corrected for by considering the same amount of disc crossings for each trajectory). As a comparison between our method and the method of Atri et al. (2019), we considered the first disc crossings of the BNS trajectories (χ₁,used to estimate τ_kin in Sect. 2.2 and displayed in Appendix A), and determined their peculiar velocities at these points. We compare these peculiar velocities to our kick estimates in Appendix D, and find that they match remarkably well. Our method does, however, differ from the method of Atri et al. (2019), meaning different assumptions were made in order to estimate the kicks. For example, our method assumes a prior distribution of objects before they are kicked (Fig. 3), whereas their method assumes that averaging peculiar velocities over disc crossings approximates the kick velocity. Also, we note that both methods assume that objects are initially located at z = 0 with no peculiar velocity, but this can be potentially be altered in order to describe objects for which this is not an accurate description.

6 Conclusion

In this work, we determined the characteristic spin-down and kinematic ages of the Galactic BNSs (Sect. 2) and examined whether their galactocentric speeds are representative of their kicks (Sect. 3). As a novel method to constrain the BNS kicks, we determined the eccentricity of their Galactic orbits (Sect. 4) and estimated their kicks through these eccentricities (Sect. 5). Below, we summarise our findings for each section:

• The characteristic ages (through Table 1 and Eq. (2)) and the kinematic ages (through the velocity vectors estimated with the method of Gaspari et al. 2024, which assumes either GC isotropy or LSR isotropy) of the Galactic BNSs indicate that they are likely to have ages ≳40 Myr (Fig. 1), even though it remains difficult to relate these estimates to the true ages of the BNSs;

• The Galactic BNSs have median galactocentric speeds of ~200 – 250 km/s (Fig. 2). However, using the findings of Disberg et al. (2024) and expanding them with a different prior spatial distribution (Fig. 3), we showed that vastly different kick distributions all lead to median galactocentric speeds in the solar neighbourhood of ~ 150–250 km/s (Fig. 4) after ~40 Myr (Fig. 5). Since the BNS ages likely exceed this value, this means we cannot constrain their kicks based on these speeds. Also, we find that GC-istotropic estimates tend to underestimate low velocities, while LSR-isotropic estimates remain accurate (Fig. 6);

• We find that despite the fact that different kick distributions result in similar galactocentric speeds, the eccentricities (Eq. (6)) of their Galactic trajectories do depend on the kick distribution, also for t ≳ 40 Myr (Fig. 7). This indicates that if we investigate the Galactic orbits of the BNSs (Fig. 8) and determine their eccentricity (Fig. 9), we can potentially infer their kick velocities;

• In order to do this, we simulated populations of kicked objects and determined the eccentricities of the Galactic orbits of the objects that are in the solar neighbourhood between 40 Myr and 200 Myr. We repeated this for different kick values, which resulted in a relationship between kick and eccentricity (Fig. 10). Combining this relationship with the BNS eccentricities we found, we constructed kick estimates for each individual BNS, as well as a total kick distribution (Fig. 11). We find that the exponential disc and Gaussian annulus prior both give similar results, whereas the GC isotropy results in higher kick estimates than the LSR isotropy counterparts. We are more confident in the LSR isotropy estimates, and we find that they are well described by a log-normal distribution peaking at ~40–50 km/s.

Moreover, we note that (1) O’Doherty et al. (2023) have found a kick distribution for neutron stars with low-mass companions that peaks at slightly higher values compared to our (LSR-isotropic) distribution (Appendix C) and that (2) applying the method of Atri et al. (2019) to the Galactic BNSs gives results similar to ours, even if we just limit this to the first disc crossing (Appendix D).

Our results are made uncertain by several factors. For instance, in our simulation, we assumed that when the BNSs are kicked, they are on perfectly circular Galactic orbits (i.e. they have no peculiar velocity) at z = 0 kpc. This means that our results are uncertain to the degree that BNSs can have nonzero peculiar velocities pre-kick and are located at |z| > 0 kpc when they are kicked. It is likely that before the final kick at the time of the second supernova, the progenitors have already acquired peculiar velocities due to the previous SN event. Moreover, the BNS kick distribution we have found is a description of eleven observed Galactic BNSs. Because of this, it is not obvious that the distribution we find is representative of all BNS kicks. It is, for example, conceivable that BNS merger times depend on the magnitude of their kick (e.g. Beniamini & Piran 2024), possibly resulting in a bias against high kicks, or are perhaps influenced by their Galactic trajectories (Stegmann et al. 2024). One could also suspect that highly kicked objects are less likely to be observed close to the Solar System. However, we argue that Fig. 10 captures kicks up to at least 500 km/s; the figure shows that objects that received kicks in between our Galactic BNS kick estimates and ≳500 km/s do cross the solar neighbourhood within the time frame of the simulations. Since these kicks do not seem to be present in our BNS sample (though it is a relatively small sample), we argue that our simulation accounts for this bias up to v_kick ≈ 500 km/s. Nevertheless, we cannot decisively exclude the possibility of an undetected BNS population that receives extremely high kicks (>500 km/s). We do note that these uncertainties apply to our method as well as the method of Atri et al. (2019).

Behroozi et al. (2014) find that in order to explain the observed offsets of SGRBs, one needs an SGRB progenitor kick distribution where one in five objects receives a kick with a magnitude greater than 150 km/s. In the log-normal fit to our (LSR-isotropic) results, the region where v_kick > 150 km/s contains approximately 30% of the probability density, which is more or less consistent with the percentage required by Behroozi et al. (2014). We therefore state that while the BNS kick distribution we find may peak at low kicks (~40–50 km/s), it still contains the amount of high kicks required by Behroozi et al. (2014) to match the galaxy-offset distribution found by Fong & Berger (2013).

Future research could expand on several aspects of our analysis, but it would be particularly interesting to expand the simulations behind Fig. 10. That is, it would be interesting to increase (1) the size of the simulations, (2) the value of t_max, and (3) the range of kicks being sampled in order to estimate the relationship between the kick and the eccentricity of the Galactic orbit more robustly. This could give a more accurate estimate for the kick corresponding to a highly eccentric orbit (ẽ ≳ 0.9), after which it could be applied to objects that are thought to receive higher kicks (such as isolated neutron stars, see e.g. Burrows et al. 2024). After all, our estimation does not use the fact that the objects we describe are BNSs, and it can potentially be applied to all kicked objects in order to shed more light on kicks and kick dynamics. Moreover, it would be interesting to investigate how the fact that BNS kicks consist of two separate kicks – due to the two SNe in the binary – affects our findings. These kicks, in turn, could possibly affect the merger times or radio lifetimes of the binaries, which may provide additional constraints on the estimated BNS kicks. This way, our general method for inferring kicks can be tailored for neutron star binaries.

Acknowledgements

We thank Fiorenzo Stoppa for useful discussions regarding this Master’s project, as well as the referee for comments that helped improve this paper. AJL was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 725246), and NG acknowledges studentship support from the Dutch Research Council (NWO) under the project number 680.92.18.02. In this work, we made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen 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).

Appendix A Disc crossings

In Sect. 2, we discuss different quantities that can tell us something about the ages of the Galactic BNSs (i.e. characteristic ages and kinematic ages), and show the estimated distributions in Fig. 1 and Table 2. The kinematic ages (τ_kin) shown in the table are determined through a Monte Carlo estimation of the BNS orbits, assuming either GC isotropy or LSR isotropy (Sect. 2.2). Figure 1 shows the distributions of τ_kin, but in Figs. A.1 and A.2 we also display the trajectories of the BNS orbits in the z dimension, in order to show how the estimated τkin and their 68% intervals correspond to the actual trajectories of the BNSs and their disc crossings. For some binaries, and mostly in the case of LSR isotropy, the orbits are well-constrained and τ_kin is a good estimate of the disc crossings. For others, the trajectories are less constrained by the isotropy assumptions, resulting in a larger uncertainty on τ_kin (on top of the methodological uncertainties in these kinematic ages). In Fig. A.1, the τ_kin ≈ 5 Myr solution for J0747–3039A/B by Willems et al. (2006) can be seen in the GC-isotropic case.

Fig. A.1 Disc crossings of the first six Galactic BNSs, estimated using the method described in Sect. 2.2, for GC isotropy (left column, through Eq. 3) and LSR isotropy (right column, through Eq. 4). The black lines show the overplotted trajectories of the BNSs, where the dark green and dark purple regions correspond to the 68% intervals of χ1, and the light green and light purple regions correspond to the 68% intervals for χ2 (as given in Table 2). The dashed lines of the same colour indicate the medians in these regions. The brown plusses show τc , where the horizontal line extends to the 68% interval (which does not tend to be visible) and the vertical line is located at the median (which are also given in Table 2).

Fig. A.2 Disc crossings of the latter five Galactic BNSs, estimated using the method described in Sect. 2.2, for GC isotropy (left column, through Eq. 3) and LSR isotropy (right column, through Eq. 4). The figure is similar to Fig. A.1. The black lines show the overplotted trajectories of the BNSs, where the dark green and dark purple regions correspond to the 68% intervals of χ1, and the light green and light purple regions correspond to the 68% intervals for χ2 (as given in Table 2). The dashed lines of the same colour indicate the medians in these regions. The brown plusses show τc, where the horizontal line extends to the 68% interval (which does not tend to be visible) and the vertical line is located at the median.

Appendix B Posteriors

In Sect. 5, we describe how we combine our estimates of the BNSs’ eccentricities (Fig. 9) and the relationship we find between v_kick and ẽ (Fig. 10), to determine the posterior kick estimates (Fig. 11). In order to infer v_kick based on ẽ, this method employs Bayes’ theorem: $$P\left(v_{\text{kick\hspace{0.17em}}}\mid \tilde{e}\right)=\frac{P\left(\tilde{e}\mid v_{\text{kick\hspace{0.17em}}}\right)P\left(v_{\text{kick\hspace{0.17em}}}\right)}{P(\tilde{e})}.$$(B.1)

Since we wanted to relate this theorem to the results of our simulations as shown in Fig. 10, we defined each cell in this figure to have a value equal to N(ẽ|v_kick) and a size equals to ẽ⋅Δv_kick while the amount of objects in each column and row equal N (ẽ) and N(v_kick), respectively. Using these quantities, we determine that P(ẽ|v_kick)=N(ẽ|v_kick)/N(v_kick). Then, we define the prior P(v_kick) as proportional to the size of each row (i.e. P(v_kick) ∝ N(v_kick)), effectively accounting for the fact that the fraction of objects that cross the solar neighbourhood differ for each v_kick, independent of ẽ. These two relations allowed us to determine ẽ: $$P(\tilde{e})={\displaystyle \int P}\left(\tilde{e}\mid v_{\text{kick\hspace{0.17em}}}\right)P\left(v_{\text{kick\hspace{0.17em}}}\right)d{v}_{\text{kick\hspace{0.17em}}}\propto {\displaystyle \sum _{v_{\text{kick\hspace{0.17em}}}}N}\left(\tilde{e}\mid v_{\text{kick\hspace{0.17em}}}\right)=N(\tilde{e}).$$(B.2)

Combining these probabilities with Eq. B.1 then gives $$P\left(v_{\text{kick\hspace{0.17em}}}\mid \tilde{e}\right)\propto \frac{N\left(v_{\text{kick\hspace{0.17em}}}\mid \tilde{e}\right)}{N\left(v_{\text{kick\hspace{0.17em}}}\right)}\frac{N\left(v_{\text{kick\hspace{0.17em}}}\right)}{N(\tilde{e})}=\frac{N\left(\tilde{e}\mid v_{\text{kick\hspace{0.17em}}}\right)}{N(\tilde{e})},$$(B.3)

which corresponds to Fig. 10 with each column normalised individually. This relation allowed us to estimate kicks based on an observed eccentricity. However, we did not have a single value of ẽ but a distribution (D_ẽ), which is why we combined Eq. B.3 with the eccentricity distribution P(ẽ|D_ẽ), which are displayed in Fig. 9, in order to determine the desired kick estimate posteriors from Fig. 11: $$P\left(v_{\text{kick\hspace{0.17em}}}\mid D_{\tilde{e}}\right)={\displaystyle \int P}\left(v_{\text{kick\hspace{0.17em}}}\mid \tilde{e}\right)P\left(\tilde{e}\mid D_{\tilde{e}}\right)d\tilde{e},$$(B.4)

for which we show the product P(v_kick|ẽ)P(ẽ|D_ẽ) in Figs. B.1 and B.2. By collecting normalised columns from Fig. 10, corresponding to the distributions in Fig. 9, integrating the posterior over ẽ and normalising the result, we effectively compute the kick posteriors following Eqs. B.3 and B.4.

Fig. B.1 Posterior distributions of vkick versus ẽ for the first six BNSs (rows), following either the GC-isotropic (two left columns) or LSR-isotropic (two right columns) estimations of ẽ, as shown in Fig. 9. For each isotropy assumption, we computed the posterior using the results from the exponential disc (left column) or Gaussian annulus (right column) prior, as shown in Fig. 10. The posteriors are shown in 2D histograms with ẽ bins of 0.05 and vkick bins of 10 km/s, and are normalised so that the (linear) colour scale shows the probability density (ρ) in units of 10−2 s/km.

Fig. B.2 Posterior distributions of vkick versus ẽ (similar to B.1) for the latter five BNSs (rows), following either the GC-isotropic (two left columns) or LSR-isotropic (two right columns) estimations of ẽ, as shown in Fig. 9. For each isotropy assumption, we compute the posterior using the results from the exponential disc (left column) or Gaussian annulus (right column) prior, as shown in Fig. 10. The posteriors are shown in 2D histograms with ẽ bins of 0.05 and vkick bins of 10 km/s, and are normalised so that the (linear) colour scale shows the probability density (ρ) in units of 10−2 s/km.

Appendix C Fits

In Fig. 11 we show our kinematically constrained kick estimates for the individual BNSs, as well as a total distribution for all BNSs. To the total BNS kicks, we fit a log-normal distribution, $$\rho (x\mid \mu ,\sigma )=\frac{1}{x\sigma \sqrt{2\pi }}\exp \left(-\frac{{(\ln (x)-\mu )}^2}{2{\sigma }^2}\right),$$(C.1)

and a Maxwellian, $$\rho (x\mid \sigma )=\sqrt{\frac{2}{\pi }}\frac{x^2}{{\sigma }^3}\exp \left(-\frac{x^2}{2{\sigma }^2}\right).$$(C.2)

We used a non-linear least-squares fit and obtain the optimal parameter values given in Table C.1. In Fig. C.1, we show the fitted distributions, as well as the results from O’Doherty et al. (2023), who estimate the kicks of neutron stars with low-mass companions, and describe this with a beta distribution. The figure shows that the log-normal distribution seems to fit our results significantly better than the Maxwellians.

Fig. C.1 Total kick distribution as determined through ẽ (dots) together with the fitted log-normal (dashed lines) and Maxwellian (dotted lines) distributions for GC isotropy (left panel) and LSR isotropy (right panel). The dots show the average of the results using the exponential disc prior and the ones using the Gaussian annulus prior, tracing the tops of normalised histograms with bins of 10 km/s. The fitted parameters of the distributions are given in Table C.1. We also show the kick distribution that O’Doherty et al. (2023) find for neutron stars with low-mass companions (black dash-dotted lines).

Appendix D Peculiar velocities

As a comparison between our method and the method proposed by Atri et al. (2019), we compare our kick estimates to the peculiar velocities at the first disc crossings of the BNSs. For each BNS, we take the disc crossings that we used in order to estimate τ_kin for χ₁ (in Sect. 2.2), and determine the peculiar velocity at that point as a potential kick velocity. Now, Atri et al. (2019) do not limit themselves to the first disc crossing, but trace back the trajectories of the black hole X-ray binaries in their sample for 10 Gyr and analyse the disc crossings during this entire period. Nevertheless, we find that even if we limit ourselves to the first disc crossing, our kick estimates match the peculiar velocities at χ₁ well (as we show in Fig. D.1). In particular, the total kick distribution for the BNS sample looks nearly indistinguishable if determined through peculiar velocities at χ₁ , compared to our estimates through ẽ. The most notable differences in the kick estimates occur for J0737–3039A/B, J0509+3801, B1534+12, and J1930–1852.

Fig. D.1 Comparison between our kick estimates as constrained through ẽ (solid lines) and peculiar velocities (vpec) at the first disc crossing (dashed lines) for GC isotropy (green lines) and LSR isotropy (purple lines). For both isotropies, we show the average between the results using the exponential disc prior and the ones using the Gaussian annulus prior (as shown in Fig. 11). The lines trace normalised histograms with bins of 10 km/s.

References

All Tables

All Figures

Fig. 1 Age-related quantities for the Galactic BNSs as discussed in Sect. 2: the characteristic age τc (brown, through Eq. (2) and the kinematic ages τkin assumming GC isotropy (dark green for χ1 and light green for χ2) and kinematic ages assuming LSR isotropy (dark purple for χ1 and light purple for χ2). The boxes extend from the first to third quartile (25% to 75% percentile) of the distributions, while the whiskers show the 5% and 95% percentiles and the different coloured lines show the medians. The figure also contains lines corresponding to τgw (dash-dotted lines), the Hubble time τH = 13.8 Gyr (dotted line, Komatsu et al. 2011), as well as a background line at t = 40 Myr (dashed line, relevant for Sect. 3). The median values of these distributions and their uncertainties are listed in Table 2. For J0737−3039A/B we show τc of J0737−3039B.

In the text

	Fig. 2 Total distributions of the magnitudes of the Monte Carlo estimated present-day galactocentric velocity vectors of the BNSs (v), shown in normalised histograms with bins of 20 km/s, for GC isotropy (green) and LSR isotropy (purple). The dotted lines show the corresponding median speeds.
In the text

Fig. 3 Priors for the initial positions in our simulation, shown by seeding 103 points and displaying their density in normalised 2D histograms with bins of 0.5 kpc. The left panel shows the exponential disc prior from Disberg et al. (2024), and the right panel shows the Gaussian annulus prior as defined in Eq. (5). The black stars show the positions of the Galactic BNSs, determined as the mean of the Monte Carlo distributions described in Sect. 2.2 (cf. Gaspari et al. 2024). The dotted lines show R = R⊙ ± 2 kpc, which traces the edges of the solar neighbourhood (as defined by Disberg et al. 2024).

In the text

Fig. 4 Evolution of the galactocentric speeds (v) of the objects in the solar neighbourhood, for the different Maxwellian kick distributions (σ = 10, 20, 50, 100, 200, and 500 km/s, for each column), shown in 2D histograms with velocity bins of 30 km/s and time bins of 1 Myr. The top row shows the results using the exponential disc prior (cf. Appendix B of Disberg et al. 2024), while the bottom row shows the results using the Gaussian annulus prior. The black lines correspond to the median speed at each time bin. We note that the colour-scale differs for each panel, because the number of objects in the solar neighbourhood differs for each kick distribution.

In the text

	Fig. 5 Medians of the galactocentric speeds (vmed) as shown in Fig. 4 for kick distributions described by σ = 10 km/s (light dotted), 20 km/s, (dotted), 50 km/s (dark dotted), 100 km/s (light), 200 km/s, and 500 km/s (dark).
In the text

Fig. 6 Accuracy of the velocity distributions determined through isotropy assumptions (viso) for the objects in our simulation using the Gaussian annulus prior (of which the speeds are shown in Fig. 4) that are in the solar neighbourhood at a certain point in time. This accuracy is estimated by dividing these distributions by the actual speeds (v) and evaluating viso/v every 50 Myr, for GC isotropy (top row, through Eq. (3)) and LSR isotropy (bottom row, through Eq. (4)), and for Maxwellian kick distributions with different values of σ (columns), shown in normalised 2D histogram with viso/v bins of 0.1 and time bins of 50 Myr, in units of (500 Myr)–1. The black lines show the medians of the distributions.

In the text

Fig. 7 Evolution of the eccentricities (ẽ, Eq. (6)) of the objects in the solar neighbourhood, for the different Maxwellian kick distributions (σ = 10, 20, 50, 100, 200, and 500 km/s, for each column), shown in 2D histograms with eccentricity bins of 0.05 and time bins of 1 Myr. The top row shows the results using the exponential disc prior, while the bottom row shows the results using the Gaussian annulus prior. The black dashed lines show t = 40 Myr, approximating the timescale after which the median velocities have been decelerated (as shown in Fig. 5). Similar to Fig. 4, we note that the colour-scale differs for each panel, because the number of objects in the solar neighbourhood differs for each kick distribution.

In the text

Fig. 8 Kinematic histories of the Galactic BNSs (columns), assuming either GC isotropy (green) or LSR isotropy (purple) in the Monte Carlo estimation, integrated backwards in time for 200 Myr. These trajectories are extended to 1 Gyr when used in Sect. 2.2 to determine the kinematic ages (and obtain z values as shown in Appendix A). The dashed circles show R = R⊙ = 8.122 kpc (GRAVITY Collaboration 2018).

In the text

Fig. 9 Eccentricities of the Galactic orbits of the BNSs determined by applying Eq. (6) to the Monte Carlo estimated orbits of each binary, assuming either GC isotropy (left panel) or LSR isotropy (right panel), and shown in 2D histograms with ẽ bins of 0.025. For each BNS, we computed ẽ for the trajectories integrated back up to tmax = 200 Myr (as shown in Fig. 8), and also for the trajectories with tmax = 1 Gyr (as used in Sect. 2.2), and compare both ẽ-distributions (tmax = 200 Myr first, and then tmax = 1 Gyr below).

In the text

Fig. 10 Number of eccentricities found in the solar neighbourhood for each vkick, shown in a 2D histogram with ẽ bins of 0.05 and bins of 10 km/s on a logarithmic colour scale. For each vkick = n ⋅ 10 km/s, we simulate 103 objects using either the exponential disc (top row) or Gaussian annulus (bottom row) prior, and evaluate their orbits in between 40 Myr and tmax = 200 Myr. The resulting eccentricities are divided into prograde (left column) and retrograde (right column) orbits. The dashed lines show 1 and 2 times the circular velocity (vLSR) at a solar radius.

In the text

Fig. 11 Posterior kick distributions resulting from our eccentricity-constrained estimates for GC isotropy (green lines) and LSR isotropy (purple lines) and for the exponential disc (dark lines) and Gaussian annulus (light lines) prior. The panels show the normalised results from integrating the posteriors in Appendix B over ẽ, meaning the lines trace the values of a histogram with bins of 10 km/s. Each panel shows the distributions for an individual BNS except the lower-right panel, which combines all distributions.

In the text

Fig. A.1 Disc crossings of the first six Galactic BNSs, estimated using the method described in Sect. 2.2, for GC isotropy (left column, through Eq. 3) and LSR isotropy (right column, through Eq. 4). The black lines show the overplotted trajectories of the BNSs, where the dark green and dark purple regions correspond to the 68% intervals of χ1, and the light green and light purple regions correspond to the 68% intervals for χ2 (as given in Table 2). The dashed lines of the same colour indicate the medians in these regions. The brown plusses show τc , where the horizontal line extends to the 68% interval (which does not tend to be visible) and the vertical line is located at the median (which are also given in Table 2).

In the text

Fig. A.2 Disc crossings of the latter five Galactic BNSs, estimated using the method described in Sect. 2.2, for GC isotropy (left column, through Eq. 3) and LSR isotropy (right column, through Eq. 4). The figure is similar to Fig. A.1. The black lines show the overplotted trajectories of the BNSs, where the dark green and dark purple regions correspond to the 68% intervals of χ1, and the light green and light purple regions correspond to the 68% intervals for χ2 (as given in Table 2). The dashed lines of the same colour indicate the medians in these regions. The brown plusses show τc, where the horizontal line extends to the 68% interval (which does not tend to be visible) and the vertical line is located at the median.

In the text

Fig. B.1 Posterior distributions of vkick versus ẽ for the first six BNSs (rows), following either the GC-isotropic (two left columns) or LSR-isotropic (two right columns) estimations of ẽ, as shown in Fig. 9. For each isotropy assumption, we computed the posterior using the results from the exponential disc (left column) or Gaussian annulus (right column) prior, as shown in Fig. 10. The posteriors are shown in 2D histograms with ẽ bins of 0.05 and vkick bins of 10 km/s, and are normalised so that the (linear) colour scale shows the probability density (ρ) in units of 10−2 s/km.

In the text

Fig. B.2 Posterior distributions of vkick versus ẽ (similar to B.1) for the latter five BNSs (rows), following either the GC-isotropic (two left columns) or LSR-isotropic (two right columns) estimations of ẽ, as shown in Fig. 9. For each isotropy assumption, we compute the posterior using the results from the exponential disc (left column) or Gaussian annulus (right column) prior, as shown in Fig. 10. The posteriors are shown in 2D histograms with ẽ bins of 0.05 and vkick bins of 10 km/s, and are normalised so that the (linear) colour scale shows the probability density (ρ) in units of 10−2 s/km.

In the text

Fig. C.1 Total kick distribution as determined through ẽ (dots) together with the fitted log-normal (dashed lines) and Maxwellian (dotted lines) distributions for GC isotropy (left panel) and LSR isotropy (right panel). The dots show the average of the results using the exponential disc prior and the ones using the Gaussian annulus prior, tracing the tops of normalised histograms with bins of 10 km/s. The fitted parameters of the distributions are given in Table C.1. We also show the kick distribution that O’Doherty et al. (2023) find for neutron stars with low-mass companions (black dash-dotted lines).

In the text

Fig. D.1 Comparison between our kick estimates as constrained through ẽ (solid lines) and peculiar velocities (vpec) at the first disc crossing (dashed lines) for GC isotropy (green lines) and LSR isotropy (purple lines). For both isotropies, we show the average between the results using the exponential disc prior and the ones using the Gaussian annulus prior (as shown in Fig. 11). The lines trace normalised histograms with bins of 10 km/s.

In the text