FROM BINARIES TO MULTIPLES. II. HIERARCHICAL MULTIPLICITY OF F AND G DWARFS

A system of three or more bodies with comparable separations is dynamically unstable (Harrington 1972). In contrast, hierarchical multiples survive for a long time because the motion in the close (inner) pair is not strongly perturbed by the outer companion(s). Stellar motions in stable hierarchies are represented approximately by Keplerian orbits. This is the definition of a hierarchical multiple system. This definition is not very sharp because some systems approach the stability boundary. Moreover, the orbital motion of very wide binaries is not observable owing to their long periods, so the borderline between the truly bound multiple systems and other stellar groups like disintegrating clusters is fuzzy.

Hierarchies can be described by binary graphs or trees, also called mobile diagrams (one example is shown in Figure 1). The position of each pair in the tree is called level and is denoted here as L1, L11, etc. The outermost (widest) pair is at the root of the tree, L1 (the root is denoted by an asterisk in the parent designation). Inner pairs associated with the primary and secondary components are at levels L11 and L12, respectively, and this notation continues to deeper levels. Triple systems can have inner pairs either at L11 or at L12. When both primary and secondary subsystems are present, we get the so-called 2+2 system.² Alternatively, a 3+1 quadruple system consists of levels L1, L11, and L111; it has three companions associated with the same primary star, resembling in this respect a planetary system. Both types of hierarchy are found in HIP 88637. The assignment of "primary" and "secondary" is traditionally based on the apparent brightness of the components. When they are comparable in brightness or mass, the distinction between primary and secondary becomes irrelevant.

Zoom In Zoom Out Reset image size

While we deal with binary and triple systems, their hierarchy is uniquely defined by the number of components, so levels are almost redundant. On the other hand, quadruple and higher-order systems can have different organizations, so their hierarchy requires additional descriptors such as levels. From the observational point of view, it is easier to explore lower hierarchical levels L1, L11, and L12 than to determine the number of multiple systems with exactly 3, 4, etc., components.

The fraction of non-single stellar systems, the multiplicity fraction f_M, equals the fraction of L1 systems in a given sample. The companion fraction f_C is the number of systems at all levels divided by the sample size, so f_C > f_M. It is easier to measure f_M, as it does not depend on the discovery of all subsystems. Multiplicity is also characterized by the fractions of systems with exactly n components, f_n. Obviously, the fraction of single stars f₁ = 1 − f_M, while f_C = f₂ + 2f₃ + 3f₄ + ... The fraction of hierarchies in a sample is $f_H = \sum _{n=3}^{\infty } f_n = 1 - f_1 - f_2$. The multiplicity statistics are defined by f_n, by the distribution of levels, and by the joint distributions of periods and mass ratios at different levels.

The system depicted in Figure 1 contains hierarchical levels most frequently found in the FG-67 sample: L11 (N = 296), L12 (N = 95), and L111 (N = 17). Here the statistical analysis of those three levels is made, because higher levels are rare, and a much larger sample would be needed for their study.

In a binary system with the primary component A and the secondary component B, the mass ratio is defined as q = M_B/M_A < 1. In a triple system where B is a close pair B,C, the mass ratio can be defined in different ways (e.g., Tokovinin 2008). The classical approach is to use the mass sum of subsystems in the system mass ratio q_sys = (M_B + M_C)/M_A. However, this definition may result in q_sys > 1 even when the component A is the most massive and brightest star in the system. Alternatively, the component mass ratio q = M_B/M_A takes into account only the masses of individual stars (primary components of respective subsystems). For example, in a triple star consisting of the close binary A,B with a distant tertiary component C, q = M_C/M_A, while q_sys = M_C/(M_A + M_B). The detection limits of imaging techniques are based on the magnitude difference, which is more tightly related to the component mass ratio q than to q_sys. The use of q avoids re-definition of "primary" and "secondary," which would be needed in many instances if we used the mass sum instead. In the following, we use consistently the component mass ratio q for multiple systems.

The data on all systems and subsystems in the FG-67 sample are assembled in the structure (table), which contains information on the hierarchical level (with links to parents and siblings), period, and component's masses. The same structure is produced for the simulated samples. As binary periods P span a huge range, we use instead the logarithm x = log₁₀(P/1day) throughout this paper.

The 22 targets containing white dwarfs (WDs) are removed, leaving 2162 systems and subsystems to consider (the sample excludes evolved primary components by design; WDs are such former primaries). The fraction of targets with WD companions is from 2% to 4%, according to the rather uncertain estimate of Paper I. As only 0.5% of known binaries are removed and most of the remaining WDs are undiscovered, the total sample size without WDs should be somewhere between 4650 and 4750. The sample size affects the estimates of f_M, f_C, etc. Here we set it conservatively to 4800 (i.e., remove 1% of targets with WDs). The multiplicity fraction derived here would be larger by ∼1% if more targets with WDs were removed. Interestingly, among the 22 binaries with WDs, 11 are at least triple, that is, their FG-type primaries (former secondaries) contain subsystems. This supports the large fraction of secondary subsystems found here in the main sample.

The statistical data are presented in Table 1. Its columns contain the HIP0, the Hipparcos number of the primary component, hierarchical level L, type of the system (see Paper I), logarithm of the period x (in days), system masses of the primary and secondary components M₁ and M₂ (in solar mass), and the component masses M_{1, p} and M_{2, p} (primaries in the inner subsystems). The two definitions of mass ratio in Section 2.3 correspond to q_sys = M₂/M₁ and q = M_{2, p}/M_{1, p}. Masses are taken from the SYS table of Paper I and are estimated from absolute magnitudes using the standard relation for main-sequence stars. Unknown masses have zero values; unknown periods are assigned x = −9. The links to children and parent (record numbers in the data structure, or −1 if absent) are generated automatically using component designations; they are not printed in Table 1.

