BROWN DWARF BINARIES FROM DISINTEGRATING TRIPLE SYSTEMS

It is an interesting fact of nature that the motion of two isolated bodies is completely deterministic, while the addition of more bodies, even just one, can render the motion completely chaotic. While the gravitational two-body problem was solved already by Newton (1687), the three-body problem has been the subject of major efforts for more than 300 yr (e.g., Euler 1772; Lagrange 1778; Jacobi 1836; Hill 1878; Poincaré 1892-1899). A formal solution was not discovered until the work of Sundman (1912), which unfortunately was of no practical use in the calculation of orbits. About a hundred years ago, the modern era began with the first explorations of orbital integration, where each body is moved in small steps. Burrau (1913) studied the famous Pythagorean problem in this manner. Orbital integration lends itself extremely well to modern computational techniques, which have revolutionized our understanding of the three-body problem. For a detailed discussion, see Valtonen & Karttunen (2006) and Aarseth et al. (2008), as well as Valtonen & Mikkola (1991) for astronomical applications.

Many subsequent studies have defined the broad characteristics of the motion of three bodies that are initially in a non-hierarchical configuration (e.g., Anosova 1986). There are essentially three states of a triple system. Initially most of the time is spent in interplay, during which the three bodies move chaotically with no periodicity. Interplay is interspersed with close triple approaches, in which the three bodies are briefly brought close together, and during which energy and momentum can be exchanged. The third state is the ejection, which can only occur immediately following a close triple approach, when one body (usually the lightest) is ejected, while the two remaining bodies provide the energy of ejection and as a result form a tighter bound binary system in an orbit that is usually highly eccentric. The ejected member may move in an approximately elliptical orbit, in which case it returns to the binary for more interplay or another ejection. Or it may move in a hyperbolic orbit, in which case the ejection leads to an escape. So, an ejection does not necessarily imply an escape, it can be into a bound distant orbit or into an escape. Figure 1 shows an example of the chaotic motions of the interplay phase, followed by a close triple encounter and subsequent ejection into an escape. The motion of the three bodies in a non-hierarchical triple system is very sensitive to even small changes in masses.

Zoom In Zoom Out Reset image size

In this paper, we explore the role of a cloud core for the evolution of a newborn triple system. In particular, we allow the three bodies to accrete as they move around the cloud core, and this has profound effects on the dynamical evolution of the system (e.g., Bonnell et al. 1997; Umbreit et al. 2005; Reipurth et al. 2010). Figure 2 shows in broad outline the possible outcomes of dynamical evolution of three identical stellar embryos moving inside and accreting from a cloud core. Either one body falls to the center of the core, rapidly growing in the process, while dynamically keeping the two other bodies toward the outskirts of the core. If a bound triple system results from this, it takes one of the two forms shown to the upper right or lower left of the figure. Or two bodies control the center of the cloud core, keeping the third body at bay, and a resulting bound triple system will then take one of the two forms to the upper left or the lower right of the diagram. For a review of multiplicity among newborn stars, see Reipurth et al. (2014).

Zoom In Zoom Out Reset image size

Finally we should emphasize that we are here studying isolated triple systems, which are thus unaffected by the presence of additional bodies. The numerical results presented here focus on internal instabilities in young triple systems, and are therefore best compared to observations of BD binary populations in T Tauri associations and young moving groups rather than in clusters. The same processes are expected to take place in young clusters, but their external perturbations can play an additional significant role, as discussed in Section 4.12.

We have used a code specifically fine-tuned to deal with the problem of three bodies moving inside a cloud core. A total of 200,000 simulations were performed, each for a timespan of 100 Myr. In order to accurately calculate the frequent close encounters, the motion of the three bodies is integrated with the chain regularization method of Mikkola & Aarseth (1993) that allows a precise treatment of the gravitational force. The three bodies are placed randomly inside a cloud core with the structure of a Plummer sphere with radius R having a potential ϕ(r) $=-M/\sqrt{{{r}^{2}}+{{R}^{2}}}$, where R is the radius of the Plummer core (Figure 3). If for a set of three bodies we define q as the ratio of separations of the outer pair (calculated from the most distant body to the mid-point of the remaining two bodies) to the separation of the two closest bodies, then one can calculate the fraction of cases that have q larger than a given number Q. This is illustrated in Figure 4, which shows that the number P(q > Q) asymptotically approaches π/Q³. In order to avoid introducing hierarchical systems as initial conditions for the simulations, the ratios of separations were not allowed to be larger than a Q-value of 5, and from the figure we see that this means less than 2.5% of the randomly chosen initial configurations were rejected. The mean initial separations were chosen randomly between 40 and 400 AU, values consistent with (but still poorly constrained by) observations of embedded young stars. The center of mass of the three bodies was then placed at the origin which also is the center of the gas cloud. In a final step, the initial three-dimensional velocity vectors were randomly chosen for each body and re-scaled so that the virial ratio was 0.5 at the beginning of the simulations, since every initially bound system is in that state at some time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cloud cores all have a radius of 7500 AU, a typical size suggested by observations (Kirk et al. 2006). Core masses were randomly chosen (i.e., probability independent of mass) in the range from 1 to 10 M_⊙. The gas is assumed at rest, so accretion onto a star causes it to slow down, because the linear momentum of the star (mass × velocity) is kept constant. The lack of angular momentum of the accreting material may affect the orbits of the close pairs, which with low angular momentum might shrink further than determined here (Bate 2000; Umbreit et al. 2005) and with high angular momentum gas could even expand (Bate & Bonnell 1997; Bate 2000). Accretion is calculated according to the Bondi–Hoyle prescription, $\dot{M}\propto \rho /{{(c_{s}^{2}+{{v}^{2}})}^{3/2}}$, where v is the velocity of the object relative to the gas, and the sound speed c_s has a value of 0.2 km s⁻¹. The cloud cores lose mass due to accretion by the stars as well as to outflow activity, which is assumed to cause loss of mass from the core by twice the amount that is being accreted. Finally, to simulate the effect of the diffuse interstellar radiation field, the remaining gas disappears linearly with time over a period of 440,000 yr, which is the duration of the Class I phase determined from Spitzer data (Evans et al. 2009). After the gas cloud has disappeared, the slowdown method was used to speed up the computation (Mikkola & Aarseth 1996). Figure 5 shows four examples of the 200,000 simulations including the initial cloud core.

