THE VELOCITY DISTRIBUTION OF HYPERVELOCITY STARS

When diffusion is the process that determines how binaries fall into the loss cone, the binary approaching rate is almost independent of the size of the loss cone itself (Lightman & Shapiro 1977).

If we fix a, a mass interval between m_c and m_c + dm_c corresponds to a given velocity interval between v and v+dv. Since the number of ejected stars in these two ranges should be the same, the probability density for a fixed semi-major axis a, f_v(v)|_a, is such that f_v(v)|_adv = f_m dm_c. The relative number of ejected stars with velocity v for a fixed a is then $v \times \left. f_{\rm v}(v) \right|_{\rm a} \propto m_{\rm c}^{-(\alpha -1)}$. Therefore,

$$\begin{equation} v \times \left. f_{\rm v}(v)\right|_{\rm a} \propto \left\lbrace \begin{array}{@{}ll@{}} v^{-2 (\alpha -1)}, & \,\; m_{\rm c} \ll m_*, \\ v^{-3 (\alpha -1)}, &\,\; m_{\rm c} \gg m_*, \end{array}\right. \end{equation} \tag{ 6 }$$

where here and in the following X|_y indicates a quantity X calculated at fix Y and we used Equation (5). In Figure 1, Equation (6) is plotted with thin solid lines for decreasing separation a toward the right. The break occurs at the ejection velocity that corresponds to an equal-mass binary (m_c = m_*), which decreases with a as a^−1/2 (Equation (5)). The probability density associated with this velocity is simply the probability density for a fixed companion mass, $f_{\rm v}(v)|_{m_{\rm c}},$ evaluated at m_c = m_*. Since $\left. f_{\rm v}(v)\right|_{m_{\rm c}} dv = f_{\rm a} da$, we get

$$\begin{equation} v \times \left. f_{\rm v}(v)\right|_{m_{\rm c}} = a \times f_{\rm a} =\, {\rm constant}, \end{equation} \tag{ 7 }$$

which is constant for our standard choice of f_a (Equation (3)).

Zoom In Zoom Out Reset image size

We now specifically consider binaries that produce the fastest ejections. The highest velocities are attained for small binary separations and large masses. However, the sum of the stellar radii—which are proportional to their masses—cannot exceed a. In our analytical treatment, unless otherwise stated, we will ask a simple relation a > R_c + R_*. This implies that for a fixed, small a, the maximum companion mass could be smaller than m_max and therefore the maximum velocity could be less than v(a, m_max).

A more realistic model requires a larger minimum separation. Stars with separations less than a few stellar radii would fill their Roche lobes and quickly merge (the exact consequence of mass transfer depends on the mass ratio, e.g., Vanbeveren et al. 1998). Considering that the less massive member has a much higher density, to avoid such Roche lobe overflow, the separation needs to be larger than ∼R_c + R_c(m_*/m_c)^1/3 if m_c > m_* or ∼R_* + R_*(m_c/m_*)^1/3 if m_c < m_*. For binaries consisting of very different mass members (e.g., m_c ≫ m_*), the simple constraint on the separation quickly becomes R_c while in the realistic model it remains as 2R_c. However, these estimates are still approximations since the stars change their shape, making it easier for mass transfer to occur (Eggleton 1983). We will impose a > 2.5max [R_*, R_c] (e.g., Vanbeveren et al. 1998) when we numerically discuss velocity distributions in Section 6. Since v∝a^−1/2, the highest possible velocity is expected to be reduced by a factor of ∼1.6 or less. As we will see below, the highest velocity is achieved for m_* ∼ m_c, the correction factor for the underestimate of the minimum separation is even smaller ${\sim} \sqrt{2.5/2}\sim 1.1$.

Under the requirement a > R_c + R_*, the maximum allowed mass for the companion, m_{c, max}, is the minimum between m_max and m_x ≈ m_*(a/R_*). The transition between m_x and m_max occurs at a ≈ R_max + R_*. For a < R_max + R_*, binaries with m_c = m_x have their separation nearly equal to the sum of their radii, a ≈ R_* + R_c. We call them contact binaries: i.e., binaries that are on the verge of starting mass transfer. They produce the highest ejection velocities for a given separation,