A significant number, 355 out of 2162 systems, have unknown periods and separations. These are spectroscopic binaries without known orbit and binaries discovered by Hipparcos accelerations (see Paper I). Masses of their secondary components are also unknown. The statistical analysis presented below deals with this missing information (32% of all binaries with periods shorter than 100 yr). The detection probability of these binaries is properly modeled, so we cannot simply ignore them. The assumption is made that all unknown periods are shorter than ∼100 yr. Random fictitious periods are generated to replace missing data when plotting period distributions, but these fictitious data are not used in any other way. Considering the detection limits of spectroscopy and astrometry, fictitious periods of spectroscopic binaries are uniformly distributed between 0.5 < x < 3.2, and acceleration binaries with constant radial velocity have 3.2 < x < 4.5, e.g., HIP 290 and HIP 359 in Table 1. Partially missing information is the weakness of the FG-67 sample, which affects the statistics at periods shorter than 100 yr.

Figure 2 shows the position of known companions to primary targets and to the secondaries in the (P, q) plane. Binaries with unknown parameters are plotted by small symbols to indicate their existence and likely location in this diagram. The contours show average detection probability. The census of low-mass companions with separations from 03 to 10'' is still rather incomplete; other regions of the parameter space are covered better. The mass ratios are distributed almost uniformly at all periods. The period distribution has a broad maximum and declines on both sides of it. Note the cluster of twin binaries with q ≈ 1 at P < 100 days (Tokovinin 2000; Lucy 2006). Figure 3 shows cumulative distributions of periods at three hierarchical levels (with fictitious periods included).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since Duquennoy & Mayor (1991) fitted a lognormal function to the distribution of periods of 82 binaries, following the earlier work by Kuiper (1935), this analytical model has been universally adopted. The distribution of the mass ratio is often approximated by the power law f(q)∝q^β (Duchêne & Kraus 2013). It is still being debated whether the exponent β depends on period and mass or is universal (Reggiani & Meyer 2013). Following the tide, we represent the period distribution by the truncated lognormal law and the q distribution by the power law with β > −1, truncated at q < 0.05 (ignore brown-dwarf companions). The analytical model is

$$\begin{equation} f(x,q) = C \; \epsilon \; q^\beta \exp [ - (x - x_0)^2/(2 \sigma ^2)], \end{equation} \tag{ 1 }$$

where x = log₁₀(P/1day). The period distribution is truncated at x < −0.3 and x > 10 and re-normalized accordingly by the constant C; without truncation, $C = (1 + \beta)/(\sigma \sqrt{ 2 \pi })$. The fitted parameters are multiplicity fraction , median period x₀, period dispersion σ, and the mass-ratio exponent β. The assumption that β is the same at all periods is tested in the following subsection. The correspondence between and f_M or f_C depends on the choice of binaries; see Section 3.3.

The parameters of the model (1) were fitted to various subsamples using the maximum likelihood (ML) method. It takes into account the detection limits and deals with the missing data. The algorithm and its implementation are described in the Appendix.

Potential dependence of mass ratio on binary period can be studied here, taking advantage of the large sample size and known detection completeness. Binaries in certain period intervals were selected and fitted with just two parameters (, β) by the ML method (the Appendix). The detection probability d(x, q) was averaged over these period intervals, leaving only its dependence on q. We did not fit the histograms of q, avoiding any binning, but they match the derived β reasonably well (Figure 4).

Zoom In Zoom Out Reset image size

The results of such calculations in the overlapping period intervals of one decade are plotted in Figure 5. Each interval contains from 70 to 300 binaries. Only companions to the main targets are used here, regardless of their hierarchical level. The secondary subsystems have different, smaller completeness, and are too few (95 total). At both short and long periods, the distribution of the mass ratio is uniform, with β ≈ 0. This matches the results of other studies (Tokovinin 2011; Reggiani & Meyer 2013).

Zoom In Zoom Out Reset image size

In the period interval x = [3, 5], the derived β > 0.5 implies preference of large mass ratios. Indeed, most known binaries with such periods have q > 0.5 because companions with q < 0.5 have a low detection probability (see Figure 2). The detection probability is taken into account in the calculation of β; however, it ignores binary systems with unknown P, q, while including their discovery by acceleration and radial velocity in the evaluated completeness. To access the importance of this effect, we assigned fictitious periods and mass ratios to the binaries with missing data and repeated the calculation. The fictitious mass ratios are uniformly distributed in the interval [0.2–0.8] at periods shorter than x = 3.5 and in the interval [0.2, 0.6] at periods 3.5 < x < 4.5 (if such binaries had larger mass ratios, they would have been resolved). By complementing the missing data in this reasonable but arbitrary way, we obtain smaller values of β in the problematic period interval 3 < x < 5 (dashed line and crosses in Figure 2). We conclude that the large values of β likely result from the missing data and are not real.

There remains a hint of a mild dependence of β on period even at x > 5, where there are no missing data. However, the range of this variation, if any, is small. We fit the parameters of the model (1) by assuming a constant β. This helps to minimize the effect of missing information at intermediate periods. The global fitting in Section 3.3 confirms the small values of β.