Zoom In Zoom Out Reset image size

Evidently a full hydrodynamic calculation of the gas dynamics would be preferable to this simplified treatment, but this would be prohibitive when performing 200,000 simulations, as we need to obtain good statistics of the complex and chaotic dynamical behavior of triple systems.

We treat the stellar embryos as point sources, with no physical extent. This is evidently a simplification, especially at the earliest evolutionary stages. The main effect, however, is that we can not quantify how often stellar collisions, which certainly occur from time to time, take place.

Stellar masses are chosen from an initial mass function defined by Chabrier (2003, 2005) which has been observationally supported by, e.g., Alves de Oliveira et al. (2012). At the lower end we have truncated this IMF at 0.012 M_⊙, a mass that we here take to represent the dividing line between planets and BDs. This is obviously quite arbitrary, but given the steep decline of the adopted IMF at this mass range, there are very few such "planetary mass objects," so our simulations are rather insensitive to this lower limit. The maximum initial mass is set at 2 M_⊙, since we are only interested in the evolution of low-mass triple systems (Figure 6). Initially all three bodies have been chosen to have identical masses. Nothing is known about the masses of multiple stellar embryos, so there are no empirical constraints on whether they are identical or they differ. However, because of the assumption that the systems are virialized, small bodies have larger velocities and are easily ejected to the outskirts of the cloud core where they have very little opportunity to grow. In other words, if we were to choose embryo masses that were different from each other, we would pre-determine the outcome of the simulations. Only by choosing three identical masses will the final masses clearly reflect differences in their motion through the cloud core (see further details in Section 3). However, if nature chooses to form multiple embryos of different masses, the same processes will take place, but the time-span over which they operate will be shortened (e.g., Anosova 1986).

Zoom In Zoom Out Reset image size

A triple system that is non-hierarchical is inherently unstable, but following a close triple approach and an ejection event it will either disrupt or it will become a bound hierarchical system, which means the system is well described as two elliptic orbits that do not cross each other. For old stellar triple systems, the main threat to their existence is generally considered to be external, from passing stars that may dislodge the outermost component, most frequently due to numerous gentle tugs rather than a single catastrophic event (Ambartsumian 1937; Retterer & King 1982; Weinberg et al. 1987; Jiang & Tremaine 2010). However, the main threat to young triple systems is in fact internal, since the majority of newly formed triple systems are internally unstable. The critical time for a triple system is always the period around periastron passage of the outer component. If at about the same time the inner system is near its apastron passage, then the three stars are closer to each other than at any other time. Let the inner binary components be denoted A and B, and their orbit be numbered 1, while the outer single body is denoted C and its orbit numbered 2. Then the periastron distance q₂ from C to the center of mass of AB is ${{a}_{2}}(1-{{e}_{2}})$. If the q₂ parameter is only a factor of a few larger than the semimajor axis a₁ of the inner system, then each periastron passage of C will cause small perturbations that after many periastron passages eventually may add up to the point that the triple system breaks apart. This may take a long time, and because "a long time" can be defined in many ways, there is no precise definition of stability. The number of systems deemed stable in a population may therefore vary slightly with time. Additional effects like the Kozai resonance further complicates the picture.

Many different formulae to predict longterm stability of a three-body system have been proposed, all of which give slightly different results (e.g., Huang & Innanen 1983; Eggleton & Kiseleva 1995; Valtonen et al. 2008). We here use the criterion proposed by Mardling (2008).

In Figure 7, we show a plot of eccentricity e₂ of the outer orbit against the q₂ parameter relative to the inner semimajor axis a₁. The larger ${{q}_{2}}/{{a}_{1}}$ is, the less dynamical effects will take place when the three stars are closest, and consequently the system is more stable. But since q₂ = ${{a}_{2}}(1-{{e}_{2}})$, then q₂ is dependent not only on a₂ but also on e₂. The larger the eccentricity e₂ is, the smaller becomes q₂, and hence the system becomes more vulnerable to break-up. Physically, this is because the binding energy of the outer orbit is small and the perturbation in energy from the periastron passage may become of the same order of magnitude, making disruption possible. We show the stable and unstable systems, classified according to the Mardling stability criterion, as red and green crosses, respectively, and we see that the distribution follows the behavior expected from the simple arguments above. A more precise statement on the stability boundary can be found in Saito et al. (2012).

Zoom In Zoom Out Reset image size

