© ESO 2020

1. Introduction

Be/X-ray binaries (BeXRBs) are the most numerous subclass of high-mass X-ray binaries in the Galaxy. They are systems containing a neutron star with a mass-losing Be-type main-sequence companion that is surrounded by a circumstellar decretion disk (e.g., Rivinius et al. 2013). These objects are typically revealed by X-ray activity that is fueled by mass accretion. Most of the mass accretion takes place during periastron passages, when the neutron star passes in the vicinity, in some cases, even through, the decretion disk of the Be star (e.g., Ziolkowski 2002; Reig 2011; Casares et al. 2017).

A supernova explosion occurring in a massive binary leads to a disruption of the system in the majority of cases (see, e.g., Brandt & Podsiadlowski 1995; De Donder et al. 1997; Eldridge et al. 2011; Renzo et al. 2019). For the systems that remain bound, BeXRBs provide a valuable but inherently biased laboratory for studying the physics of the supernova explosions that formed their neutron stars. BeXRBs provide a well-defined and simple population: each hosts a neutron star primary with a mass of ∼1.4 M_⊙ (with the notable exception of MWC 656, which hosts a black hole; Casares et al. 2014) and a secondary star from a relatively narrow spectral distribution that peaks at B0 (Reig 2011). This population nevertheless exhibits a wide variety of properties that encode the information about the past supernova event in the system and the massive binary progenitor. Moreover, the short lifetime of the BeXRB phase, typically ∼10 Myr (van den Heuvel et al. 2000), does not allow parameters such as the neutron star masses, orbital periods, and peculiar velocities to change significantly. These parameters are therefore close to their birth values just after the supernova explosion.

The original idea that there might be subpopulations in the BeXRB population was proposed by Pfahl et al. (2002), who reported a subclass of BeXRBs with low eccentricities and low X-ray luminosities. They proposed that this subclass originates from the binaries where the initially more massive star undergoing a supernova explosion has a rapidly rotating core, which results in a neutron star that has received only a small natal kick. Podsiadlowski et al. (2004) and van den Heuvel (2004) proposed that this low-eccentricity subpopulation might be explained if the neutron stars in these systems underwent an electron-capture supernova (ECSN). The ECSNe are the result of the collapse of an oxygen-neon-magnesium core of a lower mass star (possibly with an initial mass as low as 6 M_⊙ if it is in a tight binary, especially at lower metallicities; Podsiadlowski et al. 2004) as it loses pressure support owing to the sudden capture of electrons by neon or magnesium nuclei (Nomoto 1984, 1987), ejecting little mass in the supernova explosion (normally ≲1 M_⊙) and imparting little to no kick to the nascent neutron star (Podsiadlowski et al. 2004; van den Heuvel 2004). The classical high-eccentricity BeXRB population would then be a result of iron-core-collapse supernovae (CCSNe), which occur after a degenerate iron core forms inside a higher mass star (e.g., Cerda-Duran & Elias-Rosa 2018 and the references therein). CCSNe eject more mass and might impart a substantial kick to the newly formed neutron stars as well. The evidence supporting the existence of two distinct explosion mechanisms is not limited to BeXRBs. Observations of double neutron stars, the bimodal velocity distribution of young pulsars, and a high number of neutron stars retained in globular clusters can be attributed to two supernova explosion mechanisms (Beniamini & Piran 2016; Verbunt et al. 2017; Pfahl et al. 2002).

Examining the spin–orbital period diagram of BeXRBs, Knigge et al. (2011; hereafter KCP) also noted that the BeXRB population consists of two subpopulations, a short-period subpopulation with a characteristic orbital period of P_orb ≈ 40 d and spin period P_s ≈ 10 s, and a long-period subpopulation with P_orb ≈ 100 d and P_s ≈ 200 s. The histogram of P_orb and P_s was used to estimate an approximate threshold dividing the subpopulations. This threshold lies at 60 d and 40 s, respectively. Especially in the case of spin periods, the dip in the histogram is considerably wide, therefore these values need to be considered with caution. The two subpopulations are more clearly separated in the spin period than in the orbital period. KCP theorized that these two subpopulations are also a consquence of two different supernovae types that occur in these binaries, where the ECSNe reportedly produce the short-period subpopulation and the CCSNe produce the long-period subpopulation. The observed BeXRB spin periods are not a direct result of the supernova explosion itself, rather, they evolve during the BeXRB phase toward some equilibrium spin period P_eq. This P_eq is dependent on the orbital period of the BeXRBs, which produces the well-known correlation in the spin–orbital period diagram, but it also depends on other parameters that are expected to be different for ECSNe and CCSNe (e.g., resulting neutron star magnetic field and neutron star mass; Waters & van Kerkwijk 1989).

Cheng et al. (2014) proposed an alternative explanation for the two subpopulations in P_s, where the bimodality in P_s can be ascribed to different accretion modes of the neutron stars in BeXRBs. Here, the BeXRB systems that exhibit giant outbursts tend to have shorter spins. During giant outbursts, the neutron star accretes from a thin disk with a relatively long lifetime, which efficiently transfers mass and angular momentum to the neutron star, so that its spin period reaches P_s ∼ 10 s. For the BeXRBs, which undergo predominantly normal outbursts or have no outbursts at all, the accretion torques are smaller. The accretion flows around the neutron stars within these systems are in the form of advection-dominated accretion flows, meaning that the spin-up is infrequent and ineffective. The sources from this subpopulation then exhibit spin periods of about P_s ∼ 100 s. While the supernova mechanism has some effect on the occurrence and type of outbursts and thus spin periods, it is unclear how dominant it is and what other processes are relevant in this case.

These two competing hypotheses can be tested using the BeXRB parameters that do not depend, or depend only negligibly, on the accretion processes, such as the neutron star mass, orbital eccentricity, and peculiar velocity. ECSNe should produce less massive neutron stars than CCSNe (Nomoto 1987; Podsiadlowski et al. 2004). Thus, provided that the BeXRB subpopulations originate from different types of supernovae, we should be able to observe two subpopulations with different neutron star masses. Unfortunately, this cannot be investigated as the long orbital periods and scarcity of eclipsing BeXRBs means that there are only a handful of systems for which the neutron star mass can be determined. ECSNe are also expected to impart smaller kicks to the neutron stars and expel less matter in the explosion, producing systems with lower eccentricities and lower peculiar velocities (Podsiadlowski et al. 2004). It is expected that the observed eccentricities are still close to the original values just after the supernova explosion because the timescales for tidal circularization for the BeXRBs with orbital periods P_orb ≥ 10 d are significantly higher than the secondary lifetime. As a result, tidal effects should have little effect on the orbit of a typical BeXRB, which has an orbital period well above 10 d (van den Heuvel et al. 2000; Reig 2011). KCP noted that there is a trend toward lower eccentricities in the short-period BeXRB subpopulation, in line with their hypothesis, but this trend is not significant. A small number of systems for which the eccentricity values are available makes it difficult to investigate this using this parameter.

Studying the kinematics of BeXRBs has also been problematic. Many of them have only a weak optical counterpart and/or lie at a considerable distance. Thus, their distances and proper motions have been unreliable, often with different astrometric catalogs giving disparate values of the parallaxes and proper motions for the same system (e.g., Ankay et al. 2001; Gvaramadze et al. 2011). This has made the determination of their peculiar velocities difficult and viable only for the close systems. The situation has changed with the advent of the second Gaia data release (GDR2; Gaia Collaboration 2018; Lindegren et al. 2018), which contains the parallaxes and proper motions for the majority of the confirmed and candidate BeXRBs in the Galaxy. This enables us to derive the peculiar velocities of the Galactic population of BeXRBs in a consistent way. However, at the time of writing, the kinematics of BeXRBs situated in the Magellanic Clouds still has to be studied indirectly.