$$\begin{equation} v(m_{\rm x}) \propto m_{\rm x}^{-1/6}, {\rm for }\; \, m_{\rm x} \gg m_*. \end{equation} \tag{ 8 }$$

The maximum of all v(m_x), v_max, occurs for approximately equal-mass binaries: m_x ≈ m_* and a ≈ R_*, v_max ≈ v_esc(M/m_*)^1/6, where $v_{\rm esc} \approx \sqrt{G m_*/R_*}$ is the escape velocity from the surface of the HVS. The result that more massive companions (m_c > m_*) do not yield higher velocities is a direct consequence of the linear dependence assumed between the radius and the mass of a star. A somewhat shallower relation is obtained combining theory and observations: $R_* \propto m_*^{0.8}$ (Hansen et al. 2004). Even assuming this latter slope, the conclusion that the maximum velocity is for equal-mass binaries would remain unchanged, since it holds for any power-law index steeper than 2/3.

The velocity probability distribution for contact binaries with massive companions, m_c ≫ m_*, can be derived considering the stars within an area of the parameter space dm_c × da, f_v dv = f_a f_m dm_cda, under the condition a/R_* ≃ m_c/m_*,

$$\begin{equation} v \times \left. f_{\rm v} \right|_{m_{{\rm c}}=m_{{\rm x}}} \propto m_{\rm x}^{1-\alpha } \propto \begin{array}{ll}v^{6(\alpha -1)}, & \,\; m_{\rm x} \gg m_*, \end{array} \end{equation} \tag{ 9 }$$

where we used Equation (8), while for m_c = m_max, the relative number of stars is constant (Equation (7))

$$\begin{equation} v \times \left. f_{\rm v}\right|_{m_{\rm c} = m_{{\rm max}}} \propto v^{0}, \end{equation} \tag{ 10 }$$

(see the dotted lines in Figure 1).

From Figure 1, it is evident that the total rate of stars $\dot{N}_{\rm v} \propto v \times f_{\rm v}$ ejected with a given velocity v is dominated by binaries with the least massive companion possible that can give that velocity (thick solid lines),

$$\begin{equation} \dot{N}_{\rm v} \propto \left\lbrace \begin{array}{@{}ll@{}} v^{0}, & \,\; v \ll \tilde{v},\\ v^{-2 (\alpha -1)}, & \,\; v \gg \tilde{v}, \end{array}\right. \end{equation} \tag{ 11 }$$

where $\tilde{v} \sim v_{\rm esc} (m_{\rm min}/m_*)^{1/2} (M/m_*)^{1/6} = v_{\rm max} (m_{\rm min}/m_*)^{1/2}$. In (Equation (11)), the low-velocity tail is due to ejections from binaries with a low-mass companion m_c = m_min on varying separations, while the high-velocity branch is given by dissolving contact binaries with a ≃ R_* (i.e., m_c ≪ m_*). The maximum velocity v_max results from dissolving an equal-mass binary. Note that in Figure 1 there is a sharp cut-off at v_max. This is a consequence of our approximate method: a full integration over all binaries that contribute to a given velocity would result in a more gradual decrease to zero of the velocity distribution.

In the full loss-cone regime, the presence of the BH leaves the star distribution function almost unchanged and nearly isotropic at all J (Lightman & Shapiro 1977). Contrary to the empty loss-cone regime, this implies a disruption probability as a function of a where looser binaries are disrupted more easily. More precisely, the cumulative probability to have D ⩽ 1 is

$$\begin{equation} P_{\rm D \le 1} \propto r_{\rm t}. \end{equation} \tag{ 12 }$$

The velocity probability density for a given binary semi-major axis is then given by f_v(v)|_adv = f_m(m_c) P_{D ⩽ 1} dm_c, where we consider that only a fraction P_{D ⩽ 1} of binaries are disrupted. Therefore, $v \times \left. f_{\rm v}(v)\right|_a \propto m_{\rm c}^{-\alpha +1} \left(M/m_{\rm t}\right)^{1/3}$, and