There are other more subtle influences on the stability of a triple system. The four categories of triple systems resulting from triple evolution with accretion and shown in Figure 2 are not equally stable. In Figure 8(a), (b) we plot the mass of the third body M_c against the mass of the binary ${{M}_{a}}+{{M}_{b}}$ for all triple systems (containing a BD binary) that remain intact at ages of 1 and of 100 Myr, respectively. A diagonal in each figure separates systems which have a dominant single from those that have a dominant binary. In the figure are included only binaries where both components have masses lower than the hydrogen burning limit, hence the cut-off at masses higher than 0.16 M_⊙. The third body can have any mass. At 1 Myr, that is, for newly born triple systems, the bound but unstable systems (green) far outweigh the bound stable systems (red). At 100 Myr, most of the unstable systems have broken up, so the two categories of stable and unstable systems are about equal in size. This is very different from the state at 1 Myr, when by far most of the unstable systems have a dominant binary. This makes intuitive sense, because when the third body is small compared to the binary, then the binary is much more likely to be able to perturb the orbit of the third body when it passes through periastron.

Zoom In Zoom Out Reset image size

It is a rarely appreciated fact, demonstrated by numerical experiments, that the motion of a truly isolated non-hierarchical triple system eventually must lead to an escape, i.e., the end result cannot be a bound triple system. However, bound triple systems are often observed. This is only possible if (1) the bodies are still in an unstable configuration and have not yet disintegrated, or (2) initial conditions accidentally created a rare stable system, or (3) the bodies achieved a stable hierarchical configuration because they have been formed in the presence of another body. Except for special initial conditions, a stable triple system can only form in the presence of a gravitational potential, either from a cloud core or from additional bodies.

Figure 9 shows the time of binary formation for the 200,000 triple systems studied here. The moment a permanent binary forms is also the moment of the last close triple encounter (which in some cases is also the only close triple encounter) when the triple system transforms from a non-hierarchical configuration to a hierarchical one, or the third body escapes. Three curves are shown, one for bound stable triples (red), another for bound but unstable triples (green), and the last for already disrupted systems (gray). The classification is done for all systems at an age of 100 Myr at the end of the simulations. The bulk of stable triples form at ages between 10,000 and 100,000 yr, after which their formation rapidly decreases. After the cloud cores are gone at 440,000 yr (indicated by a vertical dashed line), no further stable triple systems form. The green line shows that bound triples can still form after the cloud cores are gone, but those triples are all unstable, and will sooner or later disrupt.

Zoom In Zoom Out Reset image size

Figure 10 shows the core mass that remains at the time of the last close triple encounter when a binary is formed. For each of the 200,000 simulations, an initial core mass is randomly picked in the range from 1 to 10 M_⊙. So a core may have a low mass at the time of the last triple encounter either because it was low from the beginning, or because the core lost mass prior to that moment. The highest mass cores preferentially form either bound stable triples (red) or the systems break up (gray). For lower and lower masses, more and more of the bound hierarchical triple systems that form are unstable (green), and for core masses less than ∼4 M_⊙ the formation of unstable triples dominate over the stable triples.

Zoom In Zoom Out Reset image size

It is widely assumed that BDs are able to form from low-mass, turbulent-pressure confined cores (Padoan & Nordlund 2004). Fragmentation of such cores naturally would lead to the observed VLM and BD binaries observed, and Jumper & Fisher (2013) show that many observed properties of such low-mass binaries are reproduced by the turbulent core fragmentation model. It thus seems likely that at least some BD singles and binaries are formed this way.

Whitworth & Stamatellos (2006) suggested that if BDs form from very low mass cores, then close BD binaries might form by secondary fragmentation when molecular hydrogen dissociates. This mechanism has not been demonstrated to work in numerical simulations, and while it produces isolated BD binaries, it does not account for the BD binaries found in wide orbits around main-sequence stars. The latter might then be due to capture of isolated BD binaries by stars, but Kaplan et al. (2012) show that this is highly unlikely.

Stamatellos & Whitworth (2009) suggested that BDs could form through fragmentation of extended very massive disks that might exist around Sun-like stars, and that those that remain in orbit around the primary star would be more likely to be in BD binaries than those BDs ejected into the field. Such disks, however, remain to be observed.

In summary, it is likely that BD binaries can be formed by several mechanisms, and it is possible that they can even work together, e.g., a BD binary formed through disk fragmentation can later be dynamically ejected (Bate et al. 2002a).

Larson (1972) advocated a dynamical view of star formation in which small groups of newborn stars interact chaotically. McDonald & Clarke (1993, 1995) suggested that BD binaries do not readily form from purely point-mass dynamical interactions in small gas-free N-body groups—a result confirmed by the numerical simulations of Sterzik & Durisen (1998)—but that dissipative interactions must be involved, as would happen if for example circumstellar disks would exert drag on passing stars.

Sterzik & Durisen (2003) made choices of a clump mass spectrum and a stellar mass spectrum and from a two-step process produced an IMF compatible with the then available observations distributed across few-body clusters, which were followed dynamically, and from which BD binaries were produced. Accretion and gas dynamics were ignored.

We here propose that a common pathway for the formation of BD binaries is through the early breakup of a triple system of stellar embryos, and we explore the statistics and time scales for the related dynamical pathways.

One of the key observational parameters for BD binaries is the binary separation distribution function. In Figure 15 is shown the separation distribution of all ejected BD–BD pairs at ages of 1 and 100 Myr. At 1 Myr there are 9209 ejected BD binaries, and at 100 Myr the number has grown to 15,894. At both ages, the characteristics of the distributions are the same: a steep rise at small separations, followed by a monotonic decrease, reaching zero at separations slightly larger than 200 AU. At 1 Myr, the peak of the separation distribution is at 13 AU, and the median is at 38 AU, the difference reflecting the highly asymmetric distribution seen in the figure. At 100 Myr, the peak of the separation distribution is again around 13 AU, and the median is at 54 AU. Most interestingly, very few BD binaries produced in these simulations have a semimajor axis larger than 100 AU, and none has a semimajor axis larger than 300 AU. The implication is that if a BD binary is observed with a semimajor axis larger than 300 AU, then either it formed in another way, or it must be a triple system where one of the components is an unresolved BD binary (see, however, Section 4.8). The above results compare well with the available observations, see Section 5.2.