In this work, we investigate whether kinematic subpopulations of BeXRBs exist and what their origin is. We compute the tangential peculiar velocities of Galactic BeXRBs and conduct statistical tests to determine whether they comprise two distinct subpopulations. For the population in the Small Magellanic Cloud (SMC), we use the distances from the closest star cluster that is assumed to be the birthplace of that particular BeXRB as a proxy for the tangential peculiar velocity and conduct equivalent tests.

2. Milky Way Be/X-ray binary population

To construct our sample of Galactic BeXRBs, we selected the systems classified as such in Walter et al. (2015) and the systems listed in the fourth edition of the Catalogue of High Mass X-ray binaries in the Galaxy (Liu et al. 2006). To select a population of the confirmed BeXRBs, we only selected sources with measured P_s. Sources without an optical counterpart or those without a match in GDR2 were also removed from the analysis. This selection yielded 24 systems out of the total of 32 listed in Walter et al. (2015). Another 11 sources satisfying the same criteria were added from Liu et al. (2006). We chose to discard 4U 1901+03 because its optical counterpart is unknown (Reig & Milonaki 2016; Walter et al. 2015), although it is matched to Gaia DR2 4268774695647764352 in Simbad¹. We also removed SAX J0635.2+0533 because of its rotation-powered (rather than accretion-powered) nature (La Palombara & Mereghetti 2017). The relevant properties of the 33 selected sources are summarized in Table 1 without these two systems. Because we relied on the P_s values to divide the subpopulations, we refer to them as the short-spin subpopulation and the long-spin subpopulation.

We emphasize here that we only studied the tangential (transverse) components of the peculiar velocities. To obtain complete information about the kinematics of the source, the radial component of the velocity is necessary as well. However, there are literature radial velocity measurements for only a handful of the sources from the sample, and Gaia does not provide radial velocity measurements for the early-type stars. Moreover, the measured radial velocities of OB stars are, in general, not accurate because the optical lines are formed in atmospheric layers that have outflow velocities of 20–30 km s⁻¹ (e.g., van Oijen 1989). We show below that this value is comparable to the typical tangential peculiar velocities of the BeXRB population. It is also not feasible to correct for this effect because the outflow velocities are variable: nonradial pulsations and wind fluctuations change them. This means that the measured radial velocities of BeXRB do not reflect the true radial motion of the system. We therefore did not consider the radial velocities and assumed that the peculiar velocity distribution of BeXRBs is isotropic.

Many possible methods can be employed to study the kinematics of BeXRBs. First, it is possible to measure the peculiar velocities by accounting for and removing the kinematic components stemming from the Galactic rotation and the peculiar movement of the Sun relative to its local standard of rest (i.e., relative to the expected motion in the Galaxy; e.g., van den Heuvel et al. 2000; Gvaramadze et al. 2011). This method was adopted to compute the peculiar tangential velocities in this paper. More precise and accurate results might be obtained if it were possible to establish the birthplace of the studied system, such as the parent star cluster or association. In this case, the birthplace can be used to anchor a local standard of rest of the system, and the peculiar velocity of the system can then be obtained by subtracting the proper motion of the birthplace from the proper motion of the system (Ankay et al. 2001; Drew et al. 2018; Lennon et al. 2018; Kalari et al. 2019). This normally results in more precise and accurate results than the previous method, as it is not dependent on us knowing the Galactic rotation curve, the distance to the Galactic center, and the peculiar movement of the Sun, and thus is not affected by the uncertainties in these parameters. However, the census of Galactic open clusters and associations beyond 1–2 kpc is incomplete, particularly because of high interstellar extinction in the Galactic plane where the majority of these objects are located. This is also evident considering the high number of newly discovered star clusters using the GDR2 data (e.g., Castro-Ginard et al. 2018; Cantat-Gaudin et al. 2018; Liu & Pang 2019). Because most of the sample BeXRBs lie at much larger distances, using this method was not practical for the vast majority of sources and therefore was not used. Another possible method to obtain an estimate of the peculiar velocity is to estimate the local standard of rest by averaging the proper motions of the stars that are close in projection on the sky to the studied source and approximately at the same distance as the studied source (Kochanek et al. 2019). Similarly, the resulting peculiar velocity estimate can then be obtained by subtracting this mean local proper motion from the proper motion of the studied source.

An alternative way for studying the kinematical properties of the BeXRB sample is indirectly, by studying the locations of BeXRBs and the sites of recent massive star formation, such as the young open clusters and associations. On the premise that the closest young cluster or association to the particular BeXRB system is its birthplace, the separation between the two can then serve as a proxy for the peculiar velocity. This method was used to study the kinematics of BeXRBs in the SMC (Coe 2005), but it is ill-suited for the Milky Way systems because the cluster and association catalogs are incomplete, as mentioned above.

2.1. Peculiar velocities using the Bailer-Jones scale length

Going from the noisy parallax and proper motion measurements to distances and velocities is non-trivial. The most notable problems are the nonlinearity of the transformation and the positivity constraint of the distance. The parallax measurements often exhibit high relative uncertainties and can even be negative, which is often the case when distant objects such as BeXRBs are considered. The naive methods fail and give unphysical results when these measurements are used. However, these measurements are perfectly valid and still hold informational value, therefore it would be a mistake to discard them. The only viable way to handle these measurements is to use a probabilistic analysis (see Luri et al. 2018 and Bailer-Jones et al. 2018 for a more detailed discussion).

Here we used the parallaxes and proper motions from GDR2 (Gaia Collaboration 2018; Lindegren et al. 2018) to compute the tangential peculiar velocities. We followed the approach outlined in Luri et al. (2018), using the workflow from Bailer-Jones (2017). The distances and tangential velocities were jointly estimated from the parallaxes and proper motions by Bayesian inference, with the prior scale lengths for each object adopted from Bailer-Jones et al. (2018). To obtain the peculiar tangential velocities, we adopted the solar Galactocentric distance R₀ = 8.2 kpc, the circular Galactic rotation velocity Θ₀ = 238 km s⁻¹, and the solar peculiar motion (U_⊙, V_⊙, W_⊙) = (10.0, 11.0, 7.0) km s⁻¹ from Bland-Hawthorn & Gerhard (2016).

2.2. Peculiar velocities using an empirically determined scale length

BeXRBs are predominantly discovered in X-rays, which means that it is possible to detect them at significant distances and/or obscured by several tens of magnitudes of extinction. They are then subject to deep follow-up observations to determine and characterize the stellar counterpart. This may result in their selection function being different from the other field stars. However, we do not expect this to affect our results as significantly as in Gandhi et al. (2019), who studied low-mass X-ray binaries hosting black holes. BeXRBs have significantly shorter lifetimes, therefore it is not possible for them to escape far from their parent populations near the Galactic plane, even at runaway speeds. Considering a runaway velocity of 30 km s⁻¹ ∼ 30 pc Myr⁻¹ (informed by the measured tangential peculiar velocities of BeXRBs by van den Heuvel et al. 2000 and theoretical peculiar velocity predictions from Eldridge et al. 2011 and Renzo et al. 2019) and a BeXRB lifetime after the supernova of 10 Myr, we obtain a migration distance estimate of 300 pc. There is also evidence that BeXRB progenitors tend to remain bound within their parent cluster or associations and only acquire high peculiar velocities later on after the supernova explosion of one of its components (Bodaghee et al. 2012). This further limits the possible displacement from the Galactic-plane massive star population. On the other hand, Treuz et al. (2018) found that there might be problems with the distances determined using the scale lengths from Bailer-Jones et al. (2018), where these distances appear to be an underestimation when compared to the distance values obtained using the conventional methods for sources closer than ∼5 kpc. For sources that lie farther away than this, the trend seems to be reversed.

