VERY WIDE BINARY STARS AS THE PRIMARY SOURCE OF STELLAR COLLISIONS IN THE GALAXY

We use a simple symplectic algorithm to integrate the orbital evolution of 5 × 10⁵ different very wide binary star systems for 10 Gyr. To increase computing efficiency, these systems are split into batches of 1000 binaries. In these 1000 binary batches, each binary is assumed to have the same primary star mass and whose mass is randomly drawn from a mass function of the present day solar neighborhood (Reid et al. 2002). Once the primary mass is chosen, the secondary mass for each of the 1000 binaries is drawn separately from a uniform mass distribution between 0.1 M_☉ and the primary mass (Raghavan et al. 2010). With the stellar masses of the binaries decided, we next assign orbits to each of the binaries. Semimajor axes are randomly selected from a log-uniform distribution between 10³ and 3 × 10⁴ AU (Öpik 1924; Poveda et al. 2007). The rest of the orbital elements are drawn from an isotropic distribution.

Our binaries are perturbed by the tide of the Milky Way's disk as well as impulses from other passing field stars. For the Milky Way's tide, we employ a model containing both radial and vertical terms (Levison et al. 2001). This tidal model implicitly assumes that the systems being simulated (very wide binaries in our case) are on fixed circular orbits about the Galactic center. In this tidal model, the vertical term is approximately an order of magnitude stronger than the radial term at the Sun's galactocentric distance (∼8 kpc), and the strength of this term is set by the local density of matter in the disk. In the Sun's case, the gas disk and the stellar disk of the Milky Way make roughly equal contributions to the local galactic density of matter (Just & Jahreiß 2010). However, the scale height of the gas disk is approximately three times less than the stellar disk (Just & Jahreiß 2010). As a result, most stars in the Milky Way lie farther away from the Milky Way's midplane than the Sun, and the local galactic matter surrounding them will actually be dominated by the stellar component. Because of this, we assume the density of matter is dominated by stars when calculating the strength of the Galactic tide in our binary simulations. At r = 8 kpc, we assume the total local galactic density is 0.04 M_☉ pc⁻³.

The other external perturbation we must model is the effect of field star passages. Rather than modeling these with direct integrations, we employ the impulse approximation (Rickman 1976). For our binary systems, this is a perfectly acceptable approximation since the orbital periods of even our a = 10³ AU systems are much longer than the encounter timescales of the vast majority of relevant and/or plausible stellar passages. To generate individual stellar passages, we employ a method similar to that used in studies of the Oort Cloud (Rickman et al. 2008). Stellar masses are selected from the observed local mass function (Reid et al. 2002), and encounter velocities are selected from the observed mass-dependent stellar velocity dispersions seen in Hipparcos data (García-Sánchez et al. 2001). The rate of encounters for each stellar mass category is set by the local dispersion and spatial density observed in Hipparcos for that particular mass range.

Rather than generating and applying a unique set of stellar encounters to each one of our 5 × 10⁵ binaries, we instead reuse the same set of stellar encounters 10³ times. New sets of stellar encounters are only generated for each batch of 10³ binaries mentioned in the previous section. Thus, our simulations only use 500 different sets of stellar encounters.

In addition to modeling the orbital evolution of binaries in the Sun's region of the Milky Way disk, we also simulate binaries at four other galactocentric distances: two that are closer to the Galactic center and two farther away. Each galactic position is separated by one disk scale length (∼2.6 kpc; Jurić et al. 2008) from the others. The only thing that changes in these simulations is the strength of the Galactic tide and the encounter rate of passing field stars. These are modulated by decreasing or increasing the local density according to $\rho =\rho _{0}e^{r/r_L}$, where r_L is the Milky Way disk's scale length and ρ₀ is the approximate density of stellar matter near the Sun's galactocentric distance (∼0.04 M_☉ pc⁻³). When modeling external perturbations in other regions of the Milky Way's disk we do not vary the stellar dispersion. As shown in Sections 3.1 and 4.1, as long as stellar encounters remain in the impulsive regime, the collision probability within wide binaries is not strongly dependent on field star encounter velocities. The physical reason for this is that although an increase in dispersion increases the total number of stellar encounters, they become weaker on average, since their encounter timescales are shortened by the increased encounter velocities.

When stars pass very near one another, we must account for the dissipation of orbital energy via dynamic tides. Our understanding of this process remains very uncertain. Consequently, we choose to employ a very simple treatment of this dissipation for each periastron passage (Press & Teukolsky 1977; Lee & Ostriker 1986). In this model, the secondary star deposits energy into oscillations of various spherical harmonics of index l within the primary:

$$\begin{equation} \Delta E_{l} = \left(\frac{Gm_{1}^{2}}{R_1}\right)\left(\frac{m_{2}}{m_{1}}\right)^2\left(\frac{R_1}{q}\right)^{2l+2}T_{l}(\eta), \end{equation} \tag{ 1 }$$