$$\begin{equation} v \times \left. f_{\rm v}(v) \right|_{\rm a} \propto \left\lbrace \begin{array}{@{}ll@{}} v^{-2 (\alpha -1)}, & m_{\rm c} \ll m_*, \\ v^{-3\alpha +2}, & m_{\rm c} \gg m_*, \end{array}\right. \end{equation} \tag{ 13 }$$

where we used Equation (5). As in the empty loss cone, the break occurs at the velocity for an equal-mass binary, v(m_c = m_*). The distribution is shown in Figure 2 as thin solid lines that correspond to different separations, which decrease rightward. For fixed masses, we instead have $\left. f_{\rm v}(v)\right|_{m_{\rm c}} dv = f_{\rm a}\, P_{\rm D \le 1}\, da$, and the rate at fix mass is

$$\begin{equation} v \times \left. f_{\rm v}(v)\right|_{m_{\rm c}} \propto a \propto v^{-2}. \end{equation} \tag{ 14 }$$

An example is shown by the dashed line in Figure 2 (upper panel) for m_c = m_*. Finally, we derive the relative number of ejected stars associated with the maximum allowed mass, m_{c, max} for a given separation. In the case of contact binaries, we follow the reasoning outlined in the previous section, taking into account the ejection probability. Since f_v dv = f_a f_m P_{D ⩽ 1} dm_cda, it follows

$$\begin{equation} v \times \left. f_{\rm v}\right|_{m_{\rm c}=m_{{\rm x}}} \propto m_{\rm x}^{5/3-\alpha } \propto v^{6\alpha -10}, m_{\rm x} \gg m_*, \end{equation} \tag{ 15 }$$

where we used Equation (8). If the binary is wide enough so that the more massive companion possible is m_c = m_max, then the velocity distribution is given by (Equation (14)),

$$\begin{equation} v \times \left. f_{\rm v}\right|_{m_{\rm c}=m_{{\rm max}}} \propto v^{-2}. \end{equation} \tag{ 16 }$$

Equations (15) and (16) are shown in Figure 2 as dotted lines.

Zoom In Zoom Out Reset image size

Comparing the decay indexes of Equations (13) and (15) with −2 (Equation (14)), one realizes that there are three different regimes in which different types of binaries dominate the overall distribution.

For a typical stellar mass distribution with α > 2 (see Figure 2, upper panel), the main contribution to the total rate of ejected stars with a certain velocity v, $\dot{N}_{\rm v}$, comes from binaries whose companion has the least mass possible (i.e., the most abundant). This is analogous to the case α > 1, in the empty loss-cone regime, and

$$\begin{equation} \dot{N}_{\rm v} \propto \left\lbrace \begin{array}{@{}ll@{}} v^{-2}, & v < \tilde{v},\\ v^{-2 (\alpha -1)}, & v > \tilde{v}, \end{array}\right. \end{equation} \tag{ 17 }$$

where the low-velocity tail corresponds to binaries with companion mass m_c = m_min and the high-velocity branch is due to contact binaries with m_c ≪ m_* (thick solid lines in Figure 2, upper panel).

For a power-law index 4/3 < α ⩽ 2 (see Figure 2, lower panel) the total rate (thick solid lines) is dominated at low velocities by wide binaries (with a = a_max and m_c ⩽ m_*), while at high velocities is dominated by equal-mass binaries (m_c = m_*),

$$\begin{equation} \dot{N}_{\rm v} \propto \left\lbrace \begin{array}{@{}ll@{}} v^{-2(\alpha -1)}, & v < v(a_{\rm max},m_*),\\ v^{-2}, & v > v(a_{\rm max},m_*). \end{array}\right. \end{equation} \tag{ 18 }$$

Since a_max ≫ a_min, it is in fact likely that v⁻² is the only relevant slope for HVSs.