Zoom In Zoom Out Reset image size

At 1 Myr, there are 6865 bound (stable and unstable) triple systems where all components are BDs; most of these triples are unstable (see Table 1). After 100 Myr, many of the unstable triples have decayed, and the number of stable and unstable BD triples has declined to just 457. The large majority of these BD triples have outer semimajor axes much larger than 300 AU, and will at first glance appear as wide BD binaries. Wide and very wide BD binaries are therefore not uncommon at young ages, but should be rare in more evolved populations, primarily due to internal disruptions (external disruptions from passing stars will additionally cut the number of wide BD triples as they age). Figure 16 shows the separation distribution function for the outer component of bound triple systems consisting of three BDs on a logarithmic scale, and it is seen that the distribution is flat from a few 100 to 10⁴ AU, that is, across this large interval it follows Öpik's Law, see the discussion in Reipurth & Mikkola (2012).

Zoom In Zoom Out Reset image size

The number of BD binaries with very close separations is significantly smaller compared to the peak value found between 10 and 15 AU. Dynamically, this is equivalent to the BD desert for binaries with a solar-type star and a BD companion (Delgado-Donate et al. 2004). However, the statistics for these very small separations are potentially affected by two shortcomings of our simulations. First, our simulations assume an initial mean separation of the three bodies in the range from 40 to 400 AU. The absence of very close initial mean separations could be responsible for the markedly smaller number of close BD binaries. To search for a dependence between initial mean separation of the newly formed non-hierarchical triple systems with the semimajor axes of the subsequently ejected BD binaries, we plot in Figure 17 the initial mean separation of BD triple systems against the semimajor axes of the ejected BD binaries. The figure illustrates the well known fact that binaries that formed from triple systems generally become tighter in order to gain the energy needed to eject the third body into a distant orbit or an escape. It also shows that there is indeed a statistical dependence of semimajor axes on the initial mean separation, but it is not a strong dependence, and is not likely to be responsible for the "desert" seen at close separations. Second, and perhaps more importantly, our simulations do not take into account that the numerous binaries ejected very early, during the protostellar phase, are likely to harbor circumbinary gas, which will shrink some of their orbits due to viscous interactions during periastron passages, thus producing BD spectroscopic binaries. Additionally, the binaries that remain in the core will shrink due to gas drag (e.g., Stahler 2010). Viscosity thus will move the peak of the separation distribution function to smaller values.

Zoom In Zoom Out Reset image size

An easily observed property of resolved BD binaries is the flux ratio at various wavelengths, which if observed over a large enough wavelength range can give the luminosity ratio of the components. In principle, this can be converted into a mass ratio, although in practice the degeneracy in the mass–luminosity-age relation for BDs makes this a rather uncertain step. Observations of VLM and BD binaries show a strong preference for high mass ratios (e.g., Allen 2007).

In Figure 18 we show the triple diagnostic diagram for the 15,894 triple systems that have disintegrated by 100 Myr, releasing into the field a binary with two BD components. The diagnostic diagram immediately shows two key features of BD binaries formed from ejection: the majority have mass ratios not too different from unity, and there are relatively few BD binaries that have been ejected from triple systems with dominant singles.

Zoom In Zoom Out Reset image size

If we sum across the abscissa, we get the mass ratio distribution in Figure 19. The distribution is striking, showing a preponderance of binaries with mass ratios near unity. No mass ratio is found smaller than 0.3. The mean mass ratio for the 15,894 BD binaries is 0.95, i.e., the components of ejected BD binaries are frequently very alike.

Zoom In Zoom Out Reset image size

This distribution reflects the fact that our stellar embryos start out with identical masses. If two of the embryos quickly bind together, then their motion through the core is almost identical, and they grow by similar amounts, so their mass ratio often remains near one, even if they grow substantially. In the cases where the binary components have different paths through the core and only become bound into a BD binary at a later time, the mass ratio M_B/M_A cannot be less than 0.15 in our simulations, since the lowest mass we consider is 0.012 M_⊙, and the highest mass we consider as a BD is 0.08 M_⊙.

Figure 20 plots the mass ratio as a function of total binary mass. The plot shows that the larger the binary mass, the larger is the possible spread in mass ratios. There are two reasons for this, one imposed by our assumptions, the other due to the shape of the IMF. The lack of small mass ratio binaries at low total binary masses is simple: if we have a primary with mass of, e.g., 0.012 M_⊙, there are no bodies in our simulations with a smaller mass, so in that case the default mass ratio is 1. A similar but weaker and subtler effect comes from the fact that the Chabrier IMF we are using grows rapidly to a broad peak around 0.04 M_⊙, so for masses smaller than that value there is a "deficiency" of even smaller masses to form small mass ratio binaries relative to masses larger than 0.04 M_⊙. This is reflected in the mass ratio distribution in Figure 20 (seen by cutting vertically across the figure), which is much steeper between 0.012 and 0.04 M_⊙ than between 0.04 and 0.1 M_⊙.

Zoom In Zoom Out Reset image size

Another very important observational quantity is the orbital eccentricity. In contrast to the two previously discussed parameters, i.e., separation and flux ratio, determination of the eccentricity requires painstaking observations through radial velocity measurements or high resolution imaging (e.g., Joergens et al. 2010; Dupuy & Liu 2011), and it is consequently known much less frequently.