Table 2 lists the parameters of the distributions found by the ML method, for all pairs without distinction of their hierarchy and for the subsamples. Its second column gives the number of binaries selected. The errors corresponding to 68% confidence intervals ("1σ") are evaluated by the ML method. We checked that there is no strong correlation between parameters. However, one should keep in mind that these errors do not account for any problems in the input data or their modeling. The formal ML errors give a lower bound of the true errors.

Table 2. Fitted Distribution Parameters

Case	N		x0	σ	β
All pairs	2162	0.571	4.54	2.40	0.094
		±0.012	±0.06	±0.06	±0.020
L1	1747	0.464	4.93	2.34	0.051
		±0.011	±0.06	±0.06	±0.011
L11	296	0.214	3.25	1.80	0.121
		±0.012	±0.12	±0.09	±0.026
L12	95	0.157	2.67	1.68	1.32
		±0.016	±0.17	±0.10	±0.28

Download table as:  ASCIITypeset image

When all pairs are selected, their fraction = 0.57 equals the companion fraction f_C. However, the code does not consider the reduced probability of detecting subsystems in the secondary components, hence delivering the slightly underestimated f_C; f_C  =  0.65 follows from the simulations presented below in Section 4.3. When only the pairs at L1 (outer systems) are selected, we obtain the multiplicity fraction = f_M = 0.46 ± 0.01. It agrees well with f_M = 0.46 found by Raghavan et al. (2010) for solar-type stars.

For the subsystems at levels L11 and L12, the presence of the outer pair at L1 is a necessary condition. Therefore, the resulting frequency refers to the sample of L1 systems, not to the whole sample. The ML results indicate that a fraction 0.464 × 0.214 = 0.100 of all targets have at least two companions around their primary stars. Similarly, a fraction 0.464 × 0.157 = 0.073 of targets have binary secondaries. This last estimate depends on the larger (and less certain) correction for the low probability of detecting secondary subsystems.

The total fraction of hierarchical systems f_H is less than the sum of the above numbers (17%) because those two groups overlap (see below). A fraction of about 0.6 of L12 (secondary) subsystems also have an L11 (primary) subsystem. The fraction of hierarchical systems is therefore estimated by ML as f_H = 0.100 + 0.4 × 0.073 = 0.129. A ballpark estimate of the uncertainty of this number is ±1%. Simulations presented in Section 4.3 confirm the derived fraction of hierarchies. Our initial analysis (before removing the WDs) gave f_H = 0.14.

Figure 6 shows the period histograms of all pairs and of primary subsystems at L11. Randomly assigned fictitious periods do not allow a rigorous comparison with the model in the x < 4.5 zone. The incompleteness in each period bin is corrected by averaging d(x, q) over q, which is valid for the uniform mass-ratio distribution only. Despite these shortcomings, the lognormal period distribution appears to be a good match to the data (the dotted curves were not fitted to the histograms!).

Zoom In Zoom Out Reset image size

Previously, when the samples of binaries were small, the lognormal period distribution was an adequate model. However, its wings extend to unrealistically short and long periods; for this reason we use the truncated Gaussian here. A better analytic model of the period distribution could be a triangular function with a rounded top,

$$\begin{equation} f_T(x) = \frac{1}{Cw} \left[ 1 - \sqrt{ a^2 + ((x-x_0)/w)^2 (1 - a^2)} \right], \end{equation} \tag{ 2 }$$

where |x − x₀| < w. This distribution extends over a finite range from x₀ − w to x₀ + w; the parameter w is close to the half-width. The parameter a provides for the rounded top (a = 0 means pure triangle 1 − |x/w|), and the normalization factor C depends on a: C = 1 − a²(1 − a²)^−1/2ln [(1 + (1 − a²)^1/2)/a]. The need for three parameters x₀, w, a instead of just two x₀, σ may be a disadvantage of this model, but on the other hand the triangular distribution does not require truncation.

The ML formalism was adapted to use the triangular function instead of a Gaussian. The results obtained with this alternative model are quite similar; they are not given here, for brevity. The parameter a ≈ 0.25 was found, meaning that the top of the triangle is rounded only mildly. The "shoulders" of the period histograms in Figure 6 are indeed closer to straight lines than to a Gaussian.

Usually the distributions of binary parameters were studied regardless of hierarchy, merging all pairs together (DM91; R10). However, there is a significant (factor 2.5) difference between the median period of all binaries and the median period of outer binaries at level L1.

The periods of L11 subsystems are shorter than the periods of all binaries, being restricted by the dynamical stability. The period histogram at L11 deviates from the Gaussian model, showing a depletion at x ∼ 2 and an excess at x < 1. Considering the uncertainty associated with the missing data, we cannot evaluate the significance of this mismatch. It is expected due to tidal migration of inner subsystems toward short periods (Section 4.4).

Note that the secondary subsystems at L12 are only slightly less frequent than at L11. The common belief that inner systems are preferentially found around primary components turns out to be true, but this preference is not strong. The median mass of the primary components at levels L11 and L12 is 1.20 M_☉ and 0.79 M_☉, respectively.

The subsystems at L12 show a tendency toward equal components, β = 1.3. Only 5 secondary subsystems out of 95 have unknown periods, making it difficult to explain high β by the influence of the missing data. Considering the bias toward large q in the detection of secondary subsystems (see Figure 2), the high β derived here for L12 should be verified by further work before being accepted. However, β ∼ 0 would mean that more secondary subsystems are currently missed, so their frequency would be even higher than the frequency of L11.

Figure 7 compares orbital periods P_S and P_L at the inner and outer hierarchical levels, respectively (only known periods are plotted). It shows that the lowest period ratio P_L/P_S is limited by the dynamical stability, but this ratio can also take large values; there is no typical or preferred period ratio. The points occupy almost all the space above the stability cutoff, showing that all allowed combinations of periods actually happen.