To investigate this, we determined a new scale length for the selected BeXRBs using the distance measurements collected from the literature. These measurements come from a number of sources and thus were collected using various methods, sometimes applied in conjunction. For some sources, more distance estimates exist, in which case we then preferred to use the estimates from the most recent works. These estimates, together with their errors (if available), are listed in Table 1.

These distances were then fit with the exponentially decreasing space density prior probability model,

$$\begin{array}{c}\hfill P(d)=\frac{1}{2L^3}{d}^2{\mathrm{e}}^{-d/L}\mathrm{if}d>0,\end{array}$$(1)

where d is the distance to the source and L is the scale length, as discussed in Bailer-Jones (2015).

Similarly to Gandhi et al. (2019), an unbinned maximum likelihood algorithm was used for the fit. To quantify the uncertainty of L, we generated randomized ensembles of d_lit values by resampling from a suitable distribution for each object. In most cases, for the objects with a published distance estimate d_lit and its uncertainty, a normal distribution was used for resampling, with the assumed mean and standard deviation corresponding to d_lit and its uncertainty, respectively. For four sources with distance limits, random values were drawn from a uniform distribution characterized by the lower and upper distance limit (GS 0834-430, IGR J22534+6243, AX J1700.2-4220, and 2RXP J130159.6-635806). In the case of XTE J1906+090, where only the lower limit on d_lit is known, we also drew from a uniform distribution, where we assumed the upper limit to be d_lit + 5 kpc. One object in our sample, 3U 1022-55, has no uncertainty on d_lit. In this case, we assumed the uncertainty to be 20 % of the published d_lit. We then resampled from a normal distribution as above.

A total of 100 000 ensembles were randomized, resulting in a mean value of the characteristic scale length of L = 1.74 ± 0.06 kpc (this is somewhat higher than the scale height of the thick disk of 0.7–1.2 kpc; Siegel et al. 2002), with the uncertainty quoted here being the standard deviation of the randomized ensambles. The scale length distribution is plotted in Fig. 1. We then adopted this scale length for all the objects in our sample and followed the same procedure as outlined in Sect. 2.1.

Fig. 1. Histogram of the scale lengths L resulting from the resampling.

2.3. Quality of GDR2 data and astrometric fits

While GDR2 provides astrometric measurements of unprecedented quality and quantity, which allows insight into the kinematics of many BeXRBs for the first time, it still contains some sources for which the solutions are ill-behaved. The astrometric parameters of these sources should be considered suspect and be flagged or filtered out from the analysis. Problems with the astrometry may arise for sources that are located in regions with high source densities, for instance, in the Galactic plane. Another potential caveat is that all GDR2 sources are treated as single stars in the astrometric solution, where binaries do not receive any special treatment (Gaia Collaboration 2018; Lindegren et al. 2018). Owing to the binary and early-type nature of BeXRBs, they might be affected by this problem. It is therefore necessary to examine the quality of the astrometric parameters of the BeXRB sample.

Unreliable astrometric solutions can be empirically identified by considering the distributions of the parallax and proper motion errors at the relevant magnitudes and colors and comparing them to the errors for the objects of interest, or by using the recommended astrometric quality indicators that are included or can be computed from the parameters that are part of GDR2. Informed by Gaia Collaboration (2018), Lindegren et al. (2018), and Lindegren (2018), we retained the systems that satisfied

• duplicated_source = False

• astrometric_excess_noise < 1 mas or astrometric_excess_noise_sig < 2

• ruwe < 1.4 or u < 1.2 × max(1, exp(−0.2(G − 19.5))).

The flag duplicated_source=True indicates observational, cross-matching, or processing problems, or stellar multiplicity, probably leading to problems in the astrometric solution. The astrometric_excess_noise (ϵ_i) is an angular measure of the astrometric goodness of fit, indicating the additional scatter that may arise from the movement of the emission centroid that in turn is due to the motion of the components inside a binary. Last, the cuts based on ruwe and u, which stand for renormalized unit weight error and unit weight error, respectively, ensured the removal of ill-behaved astrometric solutions. The ruwe is obtained by dividing u by a normalization factor that is a function of the source magnitude and color and is included in the GDR2 archive. The u values are obtained through

$$\begin{array}{c}\hfill u=\sqrt{{\chi }^2/(N-5)},\end{array}$$(2)

where χ² is astrometric_chi2_al and N is astrometric_n_good_obs_al, which can both be queried in the GDR2 archive. Cuts based on u were applied to the objects with no color information in GDR2, for which the ruwe values were not available.

The relevant quantities and quality flags for the objects in the Galactic BeXRB sample are listed in Table 2. Based on these quality cuts, we discarded five objects from the subsequent analysis. Table 2 shows that the unmodeled orbital motion due to the binary nature of the BeXRBs does not seem to affect the astrometric solutions for the majority of sources in a significant way. The longest P_orb in the sources is about 262 d for SAX J2239.3+6116 (in’t Zand et al. 2000), which has an astrometric solution well below all considered cut limits. This P_orb is still considerably shorter than the 22-month observing time of GDR2. For periods such as this and shorter, any orbital motion will therefore largely average out (Jennings et al. 2018).

The five discarded sources (4U 0728-25, GX 304-1, GS 1843+00, XTE J1906+090, and IGR J19294+1816) are not outstanding in the BeXRB sample in P_orb, distance, or optical brightness. The last three sources, which were discarded due to the increased ϵ_i, lie relatively close to each other in the same region of the sky. They also have lower values of visibility_periods_used, ranging from 10 to 13, while the mean value for the studied BeXRB sample is 15. We opted to list the computed velocities of these sources in the subsequent tables, but we did not consider them in the statistical analysis as their velocities cannot be considered reliable.

2.4. Kinematics of Galactic BeXRBs

We obtained the peculiar velocities for 33 Galactic BeXRBs showing pulsations (see Table 3). In addition to the five sources we discarded because the astrometry was unreliable, we also decided to remove 1H 1249-637 because of its uncertain nature as a γ Cas analog (these are systems that most likely do not host a neutron star, where the X-ray emission is generated by interactions between magnetic fields on the Be star and its decretion disk; e.g., Smith et al. 2016), yielding a final sample of 27 sources. However, this did not affect the results of the following analysis in any significant way. We split the sample according to the dip in P_s adopted from KCP, which is P_s, split = 40 s. The short-spin subpopulation comprises 7 sources, and the long-spin subpopulation is more numerous, with 20 sources. The distribution of the tangential peculiar velocities, derived using the scale lengths from Bailer-Jones et al. (2018; v_pec, BJ) and the empirically determined scale length from Sect. 2.2 (v_pec, iso), with respect to P_s is shown in Fig. 2.

Fig. 2. Peculiar tangential velocities of Galactic BeXRBs with respect to Ps. The red vertical line delineates the approximate divide (Ps, split = 40 s) between the two BeXRB subpopulations according to KCP. Top: tangential peculiar velocities computed using the scale lengths from Bailer-Jones et al. (2018). Bottom: same as above, but with the scale length L resulting from the resampling we used.

Although only a few BeXRBs fall into the short-spin subpopulation, it is apparent that this subpopulation seems on average to be moving with a higher peculiar velocity than the long-spin subpopulation. Using the means of the tangential peculiar velocity posteriors of the individual BeXRBs, we can estimate the characteristics of the two BeXRB subpopulations. The short-spin subpopulation possesses a mean tangential peculiar velocity of approximately 29 ± 11 km s⁻¹, which for the long-spin subpopulation is about 16 ± 8km s⁻¹. The reported errors are the standard deviations of the velocity mean distributions. The results obtained using the empirically determined scale length derived in Sect. 2.2 are similar, but individual velocity measurements exhibit higher uncertainties. In this case, we obtained 28 ± 11 km s⁻¹ and 16 ± 8 km s⁻¹ for the short- and long-spin subpopulation, respectively.