where m is stellar mass, R is stellar radius, q is periastron, and the subscripts refer to the primary and secondary star. T_l is a dimensionless parameter that depends on the stellar masses and radii as well as q. The values of T_l are of the same order for different l as long as the periastron and stellar parameters are fixed. Thus, as long as q/R₁ ≳ 3 there is less than a ∼1% error if we neglect the higher order terms and only calculate the l = 2 (quadrupole) and l = 3 (octopole) terms. In our simulations, T₂ and T₃ are calculated during each periastron passage and the resulting energy change is applied instantaneously at periastron passage. (The energy that the primary star deposits into the secondary is also calculated and applied of course.) Our dissipation calculations are incorrect for q/R₁ < 3, but as Figure 3 shows, none of these systems ever actually collide since tidal circularization begins near q/R₁ ≃ 5.

The values of the T_l terms also depend on the assumed internal structure of the star. The fiducial model we employ assumes all stars to be n = 3 polytropes. Because n = (3/2) polytropes are more applicable to low-mass stars, we also rerun our simulations assuming these polytropes instead (Lee & Ostriker 1986). As we will see in our results, our collision rates are nearly unaltered between these two models. The reason is that the modest variation in T_l values for different polytropes is overwhelmed by the steep drop of the quadrupole and octupole terms with q. We refer readers to the original texts describing these models for the details in calculating T_l terms (Press & Teukolsky 1977; Lee & Ostriker 1986). It should be noted that the n = 3 polytrope model in Press & Teukolsky (1977) contains arithmetic errors corrected in Lee & Ostriker (1986).

Our binaries are integrated using a simple symplectic scheme for 10 Gyr, roughly the age of the Milky Way thin disk. In this kick–drift–kick scheme, the secondary is drifted along a Keplerian orbit about the primary for a time step τ (typically τ = 550 yr). In between these drift time steps, the acceleration due to the Galactic tide is applied as well as any impulses from passing field stars. As stated previously, we run batches of 10³ binary systems, all sharing the same primary star and external perturbations. This scheme of shared primaries and stellar encounters greatly enhances our computing efficiency, with only a marginal loss in statistics since we still run 100 different batches of 10³ binaries for each galactic position.

In our simulations we define a stellar collision at any point where a binary passes through periastron, and the periastron value is below the combined physical radii of the primary and secondary stars. Although stellar radii vary with mass, metallicity, and age, there is an approximate linear one-to-one relationship between stellar mass and radius between 0.1 and 1.0 M_☉ (Bayless & Orosz 2006). Based on this, we assume that our simulated stars have radii according to the following empirical relationship:

$$\begin{equation} R_{*}(R_{\odot }) = 0.0324 + 0.9343m_{*}(M_{\odot }) + 0.0374m_{*}^{2}(M_{\odot }). \end{equation} \tag{ 2 }$$

These radii are used in calculating collision rates as well as calculating the tidal dissipation described in previous sections.

In most of our simulations we keep the local stellar dispersion fixed when generating stellar passages in our simulations. In reality, the dispersion increases nearer to the Galactic center and decreases further away, and the typical field star encounter velocity can actually vary by over an order of magnitude throughout the Milky Way disk (Veras et al. 2014). However, the exact dependence of local stellar dispersion on galactocentric distance is very uncertain and model-dependent (Brasser et al. 2010). In Section 4.1 we will demonstrate that the rate of periastron change in a very wide binary is independent of the local stellar dispersion as long as stellar encounters remain in the impulsive regime.

Here we compliment this with an additional set of numerical experiments showing that the rate of stellar collisions within very wide binaries is relatively insensitive to the local stellar dispersion and is instead mostly set by the local stellar density. Modeling the formation of the Oort Cloud at different galactocentric distances, Brasser et al. (2010) assume that the local stellar velocity dispersion scales as (8/R_G)^0.92, where R_G is the galactocentric distance measured in kpc. In our additional simulations, we employ this model, scaling the Rickman et al. (2008) stellar encounter dispersions as a function of the galactocentric distance of our simulated binaries. This scaling causes the local stellar dispersion to vary by over a factor of four between our simulated binaries closest to the Galactic center and our outermost simulated binaries.

The collision rates seen in these new simulations are compared against our fiducial simulations in Figure 4 (both sets of simulations contain the same tidal dissipation prescription). As can be seen, the collision rates as a function of galactocentric distance in both sets of simulations are very similar. The simulations with a variable local stellar dispersion may have a slightly enhanced encounter rate in the interior of the Galaxy, but this variation is dwarfed by the sensitivity to tidal dissipation seen in Figure 2. Thus, we can conclude that the stellar collision rates determined in our work are fairly independent of our assumptions of the typical field star encounter velocity.