The distribution of eccentricities for 15,894 ejected BDs at 100 Myr is shown in Figure 21(a). It forms a monotonically increasing function toward higher eccentricities, and the median value is 0.74. This is not consistent with the limited information gathered so far on BD orbital eccentricities, which shows eccentricities concentrated toward values smaller than 0.6, with only few cases of higher eccentricity (Dupuy & Liu 2011; Biller et al. 2013). Even taking into account observational biases, the simulations and observations are not consistent. However, it must be recalled that while our simulations combine orbital dynamics and mass accretion, we do not include effects of gas dynamics. These are schematically indicated by red arrows in Figure 21(a), and are further discussed in Section 4.8. The high-eccentricity binaries will decrease in number as a result of viscous evolution, leading to an increase in low-eccentricity systems. Consequently, the initially rising distribution due to orbital dynamics is transformed by gas dynamics to a decreasing function. Additionally, circularization during the pre-main sequence phase creates circular orbits for the very shortest period binaries with periods less than ∼8 days (Zahn & Bouchet 1989).

Zoom In Zoom Out Reset image size

In an attempt to understand this distribution, we examine whether the eccentricity has any dependence on the three parameters mass ratio M_B/M_A, total system mass M_A+M_B, and semimajor axis a. Figure 21(b)–(d) plot these parameters against eccentricity. As is evident from the figures, there is hardly any dependence on these parameters. The only possible exception is that BD binaries with very small mass ratios have a somewhat weaker dependence on the eccentricity than BD binaries that are twins or nearly identical, which shows up as a barely visible gentle rise in a Bézier fit to the data points (red line in Figure 21(b)).

The BD binaries that are formed through expulsion were all once part of triple systems, and it is of interest to ask what are the properties of those third bodies and whether the nature of the third body correlates with any properties of the ejected BD binaries. Figure 22 plots the BD binary semimajor axis as a function of the mass of the single star to which the binary was previously bound. The main conclusion from the figure is that ejected BD binaries are most likely to have been bound to another BD. In other words, BD binaries mostly arise from BD triple systems that disintegrated. This partly reflects the shape of the IMF plus the assumption that all bodies in the initial triple systems have the same mass. The fact that most BD binaries were bound to another BD indicates that a large number of triple systems break up before accretion can significantly alter the masses of the components. It also reflects the fact that in unstable triple systems, it is predominantly the lowest-mass body that is ejected, and so if the binary consists of two BDs, then the third body is likely to have an even lower mass.

Zoom In Zoom Out Reset image size

It is also clear from the figure that binaries ejected from stars (as opposed to BDs) tend to have slightly smaller semimajor axes (mean around 20–50 AU), while the mean value for BD binaries that were bound to another BD is ∼50–60 AU.

Figure 23 shows a plot of the mass ratio of ejected BD binaries as a function of the mass of the third body. Spectral types are indicated for different masses, based on main sequence field stars. As already known from the triple diagnostic diagram, triple interactions with accretion are unlikely to produce S-low binaries, that is, BD binaries with small mass ratios do not originate from triple systems that include a star.

Zoom In Zoom Out Reset image size

As mentioned in Section 3, the binary components are denoted A and B, with A being the more massive, while C is the third body, either bound or escaped. For a three-body system consisting of stars and BDs, there are six possible combinations: either the third component C is a star or a BD, and for each of these two cases, A and B can be both BDs or both stars, or A can be a star with B as a BD. These six categories are listed in Table 1 for stable hierarchical systems (H), unstable hierarchical systems (U), and disrupted systems (D), all at ages of 1, 10, and 100 Myr. BD binaries that originate from disintegrated triple systems may either have had a BD or a star as a third member. At 1 Myr this adds up to 4.26% + 0.35% = 4.61% of the initial population of 200,000 triple systems in the simulations. At 10 Myr, this has grown to 7.40%, and at 100 Myr it is leveling out at 7.94%. As discussed in Section 4.6, the majority of these BD binaries originate from systems where the third body was also a BD.

An important observational parameter is the fraction of BDs that are binaries. This we can derive from Table 1, and the numbers at an age of 1 Myr are given in Table 2. The percentage of the original triple systems which have ejected a single BD at 1 Myr is 6.07% and, as mentioned above, the number of ejected BD binaries is 4.61%. In total this means that 10.68% of the 200,000 systems are resulting in either a single BD or a BD binary. Of these, the fraction of binaries are 4.61/10.68 = 0.43. This is higher than observed. However, the numbers assume that every BD binary is resolved and identified as such. Given that many of these binaries are very close and will be resolvable only through adaptive optics observations or through a radial velocity study, we must expect that a fraction of the binaries are not resolved, but appear to the observer as single BDs. Table 2 includes the BD binary fraction assuming that the unresolved binaries range from 25 to 50 to 75%. The actual binary fraction of 0.43 then declines to 0.32, 0.22 and 0.11, respectively.

Similar calculations at 10 and at 100 Myr yield the same binary fractions as at 1 Myr to within the two-digit accuracy we use here. The binary fraction is thus not age dependent. This is because the six categories of triple systems listed in Table 1 decay at essentially the same rate with time.

Note that in these calculations, we are not including the BD triple systems which remain bound. If the close BD binary is not resolved, such triple systems will appear as two widely separated BDs (typically ≫200 AU).