For a very shallow slope of index α < 4/3, most of the power at high velocities comes from contact binaries with m_c ≫ m_* and the distribution follows $\dot{N}_{\rm v} \propto v^{6\alpha -10}$.

The important result here is that the high-velocity slope is very steep, with a slope of −2 or steeper. Comparing our results with those of the previous section, we note that in both regimes $\dot{N}_{\rm v} \propto v^{-2.7}$ for a Salpeter mass function.

In case binaries—especially massive ones—tend to be of equal mass, (m_c = m_*), the velocity distributions are given by the dashed lines in Figures 1 and 2. The relative number of stars in a velocity interval is constant in the empty loss cone (Equation (7)) and $\dot{N}_{\rm v} \propto v^{0}$. In the full loss-cone regime, instead, the relevant distribution is given by Equation (14), and the rate of ejected stars is $\dot{N}_{\rm v} \propto \; v^{-2}$. Interestingly, harder binaries have a smaller J_lc and could be in the full loss-cone regime, while wider binaries may not (Perets & Gualandris 2010).

This would produce a broken power-law distribution from a flat slope to v⁻², where the break is at the separation a where there is the transition between the two regimes.

Above we derived the number of ejected stars produced by binaries of all separations and mass ratios. This allowed us to deduce the total velocity distribution, the one presented in Figures 1 and 2 by a solid thick line. Here we summarize its derivation.

First, we realize that the fastest possible ejection velocity is that obtained from a nearly equal mass, contact binary (m_* ≈ m_c, a ≈ 2.5R_*),

$$\begin{equation} v_{\rm max} \approx 0.56 v_{\rm esc} \left(M \over m_* \right)^{1/6}\approx 3500 \left(m_* \over 3 \,M_{\odot } \right)^{-1/6} \,{\rm km \; s^{-1}}, \end{equation} \tag{ 19 }$$

where the numerical estimates are for v_esc = (2Gm_*/R_*)^1/2 ∼ 600 km s⁻¹ and M = 4 × 10⁶ M_☉. If the binary separation distribution is truncated at some a = kR_* (e.g., due to Roche lobe overflow), then v_max would be smaller by a factor of (k/2.5)^1/2.

Velocities lower than v_max, populating the high-velocity branch, can be obtained either by wider equal-mass binaries (e.g., Figure 2, lower panel), or by binaries with smaller companion masses but constant separation a ≈ 2.5R_* (e.g., Figure 2, upper panel), whichever dominates. In this latter case, the relative frequency of ejected stars is simply given by the higher relative frequency of the lighter companions $\propto m_{\rm c}^{1-\alpha }$. Since the ejection velocity depends on lighter companions as $v \propto m_{\rm c}^{1/2}$, the velocity distribution has a slope of v^{−2(α − 1)}. On the other hand, the relative frequency of low-velocity HVSs that result from wider dissolved binaries depends on the circumstances. We have assumed a × f_a = const., so wider binaries in general are just as common. In the empty loss-cone regime, where the rate of dissolving binaries is independent of the tidal radius, this leads to a flat distribution v × f_v = const. However, in the full loss-cone regime, where the dissolving rate is proportional to the tidal radius, this leads to v × f_v∝r_t∝a∝v⁻².

Given a mass function where⁵ α > 2, the lighter companion slope for a fixed separation a is steeper than the slope for a fixed mass m_c (Figures 1 and 2, upper panel). The high-velocity distribution is therefore given by v × f_v∝v^{−2(α − 1)}, down to the velocity produced by a binary with a ≈ 2.5R_* and a minimal mass companion m_min:

$$\begin{equation} \tilde{v} \approx 0.63 v_{\rm esc} \left(M \over m_{\rm t} \right)^{1/6} \left(m_{\min } \over m_{\rm *} \right)^{1/2} \approx 1580 \, {\rm km \; s^{-1}}. \end{equation} \tag{ 20 }$$

The numerical evaluation is for the same values of Equation (19) and m_min = 0.5 M_☉. Below this velocity, the distribution is either flat, v × f_v = const. (empty loss cone) or v × f_v∝v⁻² (full loss cone).