Zoom In Zoom Out Reset image size

Figure 1(B) shows that instead of colliding, many binary companions pass close enough to each other that their orbital energy is rapidly dissipated by dynamic tides. In this example, a binary's semimajor axis is shrunk from 8000 AU to less than 100 AU in just seven periastron passages. Unlike the illustrative case displayed in Figure 1(B), we normally stop following the evolution of binaries with a < 100 AU for computing efficiency reasons. In addition, at small a our technique of applying an instantaneous tidal dissipation at periastron becomes questionable, and binaries also transition from the dynamic tidal regime to the equilibrium tidal regime at small a. However, these binaries' semimajor axes will continue to shrink below what is modeled in our typical simulations. This raises the prospect that some very wide binaries can eventually be transformed into very tight contact binaries.

A simple estimate of the encounter distance necessary for this transformation can be made calculating the rate that dynamic tides can change a binary's semimajor axis. In our idealized simulations, the energy dissipation due to dynamic tides, ΔE(q), is applied instantaneously at periastron passage. The level of dissipation in our fiducial tidal model depends only on stellar mass, radius, and periastron and is independent of semimajor axis. The energy dissipation rate is therefore given by the orbital period of the binary if the periastron remains constant. This yields the following rate of change in semimajor axis:

$$\begin{equation} \dot{a} = \frac{1}{\pi m_{1} m_{2}} \left(\frac{m_{1} + m_{2}}{G}\right)^{1/2}a^{1/2} \Delta E(q). \end{equation} \tag{ 3 }$$

From this expression, we can then calculate the time required for dynamic tides to transform a binary from an initial semimajor axis (a_i) to a final semimajor axis (a_f):

$$\begin{equation} t_{c} = \frac{2\pi m_{1} m_{2}}{\Delta E(q)}\left(\frac{G}{m_{1} + m_{2}}\right)^{1/2} \left(a_{f}^{1/2} - a_{i}^{1/2}\right). \end{equation} \tag{ 4 }$$

If a_i ≫ a_f we can ignore the a_f term on the righthand side of Equation (4) and obtain a simple estimate of the tidal circularization timescale.

Assuming a "typical" very wide binary configuration of m₁ = m₂ = 0.4 M_☉ and a_i = 10⁴ AU, we calculate this circularization time as a function of binary periastron for our fiducial tidal dissipation model (Press & Teukolsky 1977; Lee & Ostriker 1986) in Figure 5. This plot allows us to estimate the maximum contribution that very wide binaries can make to the population of compact binary systems in our Galaxy. We see that for q beyond 10 stellar radii the circularization timescale is greater than the age of the Galaxy. This therefore requires periastron to be below ∼7 stellar radii for tidal circularization to be possible. However, this estimate also assumes that periastron is fixed. The minimal importance of a_f in Equation (4) demonstrates that most of the circularization time is spent at large a. Consequently, if t_c is large there is plenty of time for field star encounters to alter the periastron during the tidal circularization process. Sometimes these encounters could push the periastron inward, speeding up tidal circularization, but others could pull the periastron beyond seven stellar radii and freeze the process. Figure 3 shows that no stars passing within ∼5 stellar radii of each other later go on to collide. This suggests that for $q\lesssim 5\bar{R_{*}}$ tidal circularization proceeds fast enough that periastron will remain nearly constant during this process. We can therefore consider $q\simeq 5\bar{R_{*}}$ to be the periastron threshold inside which tidal circularization efficiently takes place.

Zoom In Zoom Out Reset image size

As mentioned previously, the intent of our simulations is to study star–star collisions, and they are not fully equipped to model the evolution of binaries down to very tight configurations. However, we can estimate the maximum contribution using our numerical results. To do this, we assume that every binary that tidally evolves to a ⩽ 100 in our simulations is destined to become a contact binary. Based on this, we can estimate the maximum probability that very wide binaries evolve into contact binaries.

We calculate these probabilities in Figure 6. We see that the contact binary formation probability is of course dependent on the dynamic tidal dissipation we employ in our simulations as well as the local stellar density, which controls how often binaries are cycled into very eccentric orbits. For our fiducial tidal model (Press & Teukolsky 1977; Lee & Ostriker 1986), we find that there is a median probability of 0.17% that a very wide binary will evolve into a contact binary. If we then optimistically assume that 10% of stars reside in very wide binaries, we predict that 1 in ∼6000 stars will reside in contact binaries through our mechanism. In contrast, modern surveys of contact binaries find that there is one contact binary for every 300–500 main sequence stars in the solar neighborhood (Rucinski 2002; Gettel et al. 2006). In our fiducial tidal model case, our mechanism falls about an order of magnitude short of explaining the contact binary fraction even with optimistic assumptions. In fact, even our extreme tidal model, which almost certainly overestimates the degree of tidal circularization, cannot explain the observed contact binary fraction. Although very wide binaries likely contribute to the formation of contact binaries, there must be more efficient formation mechanisms such as Kozai oscillations with tidal friction within triple star systems (Kiseleva et al. 1998; Eggleton & Kiseleva-Eggleton 2001; Fabrycky & Tremaine 2007). Nevertheless, we find it very surprising that there exists a potential pathway for the very widest known binaries to evolve into the very tightest binaries.