Zoom In Zoom Out Reset image size

One striking feature of Figure 7 is the absence of outer periods shorter than 10³ days, as noted already by Tokovinin et al. (2006). This is not caused by the observational selection, as tertiary companions with short periods are readily discovered by radial velocity variation; the detection probability at x < 3 is high (see Figure 9 of Paper I). Gies et al. (2012) independently discovered the lack of short outer periods by eclipse timing of 41 binaries from Kepler: 14 of them show trends indicative of tertiary companions, but none have tertiaries with P_L < 700 days. In the three-tier hierarchies L11+L111, the shortest period at L11 is x = 2.88, confirming the exclusion zone of outer periods found for the L1+L11 and L1+L12 hierarchies. That said, triple systems with short outer periods are known among more massive stars (e.g., λ Tau with P_L = 30 days).

The minimum period (or separation) ratio allowed by dynamical stability has been studied by several authors. For example, the stability criterion of Mardling & Aarseth (2001) for co-planar orbits can be written as

$$\begin{equation} P_L/P_S > 4.7 (1 - e_L)^{-1.8} (1 + e_L)^{0.6} (1 + q_{\rm out})^{0.1}, \end{equation} \tag{ 3 }$$

where e_L is the eccentricity of the outer orbit, while the ratio of the distant-companion mass to the combined mass of the inner binary q_out plays only a minor role. The solid line in Figure 7 corresponds to P_L/P_S = 4.7; all points are above it (with two exceptions caused by the uncertainty of period estimates from projected separations). Tokovinin (2004) studied the dependence of the period ratio on e_L for triple stars with reliable outer orbits and suggested an empirical stability limit with a stronger dependence on e_L,

$$\begin{equation} P_L/P_S > 5 (1 - e_L)^{-3}. \end{equation} \tag{ 4 }$$

At e_L = 0.67, the limiting period ratio (3) becomes a factor of 10 larger than at e_L = 0. The dynamical truncation of P_L/P_S is random and extends over at least one decade, depending on the distribution of e_L. Figure 8 shows the cumulative distributions of this cutoff Δx = log₁₀P_L/P_S if a cosine distribution of e_L between 0 and 0.8 is assumed, with average e_L = 0.4. Available data on multiple stars indicate that orbits of outer systems tend to have moderate eccentricities (Shatsky 2001); the linear distribution f(e_L) = 2e_L seems to be excluded.

Zoom In Zoom Out Reset image size

Here we use a crude model for the dynamical truncation function T(Δx) (probability of a triple system to be dynamically stable) by setting T = 0 for Δx < 0.7, T = 1 for Δx > 1.7, and linear in-between (dotted line in Figure 8). This model agrees with the stability criterion of Mardling & Aarseth for the assumed distribution of e_L, while the criterion (4) predicts a stronger cutoff. An attempt to determine the truncation empirically by the distance of points in Figure 7 above the solid line did not produce any convincing results. Most outer and some inner periods are determined from projected separations to within a factor of three, so the empirical values of Δx are approximate. For the same reason, precise knowledge of the truncation function T(Δx) is not needed here and its crude model is sufficient.

The distributions of both inner and outer periods in the allowed region of Figure 7 look similar, as though they were drawn from the same population. Are the periods at levels L11 and L12 shorter than at L1 simply because of the dynamical cutoff?

The idea that the inner and outer pairs of multiple stars are selected from the same (or similar) parent distributions, subject only to the stability constraint, merits further investigation. Several arguments in its favor are furnished by prior work. Tokovinin et al. (2006) noted that the mass ratios of close binaries with and without outer (tertiary) companions are distributed equally, although the former group has statistically shorter periods. The frequency of spectroscopic subsystems in visual binaries is similar to the frequency of spectroscopic binaries in the open clusters and field (Tokovinin & Smekhov 2002). The frequency and mass ratio of resolved subsystems in wide binaries are again comparable to the binaries of such separations in the field (Tokovinin et al. 2010). Finally, the angular momentum vectors in the inner and outer systems show only a weak correlation (Sterzik & Tokovinin 2002), and there are well-documented cases of counter-rotating multiples.

Independent multiplicity is expressed mathematically by setting the joint period distribution at levels L1 and L11 to be a product of the individual distributions at these levels, including the dynamical truncation factor T(Δx):

$$\begin{equation} f (x_{{\rm L11}}, x_{{\rm L1}}) = f_{{\rm L11}}(x_{{\rm L11}}) f_{{\rm L1}}(x_{{\rm L1}}) T(x_{{\rm L1}} - x_{{\rm L11}}). \end{equation} \tag{ 5 }$$

Figure 9 illustrates this statistical model. The generating distribution for L11, f_L11(x_L11), could be the same as the period distribution at L1, f_L1(x_L1). This logic can be applied to other levels, leading to a complete statistical description of hierarchical multiplicity. Comparison with the real data presented below shows that this model is too simplistic and needs several adjustments. However, to the first order, it works even in its simplest form.

Zoom In Zoom Out Reset image size