An implication of our results in both regimes is that, if a_min ≈ 2.5R_*, the high-velocity branch in case of α > 2 depends crucially on the existence and physical properties of contact binaries (i.e., binaries on the verge of starting mass transfer). Observationally, we lack firm constraints. It is not clear if contact binaries follow the same mass distribution as wider binaries, and if the minimal separation is indeed given by the mass transfer limit or is larger.

The initial position and velocity of the binary center of mass are chosen so that the binary approaches a massive BH in a parabolic orbit. The encounter of a circular binary with the BH is characterized by nine parameters: masses of the BH and binary members: M, m_*, and m_c, binary separation a, penetration factor D, binary plane orientation (θ, φ), and binary phase ϕ at the initial distance r₀. As long as a simulation starts at a large enough radius r₀ ≫ r_t, the results are largely independent of it. Since we consider HVSs with m_* = 3 M_☉ or 10 M_☉ ejected from a massive BH with M = 4 × 10⁶ M_☉, we have at most six random variables in our Monte Carlo simulations: m_c, D, θ, φ, ϕ, and a.

For each realization, we numerically integrate the evolution of the binary based on the restricted three-body approximation. More precisely, we use Equations (8)–(10) and (12) in SKR. The numerical code is provided with a fourth-order Runge–Kutta integration scheme. We set our initial conditions at r₀ = 10r_t. There, we assign a penetration factor D, an initial binary phase ϕ, and the binary plane direction, defined by the two angles θ and φ. There are two special inclinations: prograde binaries, whose spins are aligned with the orbital angular momentum of the center of mass, and retrograde binaries, whose spins are instead counter-aligned.

We assign D and the various angles as follows. In the empty loss-cone regime, we assume that D is uniformly distributed between zero and D_max = 2.1, where D_max is the largest value for which disruption can occur for circular binaries with any binary plane inclination.⁷ In fact, for each binary system, there is a different minimum D that a binary can reach without one of the stars being tidally disrupted, D ∼ R/a, where R is the stellar radius of the more massive star. This is generally ≪1; however, it can become of the order of ∼1 for contact binaries. Modifications of the high-speed branch due to stellar disruptions will be important if/when speeds in excess of 1000 km s⁻¹ will be measured. Lacking observational motivation, we will omit this effect in this paper. In the full loss-cone regime, we randomly choose the pericenter distance r_p, instead of D, between 0 and D_maxr_{t, max}, where the maximum tidal radius is r_{t, max} = a_max[M/(m_* + m_min)]^1/3. This can be translated into the dimensionless distance D once we assign a mass for the companion and a binary separation (see below). We also randomly choose the binary plane orientation (θ, φ) and the binary phase ϕ between 0 and π. Since a binary starting with a phase difference π has the same post-encounter energy in absolute value but opposite in sign (SKR), the absolute values of the energies provide the ejection velocities for the whole range 0 < ϕ < 2π.

If the numerical result indicates that the binary survives without disruption or that the minimal star separation during the evolution is smaller than the sum of the star radii (i.e., collision), the realization would be disregarded and another one produced until we have an ejection. The removal of colliding binaries prevents us from including spurious, very high-velocity events in the velocity distribution.

The analytical part of the ejection velocity depends solely on the companion mass and on the binary separation. We consider two mass distributions for m_c: the Salpeter mass function, $f_{\rm m}dm_{\rm c} \propto m_{\rm c}^{-2.35} dm_{\rm c}$ (0.5 M_☉ ⩽ m_c ⩽ 100 M_☉), and an equal-mass distribution: m_c = m_*. The binary separation is assumed to be distributed uniformly in the logarithmic space: f_ada∝da/a with a maximum separation given by the mean distance of stars in the sphere of influence of the BH: $a_{\rm max} \sim n_*^{-1/3} \sim 10^6\, R_{\odot } \sim 2 \times 10^{-2}$ pc where n_* ≈ 10⁵ pc⁻³ is the mean stellar number density. We impose a minimum separation roughly equal to the distance at which the Roche overflow will start a_min = 2.5max [R_*, R_c]. However, we do not directly use this distribution. To more efficiently evaluate the high-velocity tail, we first adopt a steeper distribution, f_ada∝a⁻²da (R_* < a < 10⁶ R_☉). This function produces high-velocity events with higher frequency, since v∝a^−1/2. To reconstruct a velocity distribution that reflects an f_a∝1/a function, we then correct our histogram by counting "a" times each realization in a velocity bin, produced by a certain "a."