Zoom In Zoom Out Reset image size

3.2.1. Equilibrium Tidal Evolution

As stated previously, we typically only consider the effect of dynamic tides (Lai 1997; Ivanov & Papaloizou 2004) on the evolution of wide binaries. We apply an instantaneous kick in energy at periastron passage in the relevant binaries (see Section 2.3). This approximation is valid for very wide orbits with eccentricities close to unity, but breaks down as the stars approach each other. In that situation, we must consider the regime of equilibrium tidal evolution. Here we do so for a few illustrative systems.

Equilibrium tides account for the deformation of each of the binary stars by the other's gravitational potential. The torques associated with tidal bulges, combined with energy dissipation within each star, act to change the binary orbit (Darwin 1879; Hut 1981; Eggleton et al. 1998). Equilibrium tides act to shrink and circularize very eccentric binary orbits like the ones we consider.

The total angular momentum is conserved during tidal evolution. Assuming the orbital angular momentum L to be the dominant component of the total angular momentum budget, a wide binary that is tidally interacting will be transformed into a close binary on a circular orbit with a semimajor axis of twice the periastron distance q. This is because $L \propto \sqrt{a (1-e^2)} = \sqrt{q(1+e)}$ and e ≈ 1 for these binary orbits.

We calculate the late-stage evolution of a set of tidally interacting wide binaries using an integrator employing an equilibrium tide model (Bolmont et al. 2012). The code uses the so-called constant time-lag model (Hut 1981; Eggleton et al. 1998). The binaries are composed of equal-mass components with masses of 0.4 M_☉ and corresponding radii of 0.412 R_☉. The stars are attributed dissipation values assuming pseudo-synchronous rotation with small obliquities (Hansen 2010; Bolmont et al. 2012), although we note that stellar dissipation rates are poorly constrained. We simulate nine binaries, each with an initial semimajor axis of 99 AU but initial periastron distances from 2.5 to 9 stellar radii. As in Figure 1(B), we assume that the transition from dynamic to equilibrium tides occurs at a =100 AU. In reality, there will be a smoother transition between the two regimes.

Figure 7 shows the outcome of the integrations. As expected, the orbits of the binaries are circularized at twice the initial periastron distance. The timescale for the evolution depends sensitively on q₀. Whereas the initially near-grazing binary with q₀ = 2.5R_* has its orbit circularized in ∼10⁵ yr, binaries with q₀ > 6–7R_* require hundreds of millions to billions of years for their orbits to be circularized. Of course, this evolution is strongly dependent on the stellar parameters, in particular the dissipation rates within each star (Hut 1981; Eggleton et al. 1998). Nonetheless, Figure 7 gives a reasonable order of magnitude estimate for the equilibrium tidal evolution timescales of these binaries.

Zoom In Zoom Out Reset image size

In our simulations, we model the orbital evolution of binaries with semimajor axes between 10³ AU and 3 × 10⁴ AU. Here we calculate the collision probability of a very wide binary as a function of its semimajor axis. Much of this derivation is nearly identical to the calculation of the flux of comets into the inner solar system from various regions of the Oort Cloud (Hills 1981). The reason these problems are so similar is because orbital perturbations from the local galactic environment drive both phenomena. Further, instead of the Jupiter–Saturn barrier that prevents some eccentric comets from reaching the inner solar system, our binary systems have a dynamic tidal barrier whose dissipation prevents some eccentric binaries from ever reaching a collisional state.

To begin, let us consider an ensemble of N very wide binaries. To simplify our problem, let us initially assume that these binaries all have the same age, stellar radius, and semimajor axis. There are three ways to lose systems from this ensemble: collision, tidal circularization (see Figure 1(B)), and ionization from stellar impulses (which we will ignore for the moment). We have seen from our simulations that very wide binaries are lost via tidal circularization if the stars passed within ∼5 stellar radii of each other in our fiducial model. For a thermalized eccentricity distribution, the fraction of orbits in our binary ensemble with periastron less than five stellar radii is

$$\begin{equation} F_{q} = \frac{10R_{*}}{a}, \end{equation} \tag{ 5 }$$