To test the significance of this difference between the populations, we conducted a two-sample Anderson-Darling test (e.g., Scholz & Stephens 1987). The result is that the two populations are indeed distinct. We also used bootstrap testing to estimate the velocity difference between the populations and quantified the effect size through Cohen’s d (Cohen 1988), which is the difference in subpopulation means ${\overline{X}}_1$ and ${\overline{X}}_2$, standardized by dividing by the standard deviation:

$$\begin{array}{c}\hfill d_C=\frac{{\overline{X}}_1-{\overline{X}}_2}{S},\end{array}$$(3)

where S is the pooled standard deviation:

$$\begin{array}{c}\hfill S=\sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}},\end{array}$$(4)

where n_i and s_i are the size and standard deviation of the subpopulation i, respectively. The results of this testing are summarized in Table 4.

Here and throughout this paper, the quoted p-values have the usual statistical meaning: they represent the probability of obtaining a test statistic at least as extreme as the observed one when the null hypothesis is correct. In the case of the peculiar velocities of Galactic BeXRBs, the null hypothesis is that the data are drawn from the same underlying distribution. We adhered to the classical threshold of p <  0.05 for a significant result.

The divide of P_s, split = 40 s adopted here to split the two subpopulations was estimated from the histogram of log P_s values of the BeXRB pulsar sample studied by KCP. The dip in the log P_s histogram in Fig. 1. of KCP is rather wide, ranging from approximately P_s = 20–80 s. To test the robustness of the above results, we repeated the above analysis for the P_s, split of 20, 60, and 80 s using the v_pec, BJ velocities. The recomputed p-values from the Anderson-Darling test, credibility intervals of the velocity differences between the subpopulations, and Cohen’s d for each P_s, split are listed in Table 5. While for P_s, split = 20 s and 80 s the velocity difference is still significant, at P_s, split = 60 s the significance disappeared.

We also tested an alternative approach of estimating the peculiar velocities by determining the local standard of rest through the proper motions of stars near the estimated distance of the studied system. For each BeXRB in the sample, we queried GDR2 for stars within 30 arcmin radius whose parallax was within 1 σ of the system parallax. This ensured that we obtained ≳1000 stars in the vicinity of each BeXRB that satisfied the parallax criterion. Using the GDR2 proper motions of these stars, we established a local standard of rest by computing the mean and dispersion of these values. Normally, these “field” proper motion values would be subtracted from the proper motion of the binary, and with a distance estimate, this would be used to compute a peculiar velocity estimate, such as in Kochanek et al. (2019). However, upon inspection of the field parameters, we found that in general, the dispersions of the proper motion obtained from the stars near the studied systems are too high for any velocity estimates to be meaningful. This was verified when we estimated the peculiar velocity of each BeXRB by drawing 40 000 random samples from the proper motion of its surrounding field, the BeXRB distance distribution, and the BeXRB proper motion distribution (which was assumed to be the normal distribution centered on the GDR2 proper motion value, with the standard deviation being the proper motion error) while also taking the correlations between the proper motions in RA and Dec into account. This yielded peculiar velocity estimates that also supported the hypothesis that there are two kinematic subpopulations of BeXRBs, separated by P_s = 40 s threshold (p ∼ 0.01). However, large errors on the parameters of the field population caused these velocities to be overestimated and to be affected by much larger errors than the velocities obtained using the previous method. As a result, these peculiar velocity estimates were discarded and were not considered in the analysis.

3. Small Magellanic Cloud Be/X-ray binary population

The SMC contains an unexpectedly high number of BeXRBs. Interestingly, all SMC high-mass X-ray binaries, with the exception of SMC X-1, are BeXRB systems. Therefore and because of its relative proximity, the SMC provides a unique laboratory for studying BeXRBs in a homogeneous and consistent manner (Coe & Kirk 2015; Haberl & Sturm 2016).

Still, the SMC distance poses problems for the current astrometric missions. Because of this and its nature as an irregular galaxy, it is not possible to investigate the peculiar velocities of a sufficient number of BeXRBs directly. Using GDR2, Oey et al. (2018) studied the kinematics of early-type SMC runaways, including 14 BeXRBs. However, the P_s for only seven of them is listed in the HMXB catalog of Haberl & Sturm (2016). It is interesting to note that the mean peculiar tangential velocity of these BeXRBs is v_tan,pec ∼ 30 km s⁻¹, which is higher than the observed mean peculiar tangential velocity of Galactic BeXRB (v_tan,pec = 15 ± 6 km s⁻¹; van den Heuvel et al. 2000). The same result is obtained regardless of the method that is used, and it remains the same whether all BeXRBs are considered or only the X-ray pulsars. A possible reason for this higher velocity may be that the metalicity of the SMC is lower than that of the Milky Way (Renzo et al. 2019). Unfortunately, only one BeXRB with P_s <  40 s is included in their analysis. While this system has the lowest peculiar residual velocity of all the sources they studied (not just the BeXRBs), its peculiar velocity as determined locally (using the kinematics of the nearby OB stars) is substantial. Nevertheless, using only one system to characterize a subpopulation is not meaningful.

Therefore, the velocities of a large sample of SMC BeXRBs need to be studied indirectly. In this section we investigate the mutual positions of BeXRBs in the SMC and nearby young star clusters where they might have formed and are now running away after acquiring a high peculiar velocity after the supernova explosion within the progenitor binary. This approach was adopted by Coe (2005), who computed a mean peculiar peculiar velocity arising from the supernova explosion for the BeXRBs in the SMC. Using 17 BeXRB pulsars, their value v_tan,pec ∼ 16 km s⁻¹ is in line with the mean tangential peculiar velocity of Galactic BeXRBs (van den Heuvel et al. 2000). Interestingly, the mean age of the clusters associated with the BeXRB pulsars in Coe (2005) is rather high, 130 ± 140 Myr (log(t/yr) ∼ 8.1). While no firm conclusions can be drawn from a value with such an uncertainty, at face value, this seems much higher than the main-sequence lifetime of a B0 V star (the system secondary), even when we consider that the star would be rejuvenated (i.e., its evolutionary clock will be reset) by the mass transfer within the binary, which would still only yield a maximum lifetime of ∼40 Myr. A possible explanation would be that the cluster ages are only poorly determined, which is common when young star clusters are considered (e.g., Netopil et al. 2015). This would also account for the high uncertainty of the mean cluster age derived in Coe (2005). Another possibility is that the component masses of the SMC BeXRB progenitors may initially be as low as 7–8 M_⊙. It is possible for stars in this mass range to explode as ECSNe, especially if they have lower metallicities, as is the case for the stars in the SMC (Podsiadlowski et al. 2004) During the binary progenitor evolution, the initially more massive star can transfer a substantial amount of mass to the secondary through Roche-lobe overflow, increasing its mass by several M_⊙. Its subsequent evolution will be very similar to the evolution of an isolated star with higher mass (e.g., Podsiadlowski et al. 1992; Pfahl et al. 2002). This would then account for the fact that the spectral distribution of SMC BeXRBs is consistent with that of the Milky Way (Reig 2011).

3.1. Clusters from Rafelski & Zaritsky

Recently, there has been a sharp increase in the number of SMC BeXRBs (Coe & Kirk 2015; Haberl & Sturm 2016). Therefore, it is worthwhile to repeat the analysis done by Coe (2005) on a larger sample size and also look for possible BeXRB subpopulations. We used the catalog of high-mass X-ray binaries in the SMC by Haberl & Sturm (2016), listing 147 BeXRBs, where we selected all pulsating BeXRBs with precisely determined positions for the further analysis. This selection resulted in a list of 56 sources. Similarly to Coe (2005), the SMC clusters comes from the list of Rafelski & Zaritsky (2005; hereafter RZ clusters).