The results in this paper are produced with 10⁷ realizations. A resolution study is shown in Figure 3, where the two velocity distributions computed with 10⁷ (black solid line) and 3 × 10⁶ (green solid line) are indistinguishable, especially at high velocities.

Zoom In Zoom Out Reset image size

Our results are shown in Figures 4 and 5 for the empty and full loss-cone regimes, respectively. In these figures, we plot the expected velocity distribution within 100 kpc⁸ from the Galactic center for HVSs with m_⋆ = 3 M_☉ (upper panels) and m_⋆ = 10 M_☉ (lower panels). The distributions with a Salpeter mass function for the companion star are shown with red lines, while blue lines are for the equal-mass binaries. Our computational method enhances the statistics of the high-velocity branches (see Section 6.1) and comparatively gives larger statistical errors at low velocities. Furthermore, the Monte Carlo calculations first provide the velocity distribution in terms of the ejection velocity. When this is converted to the velocity at infinity, a narrow velocity region just above the escape velocity spans a wide region in the logarithmic space. To resolve the low-velocity tail (0 < v < v_G), we would need many velocity bins in the numerical calculations and consequently a much larger number of random realizations. However, since the velocity distribution is expected to be smooth in the narrow velocity region around v_G, we have used a linear fit to the numerical data to construct the low-velocity tail. In contrast, the distribution for velocities v > 185 km s⁻¹ comes directly from our numerical results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Each part of the distribution encodes different information. The mass distribution function affects just the highest velocity branch, beyond $\tilde{v}$. This is because the main contribution to these velocities comes from binaries with same separation (a_min) but different companion masses. In our examples, a Salpeter mass function (red lines) results in a significantly steeper slope than for equal-mass binaries (blue lines). Shallower power-law slopes than the Salpeter one (but greater than 1) would result in intermediate velocity curves, between the blue and the red lines. The slope of the distribution at the right of the peak (between ∼v_G and $\tilde{v}$) is determined by the binary separation distribution, since the main contribution to this branch comes from binaries with the same star masses but different separations. As an example, a flatter distribution (e.g., f_a = const.) than a f_a∝1/a would have more power at larger separation and therefore would comparatively enhances lower velocities (v∝a^1/2), producing a steeper slope in the velocity distribution (we will use this feature when comparing our model with data in the next section). Finally, before the peak, the distribution is determined by the Galactic deceleration.

The positions of the breaks bear important information as well. The peak velocity is a measure of the Galactic escape velocity, and value of $\tilde{v}$ is related to the mass of the lighter possible companion and v_age is due to the HVS lifetime. However, while the peak velocity may be easy to determine, since we necessarily expect more events there, $\tilde{v}$ and v_age may not. For example, in our Figures 4 and 5, the break related to v_age is below 200 km s⁻¹ for a m_⋆ = 3 M_☉. At these velocities, the number HVSs constitute only a tiny fraction of the "background" halo stars, whose Gaussian velocity distribution dominates for |v| < 300 km s⁻¹. For a m_⋆ = 10 M_☉ HVS, the age break occurs at very high velocities (∼3000 km s⁻¹), where, in all cases but one, the slopes are very steep, inhibiting detections. The break associated with $\tilde{v}$ is again beyond the velocity peak for m_* ≳ 10 M_☉ and it is not present at all beyond this mass.