where a is the binary semimajor axis and R_* is the typical stellar radius. We will call this set of orbits the loss cone, since very wide binaries are lost through tidal circularization or potentially collisions if they occupy this orbital space. Again if we assume that orbits are isotropized, then two-fifths of these binaries will actually collide rather than tidally circularize (since a periastron of two stellar radii or less yields a collision). If binary orbits were distributed isotropically but did not evolve under external perturbations, this would be the total fraction of our binaries lost as a function of a. However, these binary orbits are not static. Consequently, new losses take place every orbital revolution as binaries with large initial periastron are perturbed to low periastron states.

Like the Oort Cloud, the pericenter changes in our binaries are driven by perturbations from the Galactic tide and other passing field stars. However, in our binary systems we are concerned with extremely tiny periastron values ∼3 orders of magnitude smaller than typical long-period comet pericenters. In comet dynamics, the Galactic tide plays a major role in injecting comets into the inner solar system. However, the periastron change per orbit due to the Galactic tide approaches zero as q approaches zero (Levison et al. 2006). Meanwhile, the tidal barrier requires relatively large shifts in periastron (Δq ∼ q) for a collision to occur. This renders the Galactic tide ineffective in producing stellar collisions. Thus, we only need to consider the influence of impulses from passing field stars to estimate the fraction of binaries that collide.

In the impulsive regime, the change in velocity a passing field star imparts on a binary star is approximately (Hills 1981)

$$\begin{equation} \Delta v = |\Delta v_{1} - \Delta v_{2}| = \frac{2Gm_{*}}{v_{*}}\left(\left|\frac{b_{1} - b_{2}}{b_{1}b_{2}}\right|\right), \end{equation} \tag{ 6 }$$

where v is velocity, b is the closest approach distance of the passing star, and m_* is the passing star mass. The subscripts 1, 2, and * refer to the primary, secondary and passing star respectively. Averaged over all binary orbit orientations, <|b₁ − b₂| > ∼r and $\left<b_{1}b_{2}\right>\sim b_{1}^{2} \equiv b_{*}^{2}$. Further, if e ∼ 1 then <r > ∼1.5a. Therefore,

$$\begin{equation} \left<\Delta v\right> \simeq 3Gam_{*} / v_{*}b_{*}^{2}. \end{equation} \tag{ 7 }$$

Meanwhile, the orbital velocity of the binary companion at a distance <r > is just

$$\begin{equation} v_{b} = \left(G\mu / 3a\right)^{1/2}, \end{equation} \tag{ 8 }$$

where μ is the reduced mass of the binary.

On average, a stellar encounter of a given impact parameter and velocity will change the orbital velocity of a binary by Δv. Again, on average the magnitude of Δv should be the same for all binaries if they have the same semimajor axes. The direction of Δv, however, will depend on the orientation of the binary relative to the passing star. Thus, the initial velocity vectors of an ensemble of randomly oriented binaries will be spread throughout a smear cone by a stellar passage (where the cone is centered upon the original velocity vector). The angular size of this smear cone relative to the initial velocity direction is just

$$\begin{equation} \theta = \frac{\left<\Delta v\right>}{v_{b}} = 9\frac{\mu v_{b}}{m_{*}v_{*}}\left(\frac{a}{b_{*}}\right)^{2}. \end{equation} \tag{ 9 }$$

If Δv ≪ v_b then the fraction of all velocity space taken up by the smear cone is

$$\begin{equation} F_{s} = \frac{\pi \theta ^2}{4\pi } = \frac{27}{4}\left(\frac{m_{*}}{\mu }\right)^{2}\left(\frac{a}{b_{*}}\right)^{4}\left(\frac{G\mu }{av_{*}^{2}}\right). \end{equation} \tag{ 10 }$$

Meanwhile, for an isotropic velocity distribution we have already calculated the fraction of velocity space occupied by orbits leading to a collision or circularization in Equation (5). Thus, we can now calculate the size ratio of the smear cone to the loss cone as

$$\begin{equation} \frac{F_{s}}{F_{q}} = \frac{27}{40} \left(\frac{m_{*}}{\mu }\right)^{2} \left(\frac{a}{b_{*}}\right)^{4}\left(\frac{G\mu }{R_{*}v_{*}^{2}}\right). \end{equation} \tag{ 11 }$$

This ratio tells us how much of the loss cone is refilled with new binaries from a single stellar passage.

If a loss cone is always filled (F_s/F_q = 1) for an ensemble of binaries with a given a, then the fractional loss rate within this ensemble of binaries will just be set by the orbital period of the binaries and will be

$$\begin{equation} \dot{L} = F_{q} \sqrt{\frac{G\left(m_{1} + m_{2}\right)}{4\pi ^{2}a^{3}}}. \end{equation} \tag{ 12 }$$

Thus, the loss rate (and therefore collision probability) decreases with increasing binary semimajor axis when the loss cone is always filled. As the binary semimajor axis is decreased, the loss rate increases initially since the orbital period is shorter and F_q is larger for tighter binaries. However, tighter binaries also become less sensitive to field star perturbations, which are responsible for refilling the loss cone. As a result, below some critical binary semimajor axis the loss cone will not always be filled, and this will act to decrease the loss rate and probability.