We compared the projected positions of BeXRBs and the RZ clusters. The position of every BeXRB was compared to all RZ clusters and the distance to the closest cluster was obtained. After this, we removed all BeXRBs that lay in the regions that are not covered by the RZ clusters catalog (objects with distances to the closest cluster > 25 arcmin). After this cut, we retained 53 BeXRBs for the analysis.

We caution that this method, while being simple, has significant drawbacks. First, it is in general not possible to establish whether the matched BeXRB/cluster pairs are really equidistant, that is to say, to determine the radial distance offset between them. It is also difficult to determine whether the matched cluster is really the birthplace of the BeXRB, especially when two or more clusters have similar separations from the particular BeXRB. Therefore, automatically picking the closest cluster may not necessarily be correct. Filtering out the old clusters that cannot be the birthplaces of the currently observed BeXRBs alleviates the problem somewhat, but it is apparent that the results obtained using this method need to be interpreted with caution. More precise proper motions from the future Gaia data releases may allow us to confirm the relative system-cluster positions with the peculiar velocities, which would help eliminate some spurious pairings.

Figure 3 shows the distances of the studied BeXRBs from the nearest RZ star cluster. The dashed line indicates the approximate boundary between the two subpopulations, where KCP observed a dip in the spin distribution. The subpopulations appear to be marginally kinematically distinct: the short-spin subpopulation have on average larger distances from the star clusters than the long-spin subpopulation. The mean distance from the nearest cluster of the short-spin subpopulation is approximately 6.6 ± 4.0 arcmin, while for the long-spin subpopulation, it is 4.4 ± 2.3 arcmin. The reported error values are the standard deviations of the distance distributions. The effect size is significant; the Cohen d is 0.75 between the two subpopulations.

Fig. 3. Distances of the SMC BeXRBs from the nearest cluster in the RZ catalog. The conservative estimate of the distance error is 0.7 arcmin, corresponding to twice the median SMC cluster radius listed in the catalog by Glatt et al. (2010). Error bars are omitted for clarity.

To test the significance of this split in the subpopulations, we conducted a two-sample Anderson-Darling test. The test confirmed that the populations are significantly distinct (p = 0.039). The bootstrap testing yielded comparable results: the populations are distinct at a credibility better than 95% (but lower than 99%).

The values we obtained above are sensitive to the inclusion of the two BeXRBs with the shortest P_s, which also exhibit the highest separation from the closest cluster. Excluding them would affect the significance of the subpopulation split, lowering it to p ≈ 0.1. We also used the same cluster catalog as Coe (2005), but it can be expected that the results might change if a different cluster catalog were used. Netopil et al. (2015) studied the inferred parameters of Galactic open clusters in different catalogs and found significant dispersion in the ages of open clusters, with the mean standard deviation of approximately 0.5 dex. Even though the open clusters in the SMC are easier to study in some aspects (the distance to the SMC is known, and the reddening is less variable) than the Galactic open clusters, it is likely that there are significant differences between the SMC cluster catalogs. These differences between the star cluster catalogs that are used are likely to have a significant effect on the results of this analysis. We therefore determined the reliability of this result using the more recent SMC star clusters catalogs.

3.2. Clusters from Nayak et al. (2018)

Nayak et al. (2018) estimated the parameters, including ages, of 174 SMC star clusters. They also collected the parameters of star clusters that were not included in their studied sample, producing a combined catalog of 468 clusters in total, which they used to study the spatio-temporal cluster distribution. We used this catalog to repeat the workflow outlined in the previous section. Because the catalog contains reliable cluster ages, it is possible to study the statistics of clusters/BeXRBs distributions after the age cuts are applied to exclude clusters of a particular age.

To determine the cluster age range in which it is worthwhile to study BeXRB/cluster pairings, we need to estimate the maximum age of the cluster that can be associated with a currently observed BeXRB. This age corresponds to the maximum age of the BeXRB secondary, where it is possible to make a conservative estimate. A star with a mass of 6 M_⊙, which is considered to be a low-mass limit for supernovae if the star is in a binary (Podsiadlowski et al. 2004), has a main-sequence lifetime of ∼110 Myr. If it accretes mass from the primary near the end of its lifetime, it rejuvenates, meaning that it will subsequently evolve like a more massive star, but its evolutionary clock will be reset. As discussed above, we also need to consider the potential uncertainty in estimated cluster ages of approximately 0.5 dex. Thus, in order to include as many viable clusters as possible, it is necessary to include clusters as old as log(t/yr) = 8.6. Clusters older than this age cut contaminate the analysis.

We therefore studied cluster age ranges from log(t/yr) = 8 to 9 because the number of clusters younger than this is lower than or comparable to the number of the studied BeXRBs and the clusters older than log(t/yr) > 8.6 are too evolved for a BeXRBs to be associated with them, even after possible age uncertainties are considered. However, in order to investigate the effect of the contamination by the older clusters, the considered age range was expanded somewhat. We conducted a series of statistical tests in which we applied an age cut every time, starting from log(t/yr) = 8 in 0.1 dex increments and omitted clusters with ages older than this value. As in the previous section, for each BeXRB we searched for the closest cluster within 25 arcmin, collected the distances between them, and compared the BeXRB–closest cluster distance distributions between the spin subpopulations using the P_s, split = 40 s as listed in KCP. The number of retained clusters for the analysis, the significance of the difference between the two BeXRB populations resulting from the two-sample Anderson-Darling test, and the effect size estimated using Cohen’s d are shown in Fig. 4.

Fig. 4. Top: number of retained star clusters when clusters older than a particular age limit from the catalogs of Nayak et al. (2018) and Piatti (2018) are excluded. Middle: significance of the BeXRB–closest cluster distance bimodality as obtained from a two-sample Anderson-Darling test with respect to the applied age cutoff. The horizontal dashed line indicates the classical significance threshold of p = 0.05. The p-values that exceed 0.25 need to be considered as approximate because they are computed by extrapolation. This threshold is indicated by the horizontal dotted black line. Bottom: effect size of the bimodality quantified using Cohen’s d as a function of the applied age cutoff.

3.3. Clusters from Piatti

Piatti (2018) also compiled a large catalog of SMC star clusters, 411 in total, that includes cluster ages. We repeated the same analysis using this catalog. The results are also included in Fig. 4 to facilitate comparison with the clusters from Nayak et al. (2018).

The catalogs of Nayak et al. (2018) and Piatti (2018) are not completely independent. This is shown in Fig. 4, where they follow similar trends, which is most clearly visible around log(t/yr) = 8.3, where the statistics for both catalogs degrade strongly. For both catalogs, the subpopulations appear to be the most distinct overall when only clusters younger than log(t/yr) = 8.2–8.6 are considered, which is reasonable considering the BeXRB lifetimes and likely errors associated with the cluster ages. When older star clusters are included, the statistics degrade, as is expected when clusters are introduced that are unlikely birthplaces of BeXRBs; they contaminate the sample.

The evidence for the BeXRB subpopulation difference is stronger in Piatti (2018) catalog, where the p-value is near the 0.05 level throughout log(t/yr) = 8.2–8.6, and below it at log(t/yr) = 8.2, and 8.4 to 8.5. The subpopulation split is less clear using the sample from Nayak et al. (2018), but even then, for both of the catalogs and the entire age range considered, we always observed the trend that the short-spin subpopulation has higher separations at face value, even though the difference is not always statistically significant.

4. Discussion

The hypothesis proposed by KCP predicts that the long-spin subpopulation of BeXRBs possesses higher peculiar velocities owing to a higher mass-loss and a stronger kick that is exhibited by a CCSN. On the other hand, if the spin bimodality is caused by different accretion modes, as proposed by Cheng et al. (2014), we are less likely to observe any difference in the kinematics of the BeXRB subpopulations. The observed trends, although marginal, challenge both hypotheses.