The above features were already understood from our analytical calculations. These give a good qualitative description of the velocity distribution (see Section 5) but our numerical results show differences, especially for the high-velocity branches beyond the peak. In our analytical treatment we assumed $\bar{E} =1$ and constant disruption probability for D < 1. The mean dimensionless energy distribution shows that mean energies (averaged over all angles) $\langle \bar{E} \rangle$ can instead be larger than unity for milder penetrations 0.1 < D < 1. The maximum is $\langle \bar{E} \rangle \approx 1.5$ at D ≈ 0.4. In this range of D, the energy dispersion can also be substantial, becoming of the order of unity at D ≈ 0.1. The consequence is that larger ejection velocities can be obtained and the maximum velocity is larger by a factor of ${\sim} \sqrt{1.5}$ with respect to our analytical expectations (Equation (19)). For instance, for m_* = 3 (m_* = 10), this implies v_max ≈ 4300 km s⁻¹ (v_max ≈ 3500 km s⁻¹). On the other hand, the presence of a non-negligible dispersion in dimensionless energy causes the slopes to be a bit shallower for the highest velocities, and any break beyond $\tilde{v}$ to be smoother than predicted analytically. This is shown, for example, in Figure 3 with a comparison between the solid black (full calculation) and red dashed ($\bar{E}=1$) lines. Our numerical calculations also show that the average disruption probability (averaged over all inclinations and phases) is quite constant for D < 0.1, after which it decreases. It is ≈0.5 at D = 1. This feature makes ejections with $\bar{E} >1$ less frequent, mitigating the effect described above. In other words, if a constant disruption probability was imposed for all Ds, there would be an enhancement of high velocities and the resulting velocity distribution would, in fact, be more different from the analytical expectations at the high-velocity end. Finally, we numerically took into account collisions when the distance between stars becomes smaller than the sum of the two radii. In Figure 3, we show that if not excluded, those realizations yield artificially high ejection velocities (blue dashed line).

Current HVS surveys select HVSs with masses of 3–4 M_☉. Up to now, ∼18 HVSs have been observed within ∼100 kpc, with radial velocities ranging between 300 and 700 km s⁻¹ (e.g., Brown et al. 2012b). Around 300 km s⁻¹, there is an equal number of detected stars that have been classified as bound HVSs: they are thought to share the same physical origin as HVSs, but their velocity does not exceed the local Galactic escape speed. Below 300 km s⁻¹, the population of stars bound to the halo dominates the velocity distribution and discovery of HVSs with only radial velocity measurements may become quite challenging. Beyond 720 km s⁻¹, instead, there are no detections at all, and there is no obvious observational bias that explains it. For our purposes, we will consider only the unbound sample of ∼18 stars, which has a narrow range of observed velocities (∼380 to  ∼ 720 km s⁻¹) in which the lower end tends to be more populated than the upper end (see the histogram in Figure 6).

Zoom In Zoom Out Reset image size

We should note that current observations measure only the radial component of the star velocity; thus, the true three-dimensional velocity can be larger. However, if the observed stars are actually produced in the Galactic center, their trajectories at such large distances—much larger than the Sun distance from the Galactic center— should be almost radial as seen from the Sun. Therefore, we do not expect that a large correction should be applied to our predicted velocities, in order to compare our distributions with observations.

Theoretical arguments suggest that our Galactic center favors a situation in which the loss cone is empty. Massive binaries are observed to have periods that are flat in logarithmic scale (e.g., Kiminki & Kobulnicky 2012). In these conditions, the predicted velocity distribution to compare with data is shown in Figures 6 and 7 as a red solid line (see also Figure 4, upper panel). Encouragingly, the range of observed velocities lies close to our peak velocity (v_G = 800 km s⁻¹), in correspondence with the wide peak of the distribution. However, the data are not including the peak velocity (Figure 6), while the model would predict that velocities beyond the peak should be as common as those below it. In Figures 6 and 7, we also plot our fiducial model for a 3  M_☉ HVS, in the full loss-cone regime of star replenishment (black lines, see also Figure 5, upper panel). Although in this case when the peak ≈460 km s⁻¹ is in the middle of the observed range, the model suffers from the same deficiencies, predicting 40% of HVSs with velocities higher than >750 km s⁻¹.