At first glance, this high fraction of BD binaries does not appear consistent with observations. In the most recent and most complete study, Todorov et al. (2014) found a binary fraction of 0.04$_{-0.01}^{+0.03}$ among young BDs for resolved binaries with separations >10 AU. With this resolution, the observations only probe longwards of the 11 AU peak of the BD binary separation distribution function (Figures 15 and 28), and so ignores the population of spectroscopic binaries, which are discussed in the following (Section 4.8). Although detections of spectroscopic BD binaries are still limited, the existing preliminary surveys indicate that such spectroscopic BD binaries should be rather common, Maxted & Jeffries (2005) suggest 17–30%, Basri & Reiners (2006) 26% ± 10%, and Allen (2007) about 20%. The total number of observed BD binaries is still less than the 0.43 of the simulations. But the calculated fraction assumes all BD binaries are ejected from triple systems. Ongoing simulations (in preparation) of higher-order multiples show that many more single BDs are ejected as soon as the number of bodies in the multiple system increases, and hence the BD binary fraction falls significantly. Finally, our simulations do not include the effect of viscous interactions, but if included a number of our binaries would shrink into the spectroscopic binary range. A more detailed comparison of binary statistics between observations and simulations must therefore await the study of higher-order multiples and a more complete survey for spectroscopic BD binaries.

Spectroscopic binaries are known to exist among the BD and VLM population (e.g., Basri & Martin 1999; Joergens 2008; Burgasser et al. 2012; Clark et al. 2012). Very few BD binaries are formed in our simulations with semimajor axes of only a few AU, which is the realm of the spectroscopic binaries, so evidently purely dynamical processes do not readily form such close pairs. The additional physical process missing in our simulations is the role of viscosity. The motion of a binary through the gas in a dense cloud core creates drag, and leads to an orbital decay of the binary components (Gorti & Bhatt 1996; Stahler 2010; Korntreff et al. 2012). In some cases even mergers can take place (Korntreff et al. 2012; Leigh & Geller 2012; Rawiraswattana et al. 2012). The majority of ejections and the resulting formation of a bound binary occurs during the embedded protostellar phase (Reipurth 2000; Reipurth et al. 2010) and at these early stages the stars are also surrounded by significant circumstellar disks, which interact during the periastron passages of the usually very eccentric orbits, resulting in further orbital decay (McDonald & Clarke 1995; Hall et al. 1996; Bate et al. 2002b). These processes stop to act when first the cloud core vanishes, and later the disks disperse, leaving the BD binary with whatever orbital parameters that resulted from the spiral-in process. Full smoothed particle hydrodynamics (SPH) simulations are required to document the role of these processes in forming the closest BD binaries. In consequence, the distribution of semimajor axes in Figure 15 reflects only orbital dynamics, and does not include the effects of gas-induced orbital decay. Similarly, the eccentricity of the viscously interacting BD binaries will gradually diminish, as schematically indicated in Figure 21(a).

From a practical point of view it is unfortunate that mergers are likely to occur for some of the close binaries in triple systems, since this implies that it is not possible to test the prediction made in Section 4.3 that BD binaries wider than ∼300 AU should be triple systems. If the close binary has suffered a merger, then the wide BD binary will in fact only consist of two objects, notwithstanding their origin in a triple system.

Figure 24 shows the binding energy as a function of the total binary mass for all ejected BD binaries. The general trend is that the binding energy increases with increasing mass. This is hardly surprising since the binding energy ($\propto {{M}_{1}}{{M}_{2}}/{{a}_{12}}$) is directly proportional to the individual component masses. The width of the distribution is partly due to the fact that BD binaries have mass ratios from 1 down to 0.3 (see Figure 18) plus variations in the semimajor axes. Figure 25 shows the distribution of semimajor axes as function of binary mass. As is well known from observations, and reproduced by the simulations, mean binary semimajor axes increase with binary mass. The effect is, however, much too small to counter the growth of the binding energy due to increase in binary mass. The overall behavior of the binding energy fits well with the trends seen in observations of VLM binaries (e.g., Close et al. 2007; Faherty et al. 2011). As Figure 24 shows, these BD binaries have exceedingly small binding energies, as low as 10⁴⁰ ergs, which is comparable to the values for the widest known binaries. While such wide binaries are vulnerable to breakup by passing stars (e.g., Weinberg et al. 1987), the compact BD binaries are stable against disruption, except possibly those that are born in dense clusters (Close et al. 2007).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The binding energy is not a measure of the likelihood that a BD binary could have survived an ejection event. When a triple system breaks up, a BD binary becomes permanently bound only at the moment when the third body is ejected. The binding energy of the BD binary is set at the moment of ejection and it can become high or low, depending on the specific circumstances of ejection. It is thus not related to its internal survivability during system breakup. If there are more bodies present, e.g., in a higher-order system or a cluster, then a weakly bound binary is vulnerable to external perturbations, as mentioned above.

When a triple system breaks up and a single star escapes, then the binary recoils in the opposite direction. It is of interest to ask if the recoil leaves a measurable imprint on the binary space motion that could be used to identify this violent event in its pre-history. However, numerical studies of BDs ejected from disintegrating triple systems show that while they may briefly during the actual ejection process achieve a high velocity, they very soon slow down as they climb out of the steep potential well of the cloud core and the binary (Bate et al. 2003; Bate 2009) The mean terminal velocity of all escapers in the 12,800 simulations by Reipurth et al. (2010) is $1.13\;{\rm km}\;{{{\rm s}}^{-1}}$. Since the recoil velocity of the binary is often lower because the binary tends to have a larger total mass than the single, it follows that BD binaries formed through disintegration have a velocity dispersion that is well within the turbulent velocity range of star forming clouds, and they are thus indistinguishable from BD binaries formed by other mechanisms.