We must therefore calculate how often stellar passages that are capable of refilling the loss cone for a given binary semimajor axis are expected to occur. The required stellar passage strength can be calculated by setting F_s/F_q to 1. Thus, to refill the loss cone Equation (11) requires

$$\begin{equation} (b_{*}^{2}v_{*})^{2} \le \frac{27}{40} \left(\frac{m_{*}}{\mu }\right)^{2}\left(\frac{G\mu a^{4}}{R_{*}}\right). \end{equation} \tag{ 13 }$$

Fortunately, the rate of stellar encounters with $b_{*}^2v_{*}$ less than a given quantity is very easy to calculate for a given stellar environment. This rate is given by

$$\begin{equation} f = \pi (b_{*}^{2}v_{*}) n_{*}, \end{equation} \tag{ 14 }$$

where n_* is the local spatial density of stars.

As f gets larger, the binary loss rate increases until f exceeds the orbital frequency. At this point the loss cone is "saturated," and further increases in f will not yield more collisions. Thus, binaries with the highest probability of collision are those with semimajor axes a_crit where f exactly equals the orbital frequency. Binaries with larger a will have lower collision rates because of their longer orbital period and saturated loss cones. Binaries with smaller a will have lower collision rates because their loss cone is only occasionally filled.

Setting f equal to the orbital frequency, we can now solve for the stellar encounter strength that fills the loss cone once per orbital period for binaries with critical semimajor axis a = a_crit

$$\begin{equation} (b_{*}^{2}v_{*})^{2} = \frac{1}{4\pi ^{4}n_{*}^{2}}\frac{G\left(m_{1}+m_{2}\right)}{a_{{\rm crit}}^{3}}. \end{equation} \tag{ 15 }$$

Now that we have re-expressed the stellar encounter strength in terms of the local stellar density and a_crit we can substitute this expression into Equation (13) to obtain

$$\begin{equation} a_{{\rm crit}} = \left[\frac{10}{27\pi ^{4}}\frac{R_{*}}{n_{*}^{2}}\frac{\mu \left(m_{1}+m_{2}\right)}{m_{*}^{2}}\right]^{1/7}. \end{equation} \tag{ 16 }$$

Thus, the binary separation with the maximum loss (and therefore collision) probability decreases quite slowly with the local density of field stars ($a_{{\rm crit}}\sim n_{*}^{-2/7}$). Moreover, this critical separation does not depend at all on the velocity distribution of field stars as long as stellar encounters remain in the impulsive regime.

Once a_crit is determined it is easy to calculate the maximum loss probability that a binary can have for a given stellar environment and mass configuration. This is just determined by the orbital period of a_crit, the size of the loss cone, and the age of the binary, T_age:

$$\begin{equation} L = 1 - (1 - F_{q})^{\frac{T_{{\rm age}}}{2\pi }\sqrt{\frac{G\left(m_{1} + m_{2}\right)}{a_{{\rm crit}}^{3}}}}. \end{equation} \tag{ 17 }$$

We can also calculate how the binary loss probability changes for binary semimajor axes greater or less than a_crit. When a < a_crit the limiting factor is how often the loss cone is refilled. In this regime loss cone refilling happens less than once per orbital period, so this portion of phase space will quickly empty and remain empty until another stellar encounter refills it. Thus,

$$\begin{equation} L_{a \; < \;a_{{\rm crit}}} = 1-(1-F_{q})^{T_{{\rm age}}f}, \end{equation} \tag{ 18 }$$

where f is given by Equations (13) and (14). As long as TfF_q ≪ 1 this probability can be approximated by

$$\begin{equation} L_{a \; < \;a_{{\rm crit}}} = T_{{\rm age}}fF_{q}. \end{equation} \tag{ 19 }$$

As mentioned above, two-fifths of the binaries that are lost through our loss cone will do so via collision. Consequently, we only have to multiply Equation (19) by 2/5 to obtain a collision probability P for binaries with a < a_crit:

$$\begin{equation} P_{a \; < \;a_{{\rm crit}}} = T_{{\rm age}}m_{*}n_{*}a\sqrt{\frac{54 \pi ^{2} G R_{*}}{5\mu }}. \end{equation} \tag{ 20 }$$

Meanwhile, if a > a_crit then the loss cone is always filled and the limiting factor is how often the binaries pass through periastron, i.e., the orbital period, T_orb. In this case, Equation (17) yields the fraction of binaries that collide for a given a:

$$\begin{equation} P_{a \; > \; a_{{\rm crit}}} = \frac{2}{5} [1-(1 - F_{q})^{T_{{\rm age}} / T_{{\rm orb}}}]. \end{equation} \tag{ 21 }$$