Overall, we found that the mean peculiar tangential velocity for the Galactic BeXRBs is v_tan,pec = 19 ± 11 km s⁻¹, which is consistent with but somewhat higher than the previous estimates of Chevalier & Ilovaisky (1998) and van den Heuvel et al. (2000). We observed that the short-spin subpopulation of the Galactic BeXRBs exhibits higher tangential peculiar velocities than the long-spin subpopulation. This also holds if the divide between the subpopulations P_s, split is shifted within 20–80 s, but we caution that if P_s, split = 60 s is adopted, the difference between the populations is not statistically significant. Using the distances of SMC BeXRBs from the nearest star cluster, we observed marginal evidence of the same effect in the SMC, but the size and reliability of the effect depends on the star cluster catalog that is used for the analysis. While some problems affect the method we used to investigate the kinematics of SMC BeXRBs, which makes the result less reliable, it appears possible to conclude that there is some evidence that the BeXRB subpopulations are kinematically distinct when we also consider the result we obtained for the Galactic BeXRBs,.

4.1. BeXRBs as runaway systems

We considered the origins of the peculiar velocity of a BeXRB. Two mechanisms induce this additional velocity in BeXRBs: the ejection of a high amount of mass by the supernova explosion, and the asymmetries in the explosion that cause the natal kick for the neutron star. The latter of the two is a consequence of the supernova explosion itself and is thought to be highly stochastic and impossible to describe by a simple analytical prescription. However, scaling relations exist that link the neutron star kick velocity to the mass of the neutrino-heated ejecta, explosion energy, and asymmetry, and the progenitor structure (Janka 2017; Bray & Eldridge 2016). The first of the two mechanisms is a consequence of a rapid mass loss from a binary. During the explosion, the primary loses a considerable amount of mass in a very short time, as compared to the orbital period of the binary. Even if the explosion is spherically symmetric with respect to the remnant neutron star, there is a strong asymmetry with respect to the center of mass of the system. This causes the system to recoil, and it can be described analytically (Zwicky 1957; Blaauw 1961; Cerda-Duran & Elias-Rosa 2018 and references therein).

In a symmetric supernova explosion, the post-supernova binary only remains bound if the system loses less than half of its mass, in other words, when the supernova ejecta M_ej comprises less than a half of the system’s pre-supernova mass, M_ej  <   (1/2)(M₁ + M₂), where M₁ and M₂ are the masses of the primary and secondary components, respectively (Blaauw 1961; Boersma 1961). When we consider the mass loss at the supernova to be instantaneous, the binary system gains a peculiar velocity, given by

$$\begin{array}{c}\hfill v_{\mathrm{sym}}=(\frac{G{M}_2}{a_{\mathrm{pre}-\mathrm{SN}}}{)}^{1/2}(\frac{M_2}{M_1+M_2}{)}^{1/2}(\frac{M_1-M_{\mathrm{co}}}{M_2+M_{\mathrm{co}}}),\end{array}$$(5)

where M_co = M₁ − M_ej is the mass of the compact object (a neutron star in the studied BeXRBs), G is the gravitational constant, and a_{pre − SN} is the pre-supernova semimajor axis of the binary. The last term on the right-hand side of the equation,

$$\begin{array}{c}\hfill \frac{M_1-M_{\mathrm{co}}}{M_2+M_{\mathrm{co}}}=e,\end{array}$$(6)

is the eccentricity of the remnant BeXRB (Iben & Tutukov 1997).

The acquired peculiar velocity (and eccentricity) can be higher when the supernova explosion was asymmetric and the compact object received a kick at birth. Depending on the orientation and magnitude of the kick, the acquired peculiar velocity is

$$\begin{array}{c}\hfill v_{\mathrm{asym}}=(v_{\mathrm{sym}}^2-2v_{\mathrm{sym}}{v}_{\mathrm{k}}cos\psi +v_{\mathrm{k}}^2{)}^{1/2},\end{array}$$(7)

with

$$\begin{array}{c}\hfill v_{\mathrm{k}}=\frac{M_{\mathrm{co}}w}{M_{\mathrm{co}}+M_2},\end{array}$$(8)

where v_sym is given by Eq. (5), w is the kick velocity of the nascent neutron star, and ψ is the angle between the kick velocity and the pre-supernova relative orbital velocity. Equation (7) shows that the highest peculiar velocity is attained if cos ψ = −1, that is, if the direction of the kick at the supernova explosion was opposite to the direction of the relative orbital velocity (Stone 1982; Gvaramadze et al. 2011). In this case, Eq. (7) can be rewritten simply as

$$\begin{array}{c}\hfill v_{\max }=v_{\mathrm{sym}}+v_{\mathrm{k}}=v_{\mathrm{sym}}+\frac{M_{\mathrm{co}}w}{M_{\mathrm{co}}+M_2}·\end{array}$$(9)

Equation (5) shows that the BeXRB attains the highest peculiar velocity if the binary is as tight as possible before the supernova explosion. Higher relative mass-loss from the supernova explosion will also yield higher peculiar velocities. Another way to increase the peculiar velocity is by the natal kick of the neutron star (however, these kicks can also counteract the kick resulting from a symmetric explosion; Kalogera 1998). The magnitude of this natal kick can potentially be very high (as evidenced by the high velocities of isolated radio pulsars; e.g., Hobbs et al. 2005), but given its random orientation to the pre-supernova relative orbital velocity and the canonical mass of the nascent neutron star of about M_co = 1.4 M_⊙, its impact on the attained systemic velocity is reduced because the neutron star has to drag its massive secondary along, provided that the system has remained bound after the supernova. The contribution of the neutron star kick to the systemic velocity abates with the increasing mass of the secondary component as M_co/(M_co + M₂)≪1.

We caution that BeXRBs are fundamentally a biased population and cannot, in general, be used to make inferences about the general properties of neutron stars and supernovae. Renzo et al. (2019) estimated that about ∼86% of the binary systems are disrupted at the moment of the first supernova that occurs within the system, primarily because of the high natal kick that is imparted to the nascent compact object. Eldridge et al. (2011) estimated a similar disruption rate of ∼80%. The majority of the surviving binaries comes from relatively tight pre-supernova orbits, which are less vulnerable to disruption. The majority of these sources is expected to evolve through a phase when they might be detectable as X-ray sources. When Renzo et al. (2019) also included ECSNe in the modeling, allowing for smaller natal kicks imparted to the neutron star, the fraction of bound systems increased significantly from ∼14% to ∼35%.

As outlined above, if a binary survives the supernova, then its velocity relative to the pre-supernova center of mass changes. Here we make a general and approximate brief estimate of the velocity that is expected to be imparted to the system in the case of the CCSN and ECSN. Podsiadlowski et al. (2004) postulated that the main factor determining the type of the supernova that the system will undergo is the initial orbital period (or equivalently, the binary separation) of the system at the start of its evolution (see their Fig. 2). Using the pre-supernova parameters of their prototype binaries and Eq. (5), we derive an approximate attained peculiar velocity of the surviving binary system for the CCSN scenario of ∼80 km s⁻¹ and ≲10 km s⁻¹ for the ECSN scenario. In the CCSN scenario, the attained peculiar velocity of the surviving binary system is predominantly due to the mass loss from the binary after a supernova explosion (the “Blaauw kick”). However, the impact of the SN asymmetry and neutron star kick also significantly supplements the attained peculiar velocity. If the kick velocity w is drawn from a Maxwellian with σ ∼ 265 km s⁻¹ (motivated by the high peculiar velocities of isolated pulsars; e.g., Arzoumanian et al. 2002; Hobbs et al. 2005; Verbunt et al. 2017), this would correspond to a mean kick velocity² of w ∼ 420 km s⁻¹. Even after rescaling w to v_k using Eq. (8), as the impulse of the kick is shared by the entire system (provided that the binary remains bound; van den Heuvel et al. 2000), this velocity is comparable to v_sym. However, a kick of this magnitude can lead to the disruption of the binary (Renzo et al. 2019), therefore the kick magnitude was likely far lower for the surviving BeXRBs, which in turn decreases its likely contribution to the attained peculiar velocity. In the latter case, we considered the effect of the lower natal neutron star kick of w ∼ 50 km s⁻¹ (a mean velocity resulting from a Maxwellian velocity distribution with σ ∼ 30 km s⁻¹; Pfahl et al. 2002; Renzo et al. 2019), which when “diluted” by the total mass of the system, is also comparable to the velocity imparted by the binary mass loss due to the supernova.