An attempt was made to test this idea using the ML method. The truncation function was added as a multiplier to the distribution (1), and the four parameters were found for the subsystems at L11 and L12. The resulting distributions can be thought of as the generating functions f_L11 and f_L12. The problem of this approach is that almost all subsystems have x < 6 (see Figure 7), so the data do not constrain the generating functions at longer periods. This means that the fitted parameters of the generating distributions , x₀, σ do not inspire confidence, and only the distributions at x < 6 (left halves of the Gaussians) are more or less defined. Moreover, a Gaussian distribution is not a good model for the subsystems at L11 (Figure 6). We do not list here the fitted parameters and simply plot the relevant parts of the generating distributions in Figure 10. The periods of L12 subsystems turn out to be shorter, and their frequency at x < 4 is even higher than the frequency of levels L1 and L11. However, considering the uncertainties involved in this crude model, we can only state that all three curves are not very different from each other and that selecting periods at outer and inner levels from the same distribution is not such a bad idea.

Zoom In Zoom Out Reset image size

In this subsection, we simulate the multiplicity of the FG-67 sample. A synthetic population of single, binary, and multiple stars is created by a random draw from the generating distribution (1) with parameters (x₀, σ, β) = (5.0, 2.3, 0). The period distribution is truncated at x < −0.3 and x > 10. The binary frequency is specified below. For each component of a binary pair with x > 3, the subsystems are drawn from the same distribution and retained with a probability T(Δx) defined above. The recipe is applied recursively to produce high-order hierarchies. After that, the synthetic sample is "filtered" by the average detection probability of the real sample (different for the primary and secondary components) and compared to the real sample.

Multiple systems with more than two hierarchical levels (three-tier hierarchy) pose a specific problem for their detection. Suppose that in a quadruple (((Aa, Ab), B), C) the component B is missed (un-detected). The system will be perceived as a triple ((Aa, Ab), C), i.e., we miss the intermediate hierarchical level; the pair Aa, Ab is attributed to L11, instead of L111. Now suppose that the true system is ((A, (Ba, Bb)), C) and we miss the same component B. In such a case, we can't discover (Ba, Bb), and the observed system becomes a simple binary (A, C). Our detection simulator eliminates all siblings of un-detected secondary components but keeps those of the primary and corrects the hierarchical levels to reflect the observed, rather than true, structure of each system.

Straightforward application of the above recipe with a fixed overproduces binaries and underproduces triples, while the total number of L1 systems (the multiplicity fraction) is correct. Therefore, subsystems at different hierarchical levels do have some correlation. The simulator is modified by adopting a variable binary frequency (however, is the same for each individual multiple system). When is high, many hierarchical multiples are produced, and when it is low, we get mostly single stars and simple binaries. This reflects the idea that the field sample is a mixture of "binary-rich" and "binary-poor" populations. The average is constrained by the observed multiplicity fraction f_M, while the relative proportion of hierarchies informs us on the variation of . The frequency of binary stars is proportional to , the frequency of triples to ², etc. By adopting a variable , we enhance the fraction of hierarchies until it matches the real sample (the same result could be obtained by increasing for inner hierarchical levels).

We split the simulated sample into three equal parts and apply to each group its own frequency _i. The correct fraction of hierarchies was obtained, for example, for _i = [0.05, 0.60, 0.75]. These numbers were found by trial and error; their mean 0.466 is the simulated multiplicity fraction. Compared to the constant multiplicity, the fraction of triples is increased by 〈²〉/〈〉² = 1.41. There are other ways to model variable by adopting some distribution and playing with its parameters. The distribution of chosen here is discrete with three equally probable values. Although the values of quoted above produce a good match to the data, other combinations or recipes are not excluded.

The simulator needs yet another adjustment to reproduce the large observed number of 2+2 systems. If primary and secondary subsystems at levels L11 and L12 are generated independently of each other, their common occurrence is too rare. The observed numbers of systems containing L11, L12, and both levels (i.e., 2+2 systems) are [N_L11, N_L12, N_{2 + 2}] = [296, 95, 41]. We define the correlation C = N_{2 + 2}/N_L12 = 0.43 as the fraction of 2+2 systems among the L12 systems. Let N_L1 = 1747 be the number of L1 (outer) binaries in the sample. If the occurrence and detection of L11 and L12 were mutually independent, their respective probabilities would be N_L11/N_L1 = 0.169 and N_L12/N_L1 = 0.054. The product gives an estimate of 2+2 hierarchies as N_{2 + 2} = 16, to be compared to the actual number N_{2 + 2} = 41. The excess of 2+2 systems indicates correlation between levels L11 and L12. This effect was noted by Tokovinin & Smekhov (2002): the frequency of spectroscopic subsystems in the distant tertiary components of known triples is enhanced in comparison to the field stars, as though the multiplicity syndrome were contagious. Subsystems of levels L11 and L12 correlate because they were formed close to each other, in the same environment or even by the same process.

The excess of 2+2 systems is accounted for in the simulations by increasing by ₊ = 1.2 times for the secondary subsystems at L12, if the primary subsystem at L11 is already present. Otherwise, the frequency of L12 subsystems is multiplied by ₋ = 0.5; without such reduction of , the simulator overproduces the level L12. By playing with these two numbers, we adjust the total number of L12 systems and their correlation with L11. We also adopt β = 1.0 for L12, guided by the ML results (when β ∼ 0 is used for the secondary subsystems, there are more undetected L12 pairs and the number of hierarchies becomes larger, f_H = 0.15). We tried several other methods to simulate the observed correlation between L12 and L11, for example, by modifying the periods at L12 or by generating L12 with = 1 in the presence of L11. The method chosen here gives good results, but it is not unique. Simulations give the observed correlation C = 0.45 and the intrinsic correlation (before applying the detection filter) of 0.56.

A reasonably good match between the real and simulated samples is demonstrated by Table 3. To reduce statistical errors, we simulated the sample 10× larger, N = 48, 000, and scaled down the numbers. The simulated multiples successfully mimic the period distributions of the real sample and the plot in Figure 7. The median periods at levels L1, L11, and L12 in the simulated and real samples agree well.