Reipurth & Clarke (2001) noted that the dynamical interactions involved in the ejection scenario implies that circumstellar material around BDs will be truncated, and they suggested that ejected BDs consequently should have less massive disks with weaker accretion signatures and less or shorter-lived far-infrared excess emission. Somehow, this has often been construed in the literature to mean that young BDs that are ejected should not have disk-signatures like infrared excesses. That is obviously not what was stated by Reipurth & Clarke (2001), nor what observations indicate (e.g., Monin et al. 2010; Bulger et al. 2014). While the outer regions of disks evidently are truncated in a close triple approach (e.g., Hall et al. 1996), this is a temporary situation since the viscosity of the disk material ensures that the disks spread out again, allowing a full range of temperatures, and accretion will also continue as long as disk material remains.

Umbreit et al. (2011) have studied the effect of triple interactions on the disks around BDs and the observable consequences. They find that disks after close triple approaches are mostly less massive, but compared to disks after two-body encounters they have similar or larger radii. In a novel result, they also show that part of the disk material initially lost is in fact re-captured by the ejected body (see also, e.g., Steinhausen et al. 2012). Hydrodynamic simulations like those of Bate et al. (2003), Bate (2009) are evidently important to study such truncations in detail, although as pointed out by Bate (2012), one should be aware that interactions in a large-N cluster may well be different from the low-N systems we are considering here.

The same arguments that apply to single BDs are relevant for binary BDs as well, in fact one may argue that the presence of a companion provides another opportunity to capture otherwise escaping disk material. It therefore seems that young BD binaries, just like young single BDs, are likely to harbor circumstellar and/or circumbinary disks, although loss of disk material might shorten the associated signatures of youth.

The very same processes described here for isolated triple systems are almost certainly also taking place in clusters. The main difference in the dynamical evolution of triple systems in loose associations and in dense clusters is that the distant weakly bound components which in low-density star forming regions in time will disrupt due to internal instability, in clusters may be dislodged first due to external perturbations. Parker & Goodwin (2011) have shown that VLM binaries with separations less than ∼20 AU are mostly immune to disruption in even the densest clusters, while most VLM binaries with separations of more than ∼100 AU can be destroyed in high-density clusters, but are generally unaffected in low-density clusters. Since most stars and BDs are formed in clusters, the effect is that the peak and median values of the semimajor distribution of BD binaries (Figure 15) will be lowered, in addition to the effect expected from viscous interactions. This does not mean that wide and very wide BD binaries necessarily must be rare, because triple systems that eject a third body (i.e., dynamically "unfold") after the initial clumpy substructure of clusters has been smoothed out and the cluster starts to expand will have been largely protected against disruption (Reipurth & Mikkola 2012).

The main parameters that control the survival of a wide binary is thus the unfolding time itself (which depends primarily on the semimajor axis of the system), the cluster mass and relaxation timescale, and the length of time since formation of the cluster until the triple system begins transformation into a hierarchical configuration. Each cluster will have its own signature on the binary/triple separation distribution function, which will be truncated for separations exceeding a cluster-dependent value (e.g., Kroupa 1995; Reipurth et al. 2007).

Comparison of simulations with observations is hampered by the fact that simulations frequently deal with semimajor axes a while observations commonly only measure projected separations s. Additionally, the orbital eccentricity plays a role. As demonstrated by van Albada (1968), <log s> − <log a> ranges from −0.13 for circular orbits to as little as −0.03 for highly eccentric orbits (see also Kuiper 1935; Couteau 1960; Halbwachs 1983). When interpreting observations of individual binaries, it is sometimes forgotten that for highly eccentric orbits the component separation can be almost twice the semimajor axis, and in such systems the fraction of time that a companion spends at separations larger than the semimajor axis can reach 82% (see Figures 26 and 27). This is particularly important to consider when dealing with systems originating from triple systems, which often decay to highly eccentric orbits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following we calculate the projected separations of the BD binaries resulting from our simulations at different ages, allowing a more direct comparison with observations.

Since we know the space coordinates of all our BD binaries at all times, it is trivial to derive their projected separations. Figure 28 shows projected separations of the 9209 BD binaries that have already been ejected at 1 Myr. Since none of these binaries are suffering any perturbations following their ejection, the distribution does not change with time, it is only the number of binaries that increases as more unstable triples break up. The distribution has a broad maximum between 2 and 20 AU, in which 35% of all ejected BD binaries can be found. The peak is around 11 AU, and the median is at 28 AU.

The projected distribution differs from the distribution of semimajor axes (Figure 15) in two respects. First, the random distribution in space of the semimajor axes washes out the deep dip at very small semimajor axes, which results because the simulations do not include viscous interactions due to motion inside the gas core and to disk–disk interactions at periastron. The peak at 11 AU is therefore an upper limit to the peak of BD binary separations. Observationally, it is well established that binary separations decline with diminishing stellar mass. Fischer & Marcy (1992) found that M-star binary separations peak around 4–30 AU, and Close et al. (2003) and Maxted & Jeffries (2005) both find a rather sharp peak of VLM/BD objects around 4 AU, while Burgasser et al. (2007) estimated the peak around ∼3–10 AU. This is, as expected, more compact than our peak of projected separations around 11 AU due to the lack of viscous effects in our simulations. Second, the projected distribution has a broader tail at larger separations, due to the fact that for systems with high eccentricity, observers are likely to see a system near apastron rather than periastron, as discussed above (Section 5.1).