To facilitate the comparison with the measured tangential peculiar velocities, it is necessary to project these velocity estimates onto a plane, assuming that the velocity distribution is isotropic. This can be done by multiplying the velocities by a factor of ∼π/4 (e.g., van den Heuvel et al. 2000). This resulted in approximate expected tangential peculiar velocities of ∼60 km s⁻¹ and ∼6 km s⁻¹ for CCSN and ECSN scenarios, respectively. Furthermore, it is also necessary to correct for the frame of reference because the estimated velocities are in the rest frame of the pre-supernova binaries and the observed velocities are measured using a frame that corotates with the Galactic disk. When we assume that the BeXRB progenitors have formed in OB associations and clusters, we need to account for the velocity dispersion inside these structures, which is typically ∼5 km s⁻¹ (e.g. Mel’nik & Dambis 2017). The velocity of the pre-supernova binary inside an OB association and the velocity acquired after the supernova explosion are both randomly oriented, therefore it is necessary to convolve the theoretical velocity distribution by the velocity dispersion. This does not have any significant effect except to smear out the theoretical velocity distribution (see Renzo et al. 2019) and slightly increase the expected velocity of the systems that underwent an ECSN. Moreover, Reid et al. (2014) noted that star-forming regions may be lagging by about ∼5 km s⁻¹ behind the rotation of the Galactic disk. When the frame of reference tied to the Galactic disk is used, the peculiar velocities derived in this way are therefore slightly overestimated compared to their true values. However, this systematic shift affects all studied Galactic BeXRB in the same way and is also quite small. For the purpose of our analysis, it can be neglected.

It needs to be noted that these general and approximate analytic estimates depend heavily on the previous binary evolution. The simulations reported by Renzo et al. (2019) showed that the neutron star natal kick w is the dominant factor affecting the attained peculiar velocity, not the effect of the sudden mass loss from a binary at supernova explosion. The average effective natal kick for compact objects in systems that remained bound after the supernova was w ∼ 66 km s⁻¹ for their simulation, where the kicks were drawn from the Maxwellian with σ ∼ 265 km s⁻¹. This is far lower than the average w that would be expected, even when supernova fallback is accounted for. The systems experiencing higher compact object natal kicks were more likely to be disrupted.

When we consider the points above, if the CCSNe and ECSNe both occur in the BeXRB progenitors, we expect according to our analytic estimates to observe two subpopulations with tangential peculiar velocities of about ∼60 km s⁻¹ and ∼10 km s⁻¹. Brandt & Podsiadlowski (1995) and Eldridge et al. (2011) also derived a similar value for the CCSN subpopulation, but Renzo et al. (2019) and Kochanek et al. (2019) both derived significantly lower values with a median of ∼20 km s⁻¹ (∼16 km s⁻¹ if projected). For systems that underwent an ECSN, both Renzo et al. (2019) and Kochanek et al. (2019) derived very low peculiar velocities.

A reference frame that corotates with the Galactic disk is prone to many uncertainties. It should still be feasible for modern astrometric missions such as Gaia, however, to distinguish between these two subpopulations.

4.2. Relationship between the peculiar velocity, eccentricity, and orbital period

The orbital solutions available for some of the systems from the Galactic BeXRB population enables us to investigate the relationships between the peculiar velocities, eccentricities and P_orb. The top part of Fig. 5 shows the P_orb–eccentricity dependence for the studied BeXRBs, similar to Fig. 6 of Townsend et al. (2011). The plot shows two populations, a low-eccentricity (with e ≲ 0.2) and a high-eccentricity population, which contains most of the remaining systems. Excluding the low-eccentricity subpopulation from the analysis, Townsend et al. (2011) reported a weak correlation between the log(P_orb) and eccentricity. It is not meaningful for us to perform a detailed correlation analysis because our sample is smaller because we lack the SMC systems. The middle part of Fig. 5 shows the P_orb and peculiar velocity relationship. These two parameters appear to be uncorrelated or independent of each other. This contradicts the result of the simulations by Brandt & Podsiadlowski (1995), who predicted an anticorrelation between the peculiar velocity and post-supernova P_orb.

Fig. 5. Relations between the eccentricities e, orbital periods Porb, and peculiar tangential velocities (vpec, BJ) for the sample of Galactic BeXRBs. Top: Porb vs. e plot for the Galactic BeXRB sample. The region that contains the low-eccentricity class of BeXRBs (Pfahl et al. 2002; the region adopted from Townsend et al. 2011) is marked by dashed lines. These systems are marked in red in the current and subsequent graphs for clarity. GX 304-01 and XTE J1906+090 are omitted from the subsequent figure panels because of their unreliable peculiar velocity that is a result of unreliable astrometry. Middle: Porb vs. tangential peculiar velocity. Bottom: e vs. peculiar tangential velocity.

Recently, Bray & Eldridge (2016) proposed a simple relation for the neutron star natal kick velocity, which is defined as a function of the ratio of the supernova ejecta to the neutron star mass. When the resulting neutron star mass is considered to be about identical for all supernova scenarios, this would mean that the resulting velocity, which is a function of the ejected mass, would in general be proportional to the neutron star natal kick, which affects the post-supernova eccentricity. Moreover, the scaling relations from Janka (2017) suggest the same outcome: a more powerful explosion results in higher peculiar velocities and eccentricities. Therefore, we expect to observe that in general, a high-eccentricity system also exhibits high peculiar velocity. The relation between the orbital eccentricities and peculiar velocities of the systems we studied is shown in the bottom part of Fig. 5. The expected relation is clearly not confirmed observationally, and it may even show the opposite trend to what would be expected if the data were taken at face value. This is puzzling, unless some post-supernova orbital circularization mechanism operates in these objects. This is not expected because as outlined above, the timescale for tidal circularization for main-sequence binaries with P_orb >  16 d is at least a few tens of Myr, which is significantly longer than the BeXRB phase lifetime (typically ∼10 Myr; van den Heuvel et al. 2000). However, for the systems with P_orb ≲ 16 d, the orbits might be partially circularized, especially if we observe the system near the end of its lifetime. Alternatively, it is possible that the orbit circularizes when the neutron star interacts with the decretion disk of the Be star at periastron passage (e.g., Martin et al. 2009). Therefore, it is prudent to regard the observed eccentricities as the lower limits of the post-supernova eccentricities.

However, there is also a high scatter in the data, and possible selection effects also need to be taken into account. Furthermore, we also need to consider (see Eqs. (5) and 6 and also Maccarone et al. 2014) that the eccentricity induced on a system by a spherically symmetric mass-loss event alone is normally rather low for BeXRBs (e <  0.2 and far below this value in most cases), while the effect of the neutron star natal kicks is thought to have more impact on the orbital eccentricity, if it is present (van den Heuvel 2004). When the contribution of the neutron star natal kick w to the attained peculiar velocity is significant, it can be expected for a given w that systems with higher eccentricities exhibit lower peculiar velocities and vice versa, as the energy released by the neutron star natal kick can go either into the kinetic energy of the binary system or into the orbit of the system, which would cause the orbit to become larger and/or more eccentric. Therefore, the orbital energy of the binary and its kinetic energy as a whole would be anticorrelated, meaning that the peculiar velocity, eccentricity, and orbital period are linked. This also connects to the previous part of Fig. 5, where it seems that the log(P_orb) is not correlated with peculiar velocity, therefore the orbital energy seems to be increased by pumping the eccentricity rather than widening the binary system. It is not realistic to expect that all sample systems have received a kick of equal magnitude. Moreover, the quality of the current measurements, the low significance of the correlation, and the number of the systems in the sample are not sufficient to reach any firm conclusions at this point.