Again, if T_ageF_q/T_orb ≪ 1 then this probability is well-approximated by

$$\begin{equation} P_{a \; > \; a_{{\rm crit}}} = \frac{2}{5}\frac{T_{{\rm age}}}{T_{{\rm orb}}}F_{q} = \frac{2T_{{\rm age}}}{\pi }\sqrt{\frac{G\left(m_{1}+m_{2}\right)R_{*}^{2}}{a^{5}}}. \end{equation} \tag{ 22 }$$

Equations (20) and (22) are of course equal for a = a_crit.

Thus, the collision probability of any binary can be estimated with Equations (16), (20), and (22). These expressions can be further simplified if we assume all stars (including passing stars and our binary companions) to have some characteristic mass, m_*. In this case, m₁ = m₂ = 2μ = m_*. With this assumption the expression for the binary semimajor axis with the highest collision probability becomes

$$\begin{equation} a_{{\rm crit}} = \left(\frac{10}{27\pi ^{4}}\frac{R_{*}}{n_{*}^{2}}\right)^{1/7}, \end{equation} \tag{ 23 }$$

while the probability of collision for any binary semimajor axis can be found with

$$\begin{equation} P_{a \; < \;a_{{\rm crit}}} = T_{{\rm age}}n_{*}a\sqrt{\frac{108 \pi ^{2} G m_{*}R_{*}}{5}} \end{equation} \tag{ 24 }$$

and

$$\begin{equation} P_{a \; > \;a_{{\rm crit}}} = \frac{2T_{{\rm age}}}{\pi }\sqrt{\frac{2Gm_{*}R_{*}^{2}}{a^{5}}}. \end{equation} \tag{ 25 }$$

An example of this collision probability is plotted in Figure 8. In this calculation, we set m_* = 0.4 M_☉ (the approximate mean field stellar mass), R_* = 0.4 R_☉, T = 10 Gyr, and n_* = 0.1 pc⁻³ (the approximate spatial density of stars near the Sun). We see in Figure 8 that the collision probability increases linearly with the binary semimajor axis until a = a_crit, at which point the collision probability is roughly 1 in 300. Beyond this, the collision probability falls off quickly with a^−5/2.

Zoom In Zoom Out Reset image size

This probability distribution can then be integrated to yield a median semimajor axis for binary collisions as well as a total probability of binary collisions. We have generated a distribution like the one shown in Figure 8 for each of the galactic environments we simulate. Then we integrate these distributions over log-a space (since we assume a log-uniform distribution of binary semimajor axes). The median semimajor axis of colliding binaries is found where the normalized integral equals 0.5. This predicted median semimajor axis is plotted as a function of local field star density in Figure 9. In addition, we include the median binary semimajor axis found among our numerically simulated colliding systems in each environment. As can be seen, the median semimajor axes found in both our numerical experiments and our theoretical predictions increase similarly with decreasing local stellar density.

Zoom In Zoom Out Reset image size

In a similar fashion, we also integrate our probability distributions over all semimajor axes (in log-a space) to obtain a predicted total collision probability for our binaries for each galactic environment. These probabilities are plotted in Figure 10. Here we see that our numerical results yield significantly lower collision probabilities than that predicted by Equations (23)–(25). Moreover the disagreement becomes worse with denser environments. A major reason for this discrepancy is that the above analysis neglected the ionization of binaries by field star impulses. This acts to lower the collision probability because the binary lifetime can be substantially shortened via ionization, giving the binary a smaller time window to produce a collision. Further, this effect is most pronounced in dense environments, where binaries are ionized more quickly.

Zoom In Zoom Out Reset image size

We can attempt to account for this ionization by using the approximate formula for the ionization half-life of very wide binaries subject to stellar passages (Bahcall et al. 1985)

$$\begin{equation} t_{1/2} = 0.00233 \frac{v_{*}}{Gm_{*}n_{*}a}. \end{equation} \tag{ 26 }$$

Assuming binary ionization is governed by exponential decay, Equations (24) and (25) can be modified to yield binary collision probabilities that account for the ionization of binaries. These modifications yield the following new expressions for collision probability

$$\begin{equation} P_{a \; < \;a_{{\rm crit}}} = \tau n_{*}a\sqrt{\frac{108 \pi ^{2} G m_{*}R_{*}}{5}} (1-e^{-T_{{\rm age}}/\tau }) \end{equation} \tag{ 27 }$$

and

$$\begin{equation} P_{a>a_{{\rm crit}}} = \frac{2\tau }{\pi }\sqrt{\frac{2Gm_{*}R_{*}^{2}}{a^{5}}} (1-e^{-T_{{\rm age}}/\tau }), \end{equation} \tag{ 28 }$$