The three panels in Figure 29 show the projected separations of the outer body in triple systems relative to the center of mass of the inner binary for systems in which all three components are BDs at the ages of 1, 10, and 100 Myr, respectively. The systems are divided into bound systems, which are either stable (red) or unstable (green), or unbound systems in which the third body has already been released from the system (blue). The black dotted line indicates the sum of all three types of systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At 1 Myr, BD triples with projected separations less than approximately 10,000 AU are almost all still bound, but mostly unstable. At larger separations, the BD triple population is dominated by disintegrated systems, where the third body is gently traveling away from the binary to ever-increasing separations. Stable triple systems consisting of three BDs are very rare, of the 15,376 BD triples, only 117 systems are classified as having stable orbits at 1 Myr. The majority of these are compact, falling in the bin with projected separations up to 2000 AU, and only in this bin are the stable bound systems more common than the already disintegrated systems. The vast majority of BD triples are either bound but unstable or already disrupted; out to about ∼10,000 AU the bound unstable systems dominate the disrupted systems, but at larger separations this reverses (see Figure 30). At 1 Myr the number of BD triple systems that are classified as stable, unstable, and disrupted are 117, 6748, and 8511, corresponding to 0.7, 43.9, and 55.4%. Evidently, the breakup of triple systems is very common at ages less than 1 Myr, as discussed by Reipurth et al. (2010).

Zoom In Zoom Out Reset image size

At 10 Myr, such as several populations in moving groups, significant further dynamical evolution has taken place. At this age, the number of BD triple systems that are classified as stable, unstable, and disrupted are 111, 1365, and 13,900, corresponding to 0.7, 8.9, and 90.4%. While almost all systems classified as stable at 1 Myr remains so at 10 Myr, 80% of the population of unstable triples at 1 Myr has disrupted by 10 Myr. At the same time, the population of already disrupted systems have moved even further apart, so in terms of projected separations, most systems out to ∼45,000 AU are now bound but unstable, and only at larger separations do the disrupted systems become more common (see Figure 30).

At 100 Myr, when almost all young groups have blended into the field population, the number of BD triple systems that are classified as stable, unstable, and disrupted are 110, 347, and 14,919, corresponding to 0.7, 2.3, and 97.0%. At even larger ages, it is likely that the last unstable systems will have disintegrated. In other words, in the field population only less than 1% of BD triple systems should still be bound.

It should be kept in mind that in many BD triple systems, viscous interactions in the close binary will cause a spiral-in of the components so they form a spectroscopic binary, and in some cases will even merge (see Section 4.8). Such systems will appear as wide BD binaries. It is therefore of interest to note that Close et al. (2007) have estimated that ∼6% ± 3% of young VLM objects are found in wide (>100 AU) systems, while only 0.3% ± 0.1% of old field VLM objects are found in such wide systems. This finding is supported by the more recent observations of Todorov et al. (2014).

It is well known that BD binaries in some cases are found as companions to main sequence or evolved stars. A fine example is the G2V star HD 130948, which has a L4+L4 pair at a projected separation of 47 AU from the star and orbiting the star with a very high eccentricity, with the eccentricity probability function peaking at 0.83 (Potter et al. 2002; Dupuy et al. 2009; Ginski et al. 2013), see Figure 31. The BD binary itself has a projected separation of ∼2.4 AU and an orbital period of about 10 yr (Dupuy & Liu 2011).

Zoom In Zoom Out Reset image size

We here argue that such BD binary companions are a straightforward consequence of dynamical interactions in triple systems. Basically, such stars with BD binary companions are the triple systems where the BD binary was not successfully ejected into an escape, but was ejected into a stable orbit around the main component. Our simulations show that such systems can originate from three identical stellar embryos, where one by chance has controlled a position near the center of the molecular core and rapidly grew to stellar mass while dynamically keeping the other two components in the outskirts of the core, where they were unable to gain much mass.

Such systems are rare, and become even rarer as the systems evolve, since the majority are unstable and therefore much dynamical evolution takes place at early stages. At 1 Myr, there are 1325 BD binaries bound to a star among our 200,000 simulations, corresponding to 0.66%. By 10 Myr less than half remain (582, corresponding to 0.29%). At 100 Myr only 377 are left (0.19%), and of these 253 (0.13%) are stable and are likely to remain so on much longer timescales.

The number of BD binaries bound to a star decreases with time, while the number of ejected BD binaries increases with time. Hence the ratio of BD binaries bound to a star relative to the number of free-floating already ejected BD binaries varies with time, at 1 Myr it is 1325/9209 (∼1:7), at 10 Myr 582/14,814 (∼1:25), and at 100 Myr 377/15,894 (∼1:42).

The large majority of BD binaries bound in a triple system are bound to another BD, not to a star. So if we ask the more general question of what fraction of all BD binaries that are bound (to either a star or a BD) relative to all BD binaries (bound plus unbound) we get 47.1% at 1 Myr, 12.2% at 10 Myr, and 5.0% at 100 Myr.

The semimajor axis distribution of the BD binaries bound to stars is seen in Figure 32. The peak of the distribution is around 25 AU, and the median value is about 59 AU. The large majority of these bound BD binaries are orbiting the central star on wide orbits, with semimajor axes of typically thousands of AU (Figure 33).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Only a small number of BD binaries bound to stars have been discovered so far, and with such small-number statistics, it is difficult to make a meaningful comparison between observations and simulations. But the existence of a BD binary associated with $\epsilon$ Indi (e.g., King et al. 2010), which is located at a distance of 3.6 pc and is one of only about 25 stars this close, may suggest that such BD binary companions are not extremely rare. Comparison with observations is complicated by the fact that many of these bound BD binaries require adaptive optics imaging to be resolved, which in many cases limit discoveries to binaries with projected separations to their host star of only a few arcseconds. As widefield adaptive optics systems become more common, more wide systems are expected to be discovered, enabling meaningful comparisons between simulations and observations. From the limited observations available, Burgasser et al. (2005) and Faherty et al. (2010, 2011) suggest that stars are more than twice as likely to have a BD binary rather than a single BD as a companion, which strongly points to a dynamical origin of BD binary companions to stars.