While the majority of BeXRBs attained their peculiar velocities at the supernova explosion, as indicated by their young kinematic age compared to the OB runaway stars (Huthoff & Kaper 2002; Bodaghee et al. 2012), it is possible that the sample contains some systems that attained some peculiar velocity prior to the supernova explosion, when their progenitors were dynamically ejected from their parent clusters or associations (Poveda et al. 1967). This two-step ejection scenario (Pflamm-Altenburg & Kroupa 2010) may be responsible for the higher observed peculiar velocities of some systems, possibly accounting for some outliers.

4.3. Origins of the kinematic bimodality

When we consider the points above, there are two possible ways to explain the higher peculiar velocities of the short-spin subpopulation, still within the framework of the different supernovae hypothesis. They might arise from tighter binary progenitors or experience greater relative mass-loss and kicks during the supernova explosion than their counterparts in the long-spin subpopulation. However, this is at odds with the current supernova models, where the ECSNe are not expected to cause high mass-loss and neutron star kicks.

In the previous sections, we assumed that the short-spin subpopulation arises from the ECSNe and the long-spin subpopulation arises from systems that underwent a CCSN, which is the hypothesis put forward by KCP. This disagrees with the scenario proposed by Podsiadlowski et al. (2004), where the short-orbit subpopulation (therefore the short-spin subpopulation because P_orb and P_s are correlated, see, e.g., Fig. 1 in KCP) is instead expected to arise from CCSN. The higher peculiar velocities of the short-spin subpopulation would therefore be explained if this model were adopted instead: the opposite of the hypothesis put forward by KCP. However, we note that the predicted peculiar velocities for the CCSN BeXRB subpopulation, derived from the prototype binaries of Podsiadlowski et al. (2004; which are simplifications of the reality), seem to be notably higher than what is observed.

The implications are interesting when we assume that the bimodality in the spin period is caused by different accretion modes of the neutron stars in BeXRBs. According to Cheng et al. (2014), the bimodal distribution of the spin period is not directly linked to the two supernova channels, but the supernova mechanism can still be one of the factors that can modulate it, through its influence on the orbital period, eccentricity, and misalignment of the Be star disk, which are some of the factors that then probably affect the occurrence of the giant outbursts. Haberl & Sturm (2016) used the SMC BeXRB census to investigate which mechanism is responsible for the bimodal spin distribution. They reported higher long-term variability for the short-spin subpopulation, favoring the accretion mode model of Cheng et al. (2014) as the mechanism behind the bimodality. For a BeXRB to show giant outbursts, the right combination of eccentricity and orbital period is apparently needed (Sidoli & Paizis 2018; Cheng et al. 2014). The relative contribution of other factors is also unclear, such as the activity of the Be star, to the frequency of giant outbursts (e.g., Ziolkowski 2002; Reig 2011). Because the effect of the supernova explosion mechanism is only indirect, it is probably unlikely that any kinematic bimodality is observable at all. If giant outbursts (producing the systems with shorter spin periods) dominate in low-eccentricity systems, which are thought to have an ECSN origin, it is expected that the short-spin subpopulation also has a lower peculiar velocity than the rest of the population, which is not observed. However, this changes when the contribution of the neutron star natal kick to the peculiar velocity is significant because we would expect the systems with low eccentricities to exhibit higher peculiar velocities in general.

5. Summary and conclusions

We investigated the kinematics of BeXRB subpopulations that arise from the neutron star spin bimodality in order to test the supernovae hypothesis proposed by KCP and the accretion mode hypothesis of Cheng et al. (2014). We used the GDR2 astrometry to derive the tangential peculiar velocities for the Galactic BeXRBs. In the SMC we used an indirect approach, where we analyzed the distances of the individual BeXRBs from the nearest star cluster as a proxy for the tangential peculiar velocities.

KCP predicted that the short-spin subpopulation, which arises from ECSNe, has systematically lower peculiar velocities than the long-spin subpopulation, which is produced by CCSNe. However, we found some evidence that the subpopulations are kinematically distinct, but in the opposite way as predicted by KCP. The kinematics of the BeXRB subpopulations is difficult to explain within the supernova hypothesis framework because all possible explanations for the increased peculiar velocity of the short-spin subpopulation are hard to reconcile with the current understanding of the ECSNe and CCSNe. Alternatively, adopting the scenario proposed by Podsiadlowski et al. (2004), which reverses the hypothesis of KCP, appears to resolve the discrepancy between the BeXRB kinematics and the supernova mechanisms.

Although the results derived here are statistically significant, they should be regarded with caution because the sample sizes and the methods we used are limited. Clearly, more reliable determinations of the star cluster parameters in the SMC will provide a more solid base for the type of analysis conducted in this paper. Furthermore, the future Gaia data releases will provide an improved astrometry that can further constrain the kinematic properties of Galactic BeXRBs.

Acknowledgments

I would like to thank the anonymous referee for their helpful comments that led to a significant improvement of the paper. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of the Python language (van Rossum 1995) and the following packages: NumPy (Oliphant 2006); SciPy (Jones et al. 2001); Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018) and matplotlib (Hunter 2007). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

References

All Tables

All Figures

	Fig. 1. Histogram of the scale lengths L resulting from the resampling.
In the text

	Fig. 2. Peculiar tangential velocities of Galactic BeXRBs with respect to Ps. The red vertical line delineates the approximate divide (Ps, split = 40 s) between the two BeXRB subpopulations according to KCP. Top: tangential peculiar velocities computed using the scale lengths from Bailer-Jones et al. (2018). Bottom: same as above, but with the scale length L resulting from the resampling we used.
In the text

	Fig. 3. Distances of the SMC BeXRBs from the nearest cluster in the RZ catalog. The conservative estimate of the distance error is 0.7 arcmin, corresponding to twice the median SMC cluster radius listed in the catalog by Glatt et al. (2010). Error bars are omitted for clarity.
In the text

Fig. 4. Top: number of retained star clusters when clusters older than a particular age limit from the catalogs of Nayak et al. (2018) and Piatti (2018) are excluded. Middle: significance of the BeXRB–closest cluster distance bimodality as obtained from a two-sample Anderson-Darling test with respect to the applied age cutoff. The horizontal dashed line indicates the classical significance threshold of p = 0.05. The p-values that exceed 0.25 need to be considered as approximate because they are computed by extrapolation. This threshold is indicated by the horizontal dotted black line. Bottom: effect size of the bimodality quantified using Cohen’s d as a function of the applied age cutoff.

In the text

Fig. 5. Relations between the eccentricities e, orbital periods Porb, and peculiar tangential velocities (vpec, BJ) for the sample of Galactic BeXRBs. Top: Porb vs. e plot for the Galactic BeXRB sample. The region that contains the low-eccentricity class of BeXRBs (Pfahl et al. 2002; the region adopted from Townsend et al. 2011) is marked by dashed lines. These systems are marked in red in the current and subsequent graphs for clarity. GX 304-01 and XTE J1906+090 are omitted from the subsequent figure panels because of their unreliable peculiar velocity that is a result of unreliable astrometry. Middle: Porb vs. tangential peculiar velocity. Bottom: e vs. peculiar tangential velocity.

In the text