The last column of Table 3 gives the system count in the simulated sample before applying the detection filter. They represent the true multiplicity corrected for observational selection. Remember, however, that this method of correction depends on the simulation recipe used here, which is not unique. Although the simulated and real statistics match quite well, this does not mean that another equally good simulator leading to somewhat different predictions cannot be devised. Interestingly, the apparent number of pure binaries remains almost unaffected by the observational selection, while the numbers of true and observed hierarchies differ substantially, especially at level L12. The simulations predict that a large number (∼150) of 2+2 systems in the FG-67 sample remain to be discovered.

The simulator overestimates the number of subsystems at high hierarchical levels L111, L121, L112, etc., because it does not account for tidal evolution in close inner binaries and selects their periods from the same lognormal distribution. We do not attempt to correct this effect and simply take it into consideration. Because of this, the simulated number of triples is smaller than observed (274 versus 290), while the numbers of high-order hierarchies are larger. The companion frequency f_C is slightly overestimated. However, the total number of hierarchies, which depends only on levels L11 and L12 and on their correlation, is simulated correctly.

The total number of hierarchical systems in the simulated sample is 630, leading to f_H = 0.132. It agrees with f_H = 0.129 estimated by the ML method (different experimental versions of the simulator produced f_H between 0.12 and 0.15). The fractions of L11, L12, and 2+2 systems in the simulated sample are 0.100, 0.073, and 0.041, respectively. The estimated fractions of single, binary, etc., systems (last column of Table 3) are f_n  =  [0.543, 0.323, 0.076, 0.043, 0.013], where the last number is the fraction of systems with five or more components. The calculated companion fraction is then f_C  =  0.65. As expected, it is larger than its ML estimate = 0.57 (the first line in Table 2). The actual companion fraction must be about 0.6.

It is truly remarkable that complex statistics of hierarchical multiple systems can be reproduced by the simulation recipe with only five adjusted parameters _i, ₊, and ₋, while the remaining parameters of the generating distribution x₀, σ, β are determined by the ML method and kept fixed. This supports the underlying idea of independent multiplicity. Therefore, the shorter periods of subsystems at levels L11 and L12 (compared to all binaries) can be explained by the dynamical truncation alone.

Figure 11 illustrates one aspect of the data that is not yet captured by the simulator. All binaries up to a certain period P are selected, and the observed fraction of inner pairs in multiples (levels higher than L1) among those binaries is plotted. The fraction of inner systems reaches ∼0.6 at P < 3 days and slowly drops to 0.18 with increasing period. The fraction of inner pairs in the simulated sample (with the detection filter applied), plotted as the dashed line, is quite similar to the observed one, except at P < 10 days, where it is much lower.

Zoom In Zoom Out Reset image size

Inner binaries in triple systems with mutually inclined orbits experience periodic modulation of their eccentricity by the Kozai–Lidov cycles. When the separation at periastron becomes comparable to the stellar radii, tidal friction absorbs the orbital energy and the period of the inner binary shortens to P_S ≲ 10 days. This mechanism, known as Kozai cycles with tidal friction (Eggleton 2006), works well when the period ratio P_L/P_S is not too high and the initial inner period P_S is not too long. It causes migration of inner periods from P_S ∼ 100 days to P_S ⩽ 10 days (Fabrycky & Tremaine 2007), while some inner binaries can merge.

Tidal migration of inner periods reflects in their period distribution (lower panel in Figure 6), although the addition of unknown fictitious periods to the histogram partially masks the effect. See also the cumulative period histogram at L111 in Figure 3, where 40% of periods are shorter than ∼3 days. The period distribution of spectroscopic binaries in the Hyades (Griffin 2012, Figure 68) clearly shows the depletion at P > 10 days and the peak at shorter periods. The non-monotonous period distribution was also found by Halbwachs et al. (2003). Tidal evolution explains why the median period at L111 is shorter than its simulated value and why the simulator overestimates the number of high hierarchical levels. Apparently, some close inner binaries merged.

We examined correlations between the mass ratios in triple systems. The 199 simple triples (levels L11 or L12 without further inner subsystems, with known periods and masses) are used. The data are affected to some extent by the observational selection, as no simple correction can be suggested in this case. We found no correlation between the mass ratios in inner and outer systems and between the mass ratio and the period. A similar conclusion was reached by Tokovinin (2008) from analysis of the biased catalog. For example, Figure 12 plots the ratio q_sys of the distant companion's mass to the mass sum of the inner binary. The median value of q_sys for outer binaries is 0.39, with no obvious dependence on the outer period P_L (when the sample is split into two halves at P_L = 10^6.7 days, the q_sys medians are 0.40 and 0.39). The median mass of the tertiary components is 0.73 M_☉; the median mass of the secondary components in all binaries is 0.70 M_☉. Therefore, there is no evidence that the distant tertiary components are less massive than the components of the inner binaries.

Zoom In Zoom Out Reset image size

Tokovinin (2008) noted that in 2+2 quadruples with P_L < 10^5.4 days both inner binaries have similar total masses and periods. The FG-67 sample contains only five such quadruples, all with q_{sys, out} > 0.5. The statistics are too small to make any further conclusions.

We checked whether twin binaries have any preference to be found in multiple systems. Half of the 96 binaries with P < 10 days are inner subsystems (see Figure 11). Among those 96 binaries, 18 are twins with q > 0.95, and 8 of those twins (also a half) are in inner subsystems. Therefore, the phenomenon of twin binaries is not directly related to the hierarchical multiplicity. The mass ratios of all binaries, including twins, seem to be independent of their hierarchical level.