where τ is the decay time constant (t_1/2/ln 2). In Figure 10, we assume v_* = 40 km s⁻¹ for the typical stellar encounter velocity and overplot the total collision probabilities predicted by our newly modified expressions. There is a marked improvement in our predicted collision probabilities, especially for dense stellar environments. The new predicted collision probabilities are now all within a factor of two of our numerically measured collision probabilities. However, even this treatment does not fully account for the ionization process. Since it is a slow diffusive process, binaries can diffuse out to larger a over the course of our simulations without actually being ionized, and Figure 8 indicates that this diffusion to large a will significantly decrease the collision probabilities of binaries that remain bound. Nonetheless, we consider Figure 10 to be a very good match between simulations and theory given that all binary and passing stars are assumed to have the same fixed mass and that the stellar encounter velocity is also assumed to be fixed.

As mentioned previously, the Milky Way's tide is ultimately ineffective in driving the periastron of binaries to collision. The final decrease in periastron during a collision is always caused by a stellar encounter. However, the Galactic tide can still play a role in pushing binaries to the brink of collision before stellar encounters take over. The direction of periastron change due to the Galactic tide is set by the value of the argument of pericenter, ω, that the binary has with respect to the Milky Way midplane. For most galactocentric distances (including the ones modeled here), the vertical disk term dominates over the radial tide component (Brasser et al. 2010). For binary eccentricities near 1 (as in the case of our colliding binaries), the precession rate of ω under the vertical component of the Milky Way's tide in a disk-only potential is approximately given by

$$\begin{equation} \dot{\omega }_{GT} \simeq -\frac{5\sqrt{2G}\pi }{q\left(m_{1} + m_{2}\right)}\rho _{0}a^2 \sin ^{2}{\omega }\cos ^{2}{i}, \end{equation} \tag{ 29 }$$

where q is the periastron and i is the inclination relative to the Galactic midplane (Heisler & Tremaine 1986; Levison et al. 2006).

In addition to the Galactic tide, there are other sources of orbital precession. These arise from the stars' tidal bulges as well as general relativistic effects. If these sources of precession become comparable to that due to the Galactic tide, the dynamics of very wide binaries can diverge from that seen in our numerical simulations, since these additional precession sources are not included in our numerical work. The precession of ω due to general relativity is given by the following expression (Fabrycky & Tremaine 2007)

$$\begin{equation} \dot{\omega }_{GR} = \frac{3G^{3/2} (m_{1} + m_{2})^{3/2}}{a^{5/2}c^{2} (1-e^{2})}. \end{equation} \tag{ 30 }$$

Meanwhile, the precession due to the tidal bulges can be estimated from equilibrium tidal models to be (Fabrycky & Tremaine 2007)

$$\begin{eqnarray} \dot{\omega }_{{\rm tide}} &=& \frac{15 [G(m_{1} + m_{2})]^{1/2}}{8a^{13/2}}\frac{8+12e^{2}+e^{4}}{(1-e^{2})^{5}} \nonumber\\ &&\times \left(\frac{m_{2}}{m_{1}}k_{1}R^{5}_{1}+\frac{m_{1}}{m_{2}}k_{2}R^{5}_{2}\right), \end{eqnarray} \tag{ 31 }$$

where k is the Love number, R is the stellar radius, and the subscripts 1 and 2 refer to the primary and secondary stars as usual.

To estimate the accuracy of our simulations we can now compare the magnitudes of the three precession rates given in Equations (29)–(31). This is done in Figure 11. To generate this plot we assume m₁ = m₂ = 0.4 M_☉, R₁ = R₂ = 0.4 R_☉, and k₁ = k₂ = 0.014 (a typical estimate for stars (Fabrycky & Tremaine 2007)). In addition, we set q = 5R_* = 2 R_☉, since we already know from our numerical work that dynamic tides dominate over Galactic tidal effects inside of this periastron distance. Thus, Figure 11 measures whether binary orbital precession is dominated by the Galactic tide right up until binary companions collide or tidally circularize.

Zoom In Zoom Out Reset image size

Figure 11(A) demonstrates that precession due to the Galactic tide is at least an order of magnitude faster than that due to tidal bulges for almost any very wide binary orbit. The exception is for binaries with a ≲ 3000 AU that are located in very low density environments. Here the two precession rates are comparable when binary periastron is at 5 R_* (near the dynamic tidal barrier). Figures 8–10 clearly show, however, that this region of parameter space accounts for only a very small fraction of the binary collisions seen in our simulations. Similarly Figure 11(B) shows a nearly identical situation when comparing general relativistic precession with galactic tidal precession. Only at very low galactic densities and a near 10³ AU are the two precession rates comparable. Again, this is only considering binary orbits with q = 5 R_*. For more circular orbits, the Galactic tide will be much more dominant. Thus, we can conclude that for the vast majority of the parameter space probed in our numerical simulations it is perfectly acceptable to ignore the orbital precession due to tidal bulges and relativity.