Zoom In Zoom Out Reset image size

We must add here that observed HVSs may still be decelerating in the Galactic halo at the observed location, while our model calculates the final coasting velocities. Using various Galactic potential models from Kenyon et al. (2008), we found that a few (up to half) of the HVSs can indeed have lower final velocities, but that they remain distributed in the same range 300–700 km s⁻¹. Lower velocities, however, can only strengthen our conclusion that both our fiducial models fail to account for current observations.

Analytically, the shape around the peak can be reproduced by the function $\propto v/(v^2+v_{\rm G}^2)^{-(\beta +2)/2}$, which has a peak at $v_{\rm G}/(\sqrt{\beta +1})$. In our fiducial models, which overpredict high velocities, β = 0 or β = 2. We therefore need steeper velocity distributions in the observed range: i.e., larger β values. We obtain possible β values by comparing our model with the observed cumulative distribution (Figure 7). The blue line is for β = 8 and the green line is for β = 3 plus a sharp cut-off at v ≈ 800 km s⁻¹. They both have a probability ⩽10% for velocities higher than >750 km s⁻¹. The corresponding differential distributions are shown in Figure 6, with the same color scheme. The β = 8 curve has a peak around 270 km s⁻¹, fully embedded—and therefore hidden—in the halo stellar distribution.

What might these distributions mean for our binary/Galaxy physical parameter? As discussed in Section 5, the shape around the peak of the distribution for m_* = 3 m_☉ depends on three physical ingredients. (1) The value of the escape velocity, (2) the binary separation distribution, and (3) the injection mode of binaries into the tidal sphere of influence. Since the escape velocity does not influence the value of β, we may assume that a steeper slope may be given by a combination of the last two points. Both an empty and full loss cone require a binary separation distribution that rises toward larger separations or at most decreases slowly to give more weight to lower velocities than to higher ones. In particular, an empty loss cone requires a binary separation distribution that rises as f_a∝a³ (for β = 8), while the green line in Figure 6 can be reproduced with f_a∝a^1/2 and a minimum separation for the binaries of a_min ≈ 10R_* (equal to 30 R_☉ in this case). A full loss cone already favors wide binaries over tight ones, and therefore requires shallower separation distributions: f_a∝a² and f_a∝a^−1/2. We may speculate on the origin of these distributions. If in the Galactic central bulge binaries had a separation distribution that peaks ≫ 1 AU (or, equivalently, a period ≫30 days), then the shorter period binaries may indeed be described by a rising power-law distribution in separations. To cite a well-known example, binaries in the solar neighborhood have a log-normal distribution in periods that peak at 10⁵ day (Duquennoy & Mayor 1991). A possible cut-off at ∼10R_* may instead be related to the binary formation/evolution process or to the injection mechanism which somehow filters out tight binaries.

The velocity distributions of HVSs in the Galactic halo have been investigated by a few groups through Monte Carlo simulations based on full three-body calculations (e.g., Bromley et al. 2006; Kenyon et al. 2008; Zhang et al. 2010, 2013). Zhang et al. (2010) is especially relevant to our study, and they have shown that the velocity distribution is sensitive to how binaries move into the BH loss cone, and that if binaries approach the BH in parabolic orbits, a significant fraction of the resulted HVSs have velocities larger than the detected values. They also study the interactions between the BH and binary stars bound to it. In this case, the binary would experience multiple encounters with the BH. Due to the cumulative tidal effect, the probability that the binary will be broken up becomes substantial even at a large penetration factor D. The observed distribution would be reproduced if binary stars slowly diffuse onto low angular momentum orbits and most of them are broken up at a large distance with small ejection velocities. This could provide an alternative interesting solution to the problem.

In conclusion, data at this stage seem to suggest the need for more power to wider binaries to account for the paucity of high-velocity stars. The above discussion is, of course, not conclusive or exhaustive, given the quality and quantity of the data. Nevertheless, it shows the great potentiality of HVSs to constrain physical properties of the Galaxy and the Galactic Bulge, once a larger data set of HVSs is collected.