One notable exception to the above statement concerns secondary subsystems at L12, which show a tendency to equal masses, β  =  1.3, according to the ML parameter fitting (Section 3.3). About half of those binaries (41 out of 95) also have primary subsystems at L11. The median mass ratio in those 41 primary subsystems is 0.75, while the remaining majority of L11 subsystems have median mass ratio of 0.65. This result does not rely on the ML analysis or completeness correction; it is a straightforward comparison between the primary subsystems whose outer companions are themselves either binary or single. Preference of equal masses is thus typical for both inner pairs in 2+2 hierarchies. This suggests that they could have formed differently, compared to other binary and triple stars.

The large sample of solar-type binaries with known detection limits confirms and strengthens known facts about their statistics (Duchêne & Kraus 2013). The parameters of the lognormal period distribution determined here (x₀ = 4.54, σ = 2.40) match well the results of R10 (x₀ = 5.03, σ = 2.28). The triangular function (2) fits the distribution of x = log P even better than a Gaussian. However, the period distribution of the inner subsystems is not lognormal; it has an excess at P < 10 days.

We find a nearly flat distribution of the mass ratio with the power index β ∼ 0 and show that it does not depend on the orbital period (with a caveat related to partially missing binary parameters). Duchêne & Kraus (2013) fitted β = 0.28 ± 0.05 to the data of R10, but found that β = 1.16 ± 0.16 for periods shorter than 10⁵ days. Considering the sensitivity of derived β to the treatment of selection effects (see Section 3.2), we believe that this claim is premature. Reggiani & Meyer (2013) derive β = 0.25 ± 0.29 for the two combined samples of dwarf binaries and argue for the independence of β on period and primary mass. The flat f(q) was found by Kraus et al. (2011) for the pre-main-sequence stars in Taurus more massive than 0.7 M_☉.

This study gives the first look at unbiased statistics of hierarchical systems with three or more components. Multiples are frequent (one in eight stars). The 13% frequency of hierarchies estimated here is essentially the same as the 12% found by Raghavan et al. (2010) in the 25 pc sample. However, they focused on companions to the primary stars and missed some secondary pairs. A few of those secondary subsystems are recovered here; more are yet to be discovered. On the other hand, these authors may have overestimated the number of high-order multiples. They list three multiples with five or more components (their Figure 24), while we found only five such systems in a sample 10× bigger. This discrepancy is explained by the better companion census of nearby stars, by the inclusion of uncertain subsystems (such as HD 186858 Aa, Ab) in the 25 pc sample, and by the small-number statistics.

Note that among the 11 quadruples in the 25 pc sample, 9 (or 82%) have the 2+2 hierarchy and only 2 have the 3+1 hierarchy. The last column of Table 3 predicts a 74% fraction of 2+2 quadruples, in agreement with the small statistics within 25 pc.

Fuhrmann (2011) gives estimates of f_M and f_H among F-type stars within 25 pc that agree well with our results. He claims that both fractions rise sharply with increasing primary mass (although with a low significance limited by the sample size). We do not confirm this trend: the observed (not selection-corrected) multiplicity fraction of 845 stars with M > 1.3 M_☉ within 67 pc is 38 ± 2%, not significantly different from the 36% in the full FG-67 sample.

The statistics of hierarchical multiples are well modeled by independent selection of subsystems from some generating distribution of periods, provided that the binary frequency is variable. This hints that the field stars were formed in different environments, binary-rich and binary-poor. Indeed, the binary frequency differs between star-forming regions (King et al. 2012). Moreover, some fraction of single stars come from disintegrated multiples and merged binaries, so the intrinsic frequency of the formation process is always larger than the multiplicity of a mature field population.

Variable multiplicity naturally leads to the enhanced probability of finding subsystems in wide binaries, in comparison to all targets. This tendency was noted by Makarov et al. (2008). The simulations faithfully reproduce the fraction of close binaries belonging to higher-order systems for the inner periods longer than 10 days (Figure 11); at shorter inner periods, the observed fraction is much higher than predicted, because of tidal evolution.

Kraus et al. (2011) already tried to model hierarchical multiplicity in a Monte Carlo simulation by producing inner subsystems with the generating period distribution and applying the stability cutoff. A similar approach was taken by Santerne et al. (2013) in their modeling of hierarchical multiples for evaluation of false positives in exoplanet detection. The multiplicity simulator developed here is anchored to the data and can be used in such studies with confidence.

Correlated presence of subsystems in both components of a wide binary is the new fact found in the course of this project (Riddle et al. 2014, in preparation). To reproduce this correlation by simulation, the frequency of secondary subsystems is enhanced in the presence of the primary subsystem and decreased otherwise. As a result, the majority of quadruple systems have the 2+2 hierarchy (67% observed, 74% predicted after correction of observational selection, 82% in the 25 pc sample). The fraction of quadruple stars is then larger than predicted by the independent multiplicity and comparable to the fraction of triples (4% and 8%, respectively). The multiplicity fractions f_n do not decrease in geometric progression of n, as suggested in the early works (Batten 1973) and repeated by Duchêne & Kraus (2013). Poor detection of subsystems in secondary components was not fully appreciated until recently.

The relatively large fraction of 2+2 systems suggests that they were formed by some special process. Other properties of 2+2 quadruples such as the similarity of components' masses and inner periods point in this direction, too (Tokovinin 2008). Both inner subsystems in the 2+2 hierarchies differ from the rest of the inner subsystems by their larger mass ratios (Section 4.5). On the other hand, when only one inner subsystem is present, it is preferentially associated with the most massive star, and its mass ratio is distributed uniformly, as in normal binaries. Therefore, there must be another, dominant formation channel that produces most binary and triple stars and, by extension, 3+1 quadruples.

Statistical work on multiple stars, including this one, is motivated by the desire to understand their formation. What have we learned from this study?

Hydrodynamical simulations (Bate 2012) reveal a complex process of binary formation and early evolution that shapes the multiplicity statistics. Stellar pairs form by fragmentation with initial separations larger than 10–100 AU (so-called opacity limit) and small masses. The components grow by accreting gas and, at the same time, migrate to smaller separations. The depletion of long periods (the right side of the "Gaussian") by dynamical interactions with the environment, combined with migration, shapes the observed period distribution (Korntreff et al. 2010). If the fragmentation generates periods x = log P longer than, say, x = 4 (a step function in x) and the period decay Δx is a random function (e.g., uniformly distributed), the convolution of those two distributions results in the linear decrease of binary frequency at x < 4, as actually observed (Figure 6). The approximately linear decrease of the period distribution at x < 4 is produced under a wide range of assumptions, provided only that the initial period distribution has a lower cutoff and that the migration is random.

Migration is essential in understanding the origin of tight hierarchical systems where three or four stars are packed in a small volume. If all inner binaries formed with initial periods x ⩾ 4, the initial outer periods must be even longer, whereas the present-day outer periods can be as short as x = 3 (Figure 7). Migration occurs therefore at all hierarchical levels. The inner binaries must migrate faster or form earlier than the outer binaries; otherwise, the multiple system becomes dynamically unstable and decays.

Most hierarchies seem to be assembled "from inside-out," by bringing together prefabricated inner pairs (or stars) to make the outer system. The success of modeling the observed statistics by independent multiplicity strongly supports this view. The mass ratios of inner and outer systems do not correlate; the masses of outer (tertiary) components are comparable to the masses of the inner components. The inner and outer periods are selected independently from the same distribution and are mutually related only by the dynamical stability (setting aside the tidal evolution). If hierarchies formed "from outside-in" by successive fragmentation of outer binaries, as envisioned, e.g., by Kraus et al. (2011), it would be difficult to explain the independent multiplicity.

Considering that many hierarchies are close to the stability limit, dynamical decay of multiple systems must occur sometimes. It can happen when the outer binary migrates faster than the inner one or when other stars in the cluster disrupt the multiple or exchange components with it. However, multiple systems produced by N-body dynamical interactions (or surviving them) have characteristics that do not match real hierarchies: high eccentricities, low masses of distant components, and moderate ratios P_L/P_S. The observed multiplicity statistics indicate that dynamical processes play only a secondary role. Most hierarchical stellar systems in the field are not the surviving remnants of primordial clusters.

The 2+2 hierarchies could be formed in a special way, e.g., by inelastic collision of two cores that prompts their fragmentation into subsystems and, at the same time, keeps these two pairs together on a common wide orbit by dissipating the energy. Such impulsively triggered multiple-star formation was envisioned by Whitworth (2001); see his Figure 2. Further evolution will depend on the relative orientation of the inner and outer orbits. If they are not aligned, the inner binaries will accrete gas with a misaligned angular momentum, which will accelerate their migration to short periods. In the opposite case of co-aligned orbits, the accretion will keep the inner binaries wide, and a quadruple with a moderate P_L/P_S ratio and nearly co-planar orbits, resembling ε Lyr, can be formed. Another striking aspect of this prototypical 2+2 quadruple—similar masses of all four stars—matches the tendency to equal component masses in both inner subsystems of 2+2 hierarchies, emerging from this study.

On the other hand, the majority of the inner subsystems have β ∼ 0, like all binaries, while their tertiary components are single. Their formation process should therefore be different from the formation of 2+2 hierarchies. The tertiary component could have formed around the initial binary (or joined it) at a later stage. The cascade can continue further by adding yet another star in a 3+1 hierarchy.

The 2+2 formation channel probably produced more quadruples than actually observed because some subsystems merged, leaving triples and binaries, while some quadruples decayed dynamically. A certain fraction of binaries and triples in the field could be born as 2+2 systems. Also, merging of the inner binaries (natural extension of migration) can possibly help in understanding why the outer periods in triple stars are longer than 10³ days.

It should be stressed that not all close binaries belong to multiple stars. This is true only for pairs with P < 3 days, which could not be produced without tidal migration during their main-sequence life. At longer periods P ⩾ 10 days, the frequency of wide tertiary companions is comparable to the frequency of such companions in the full sample (Tokovinin et al. 2006).

Figure 14 illustrates the formation channels and the fractions of different hierarchies in the field. The fragmentation produces single and binary stars. Cascade "inside-out" fragmentation generates 2+1 triples and 3+1 quadruples. Unstable hierarchies decay dynamically into single and binary stars; some close binaries merge. The alternative 2+2 process creates quadruple stars with components of similar mass. Merging of inner pairs in the 2+2 quadruples contributes to the populations of triples and binaries. About 1% of targets with more than four components could be produced by the combination of these elementary processes.

Zoom In Zoom Out Reset image size

There are various ways to improve and continue the present study. The determination of unknown periods and mass ratios of close binaries is one obvious project to pursue. The prediction of the large number of undiscovered subsystems in the secondary components should be tested observationally. It will be interesting to study relative motions in resolved triple and quadruple systems to get a better idea of the relative orbit orientation.