Architecture of Hierarchical Stellar Systems and their Formation

Accumulation of new data on stellar hierarchical systems and the progress in numerical simulations of their formation open the door to genetic classification of these systems, where properties of a certain group (family) of objects are tentatively related to their formation mechanisms and early evolution. A short review of the structure and statistical trends of known stellar hierarchies is given. Like binaries, they can be formed by the disk and core fragmentation events happening sequentially or simultaneously and followed by the evolution of masses and orbits driven by continuing accretion of gas and dynamical interactions between stars. Several basic formation scenarios are proposed and associated qualitatively with the architecture of real systems, although quantitative predictions for these scenarios are still pending. The general trend of increasing orbit alignment with decreasing system size points to the critical role of the accretion-driven orbit migration, which also explains the typically comparable masses of stars belonging to the same system. The architecture of some hierarchies bears imprints of chaotic dynamical interactions. Characteristic features of each family are illustrated by several real systems.


INTRODUCTION
Formation of stars and planets is at the forefront of astronomical research, being driven by the need to understand our origins. Despite tremendous progress, modern theory does not yet have enough predictive power to explain statistics of nascent stars from the first principles. Star formation is a complex and chaotic phenomenon involving a wide variety of inter-wined factors and physical processes. Stars are born in groups, and many stars are bound in binary and multiple systems; properties of these systems are related to their formation and early evolution, and their study advances the general knowledge of star formation. In this review, I focus on multiple systems containing three or more stars and attempt to link their architecture to potential formation scenarios. While the distributions of periods, eccentricities, and mass ratios of binary stars are generally recognized as important tracers of the star formation mechanisms, hierarchical systems bring additional information encoded in the period ratios, mutual orbit orientation, and masses of the components.
The observational knowledge of multiplicity has progressed from a crude description of small and complete samples of nearby solar-type stars (Duquennoy & Mayor 1991;Raghavan et al. 2010) to a more detailed assessment of multiplicity statistics in different mass ranges and environments (Duchêne & Kraus 2013;Moe & Di Stefano 2017). Large multiplicity surveys, focused mostly on binaries, have revealed recently a strong dependence of multiplicity on mass, and also discovered that multiplicity statistics depend on the density and temperature of the formation environment and on its metallicity. The idea that statistics of binaries (and, by extension, stars and planets) are not universal is firmly taking root in astronomy and guides the research in new Electronic address: andrei.tokovinin@noirlab.edu

directions.
Theoretical studies of star formation also advance rapidly. Modern hydrodynamical simulations with increasing resolution and improved treatment of radiative physics and magnetic fields provide an unprecedented and detailed view of the formation and early evolution of stellar systems. For example, Bate (2019) created animations of multiple-star formation in a massive collapsing cloud and pointed out that most close binaries result from dissipative gas-assisted capture of two stars that formed either independently in the same cluster or in close proximity (in the same filament). Other mechanisms (disk instability and dynamical interactions) also play a substantial role, and several mechanisms may combine to create a given system. This emerging picture is further confirmed in the simulations by Lee et al. (2019) and Kuffmeier et al. (2019) who studied smaller clouds with a higher resolution. In these simulations, companions form sequentially, approach the main star, and migrate to closer orbits owing to the continuing gas accretion and dynamical friction.
New knowledge is usually produced by the joint advances of theory and observations, their interaction and confrontation. However, the complexity of star formation, its dependence on the conditions, and the multitude of processes involved, stand in the way of direct comparison between observations and theoretical predictions. Classification of objects helps to systematize available data and is a necessary step to their interpretation.
Its textbook examples in astronomy are galaxy types and variable stars. Classification is also widely used in biology, relating external (observable) characteristics of species to their origin. Biological classification was developed empirically well before the genetic code became readable. Here an attempt is made to classify stellar hierarchies based on their architecture and to relate these families to the formation processes, adopting the 'genetic' approach. However, reliable models of multiple-star formation are not yet available, their predictions are mostly qualitative. Classification of stellar hierarchies proposed here remains therefore intuitive, reflecting the current incomplete understanding of the formation mechanisms. I hope that it will be revised and improved in the future.
Historically, stellar hierarchies were first studied on the individual, case-by-case basis because only a few close systems were known, while wide visually resolved triples move too slowly on human time scale. In his book, Batten (1973) recognized the potential importance of hierarchies for understanding formation of close binaries. This idea was based on the frequent occurrence of close binaries in higher-order hierarchies, noted by the observers. The next step was made by Fekel (1981) who assembled the first list of 35 spectroscopic binaries with tertiary companions in relatively tight (outer periods <100 yr) hierarchies and attempted to connect their statistical properties with the formation mechanisms known at that time. One of the best historic ways of finding hierarchies was by detecting eclipses in visual binaries; a list of 80 such systems was published by Chambliss (1992). On the other hand, several wide physical hierarchies were documented in the catalogs of nearby stars.
Trying to measure stellar masses by a combination of visual and spectroscopic orbits in the early 1990s, I accidentally discovered several triples and started to collect data on such systems, being convinced of their astrophysical importance by the works of Batten, Fekel, and others. The result was the first compilation of known physical hierarchical systems covering the full range of periods -the Multiple Star Catalog, MSC (Tokovinin 1997). In parallel, systematic spectroscopic monitoring was initiated to discover more hierarchies; its results are summarized by Tokovinin & Smekhov (2002). I extended the 25-pc sample to 67 pc with the aim to study unbiased statistics of hierarchies with solar-type components (Tokovinin 2014). After completion of the 67-pc project, I continued to update the MSC.
It is instructive to see how the knowledge of hierarchical multiplicity in the solar neighborhood progressed with time. In the 22-pc sample of Duquennoy & Mayor (1991), the estimated fraction of hierarchies among all systems was 5%. The more detailed survey of the 25-pc sample by Raghavan et al. (2010) boosted this fraction to 13%, while the total multiplicity did not change. This increase was due to the discovery of additional subsystems in known binaries. The latest update on the 25-pc sample by Hirsch et al. (2021) gives the 17% fraction of hierarchies, larger than in the earlier studies (new subsystems were discovered mostly by high-resolution imaging). According to these authors, only 2/3 of non-single stars are pure binaries, while the remaining 1/3 of the systems host three or more stars. Hierarchical multiplicity is by no means rare among solar-type stars, and it increases with mass; almost all massive stars are at least triple (Sana 2017).
I begin the review by introducing hierarchical systems in section 2, defining their architecture and the corresponding parameters that serve for classification. Known statistical trends are briefly covered. Then in section 3 the formation mechanisms of stars, binary stars, and multiple systems are schematically outlined, focusing on the relation between formation scenarios and the archi-
tecture of the resulting products. The proposed classification of hierarchies is presented in section 4, where each group is illustrated by real systems. Summary and outline of future progress in section 5 close this review.

PROPERTIES OF HIERARCHICAL SYSTEMS
2.1. The Multiple Star Catalog This study uses data on known hierarchical systems collected in the MSC; its current update is described in Tokovinin (2018b). The catalog is a compilation of random discoveries and surveys. It is heavily affected by observational selection which distorts relative frequency of observed hierarchies compared to their intrinsic distribution. For example, systems containing eclipsing binaries are over-represented, being easier to discover. Nevertheless, the MSC holds a large (almost 3000) sample of real multiple systems, allowing to distinguish characteristic patterns in the multi-dimensional space of their parameters.
Periods and separations of binary and multiple systems span many orders of magnitude, and their other properties such as masses and structure are also very diverse. The diversity of stellar hierarchies surpasses the diversity of exoplanetary systems. Nowadays, exoplanets are successfully classified into several groups (e.g. He et al. 2020). Here a similar effort to classify the diverse population ('zoo') of stellar hierarchies is undertaken.
The MSC is regularly updated. Its latest version is posted online 1 and is available as the Vizier catalog J/ApJS/235/6 2 . Full information on the multiple systems mentioned in this paper can be retrieved from the MSC. Each system has a unique code based on its J2000 coordinates, e.g. 11551+4629. Similar codes are adopted in the Washington Double Star Catalog, WDS (Mason et al. 2001), and I call them 'WDS codes' for brevity. Each system in the MSC has a grade; grades 3, 4, and 5 have at least three reliably identified components, while grades 1 and 2 (questionable and rejected systems) are not considered here.

Types of Hierarchy
A dynamically stable hierarchical system can be decomposed into a combination of binaries. Its structure is described by a binary-tree graph, also called mobile diagram. Figure 1 introduces basic types of hierarchy. A triple system consists of the outer (wide) pair which is at the root of the hierarchy, and an inner close pair. The hierarchy can be represented in various equivalent ways, e.g. by a graph, by numbers or levels (outer binary is level 1 and the inner binary is level 11 or 12 if it belongs to the main or secondary component, respectively),  by nested brackets (Eggleton & Tokovinin 2008), or by reference to the parent component (a subsystem Aa,Ab has component A as parent). The first MSC version (Tokovinin 1997) coded the hierarchy by levels, and the current MSC uses references to parents, coding the hierarchy by the triads (primary, secondary, parent). The outermost subsystem (root of the tree) has asterisk (*) as parent. Links to parent is a flexible scheme that, on the one hand, defines the hierarchy and, on the other hand, allows easy modifications in response to discoveries of new subsystems.
Quadruple systems in Figure 1 can have two possible hierarchies. The 2-tier hierarchy (two close binaries orbiting each other) is called a 2+2 quadruple. The 3-tier hierarchy where a close binary Aa,Ab has a tertiary component B and this triple AB is orbited by another more distant star C, is called a 3+1 or 'planetary' hierarchy. Indeed, it resembles a planetary system if the central star has a substantially larger mass than its companions. Among solar-type stars within 67 pc, 2+2 quadruples are ∼4 times more frequent than 3+1 quadruples (Tokovinin 2014). The proportions of single, binary, triple, quadruple (and higher-order) systems in this sample are 54:29:12:5, so a 0.17 fraction of the total population are hierarchies. Hirsch et al. (2021) confirm this fraction for the small but well-studied 25-pc sample Complex hierarchies containing more than four stars can be viewed as combinations of elementary binaries, triples, etc. Figure 2 shows the structure of the unique 5-tier system 65 UMa (11551+4629); no other 5-tier hierarchies are known so far. It has an almost planetary-type structure, except that the outermost component D is itself a binary. A hierarchical system with N tiers can have no more than 2 N components if all available levels are filled. This condition is fulfilled in 2+2 quadruples. The 65 UMa system has only 7 components, while the maximum possible number of stars in a 5-tier hierarchy is 2 5 = 32. This situation is typical: only a fraction of available levels are usually filled.
With a few exceptions, we cannot be sure that all components and subsystems in a given hierarchy are discovered. For example, in 65 UMa the visual components B or C can themselves be yet undetected close binaries or even triples. This means that the current description of the hierarchy may be incomplete and it will change with new discoveries. A close spectroscopic binary Aa,Ab with a distant visual companion B looks like a typical triple system with a 2-tier hierarchy. However, further study might reveal additional component Ac revolving around A at close separation, converting this triple into a 3+1 quadruple. After this discovery, the inner binary Aa,Ab moves from level 11 to level 111, and the new subsystem Aab,Ac occupies the intermediate level 11.

Main Parameters of Hierarchies
A hierarchical system can be decomposed into several nested binaries. Motion in each of those binaries is approximately described by a Keplerian orbit. Dynamical interactions in the system change parameters of those instantaneous (osculating) orbits with time, but usually these changes are small and slow. Therefore, the standard orbital parameters of a binary -period P , semimajor axis a, eccentricity e, and masses of the components M 1 and M 2 -also serve to characterize hierarchies. A hierarchy with N components has N −1 subsystems and, hence, orbits. The masses in the outer orbits refer to the combined masses of the inner subsystems. For example, a triple system containing three equal stars will have the inner mass ratio q in = M 2 /M 1 = 1 and the outer mass ratio q out = M 3 /(M 1 + M 2 ) = 0.5.
Apart from the two sets of individual orbital elements, a hierarchical triple has additional parameters, namely the period ratio P out /P in and the mutual inclination Φ -the angle between the angular momentum vectors of the inner and outer orbits ranging from 0 • for coplanar and co-rotating systems to 180 • for coplanar and counterrotating ones; perpendicular orbits have Φ = 90 • . Unfortunately, mutual inclination is measured only for a small subset of known hierarchies owing to the observational limitations. These mutual parameters define the degree and character of dynamical interaction between the inner and outer orbits. More complex hierarchies have even more such mutual parameters. Figure 3 compares the inner and outer periods at adjacent hierarchical levels. A triple system gives one point in this plot, a quadruple system two points, etc. To reduce somewhat the observational selection, only 1820 systems within 200 pc are selected from the current MSC. The symbol colors in the panels distinguish the systems by the inner mass ratio q in and by the total system mass. The solid line marks the limit of dynamical stability P out /P in = 4.7 (see section 3.7) and the dashed line marks P out /P in = 100. Periods of wide binaries are estimated approximately (within a factor of ∼3) from their projected separations using the Kepler's third law, so some long-period triples appear below the stability limit owing to inaccurately known periods. This figure illustrates the wide range of the periods and their ratios; no preferred ratio is apparent.

Statistical Trends
Roughly speaking, a triple system is just a combination of two binaries, inner and outer. The binary statistics are relatively well established, at least for solar-type stars. Can the statistics of triple systems be derived from the binary statistics by assuming a random choice of the inner and outer pairs from the general binary population and keeping only dynamically stable combinations? This simple approach implies that the outer and inner pairs in a hierarchy could be formed independently of each other, and it is called independent multiplicity model (IMM). Indeed, the properties of all binaries and of the binaries belonging to hierarchies are similar. For example, the frequency of wide visual companions to spectroscopic binaries (except the closest ones) is not very different from this frequency for single stars. Conversely, the rate of close binaries among components of wide visual pairs is comparable to the overall rate of close binaries in the field. With some caveats and exceptions, the IMM roughly matches the statistics of solar-type hierarchies (Tokovinin 2014), offering a crude but convenient mathematical recipe to compute the distribution of hierarchies over periods.
Let ǫ be the probability that each star is a binary. The IMM postulates that the frequency of triples should be proportional to ǫ 2 , quadruples to ǫ 3 , etc., because subsystems are chosen independently from the same parent distribution. However, the value of ǫ that matches the binary fraction in the solar neighborhood predicts fewer triples than actually observed. This discrepancy can be fixed by assuming that the field population is a mixture coming from different environments with different binary frequencies ǫ. Then the relative proportions of single, binary, triple, etc. stars can be successfully reproduced by the IMM with a suitable choice of the mean ǫ and its variance (Tokovinin 2014). In this paradigm, the binaryrich constituents of the field contribute most binaries and hierarchies, while the binary-poor environments supply to the field mostly single stars and a minor fraction of binaries. The idea of variable ǫ is supported by the observational evidences of environmental multiplicity depen-dence. To give a few examples, Deacon & Kraus (2020) found a measurable difference in the wide-binary frequency between open clusters and moving groups. In the globular cluster NGC 3201, Kamann et al. (2020) measured a binary frequency of 23.1±6.1% in the first population of stars and 8.2±3.5% in the second-generation population that formed later in a more dynamically active environment.
The IMM implicitly suggests (although, strictly speaking, does not require) independent formation of the subsystems, therefore relative orientation of the inner and outer orbits in a triple should be random; there should be as many co-rotating triples as counter-rotating ones. This assertion can be verified by comparing the numbers of apparently co-and counter-rotating triples, assuming their random orientation with respect to the observer and equal chances of discovering co-and counterrotating systems. This method works only for resolved (visual) triples in the solar neighborhood, where the sense of orbital motion in both subsystems can be determined. This is possible even when only a short fraction of the orbit is observed, i.e. for periods of hundreds and thousands of years. An excess of apparently co-rotating triples was found by Worley (1967) and later confirmed by Sterzik & Tokovinin (2002) and Tokovinin (2017b). Therefore, the orbits of resolved visual triples are not randomly oriented. A more detailed study shows that the degree of orbit alignment depends on the outer projected separation s out (in au) and on the ratio of periods or separations. Compact visual triples with s out < 50 au have approximately aligned orbits; the degree of alignment decreases with separation, and at s out > 10 3 au the relative orbit orientation becomes random.
Yet another deviation from the IMM in the 67-pc sample of solar-type stars is the excess of 2+2 quadruples compared to their fraction predicted by the model (Tokovinin 2014). This means that the occurrence of inner subsystems in both components of a wide binary is correlated. Such a correlation was suspected from the study of triple systems consisting of a wide binary and an inner spectroscopic pair (Tokovinin & Smekhov 2002). Radial velocity (RV) monitoring of the tertiary components in these triples revealed that many tertiaries have a variable RV, i.e. they also contain subsystems and actually these systems are 2+2 quadruples. In other words, presence of a subsystem in one component of a wide binary enhances the chance of finding a subsystem in the other component. To account for this observational fact, the multiplicity model of the 67-pc sample postulates correlated occurrence of the inner subsystems in both components of a wide pair.
Finally, the fact that very close (periods under 10 days) binaries are more often found within hierarchies (as inner subsystems) than among other stars obviously contradicts the IMM. The relation between close binaries and hierarchical systems has been suspected a long time ago, e.g. by Batten (1973), and now it is well documented. Tokovinin et al. (2006) searched for tertiary components to spectroscopic binaries and found that their fraction anti-correlates with the binary period: it is close to 100% for binaries with periods under 3 days and drops progressively to 40% at periods of 10 days and longer; the latter is similar or even lower than the average fraction of wide companions to all stars and proves that close binaries without tertiary components certainly exist. On the other hand, very close binaries prefer to be accompanied by other stars, rather than live alone. Hwang et al. (2020) found that 14.1±1% of contact binaries have wide companions, while their fraction for all stars is 4.5%. They used the Gaia catalog and made a simplifying assumption that stars with variable fluxes are contact binaries (most of them are).
Inverting the argument, we may say that components of wide binaries contain close pairs more frequently than other stars in the same population. This trend possibly extends to subsystems with periods longer than 10 days. Deacon & Kraus (2020) studied wide binaries in the field and in several nearby open clusters and found that the probability of finding a close subsystems is two times larger for the components of wide binaries than for the average stars in the same cluster. This effect can be partially explained by the variable ǫ: if wide binaries come preferentially from the binary-rich population, the fraction of subsystems in their components should be also larger than on average. Relation betwen close and wide binaries in the Taurus association has been studied by Joncour et al. (2017).
To take a closer look at the inner subsystems, I selected from the MSC 425 inner subsystems with primary components of less than 1.5 M ⊙ mass within 200 pc distance and with known orbits (mostly spectroscopic). Their statistics are shown in Figure 4. The period-eccentricity plot in the left panel resembles similar plots for all spectroscopic binaries (e.g. figure 14 in Raghavan et al. 2010). Most pairs with inner periods shorter than 10 days have circular orbits owing to the dissipative effect of tides, although there are exceptions. On the other hand, at longer periods all orbits are eccentric (at P in ∼ 1000 days, circular orbits are found again). Twin subsystems with q in > 0.95 (green symbols) do not stand apart from other pairs in this plot. The dotted line indicates the regime where separation at periastron is small enough for tidal interaction; binaries located around this line will evolve to circular orbits with shorter periods (Meibom et al. 2006;Moe & Kratter 2018).
The right panel of Figure 4 illustrates the lack of correlation between period and mass ratio in the inner subsystems. The mass ratios are measured directy for doublelined binaries and are the minimum values for single-lined binaries. Most twins with q in > 0.95 have P in < 30 days (the same trend exists for all spectroscopic binaries, not just for the inner subsystems). The blue line shows the cumulative distribution of periods. One notes many contact systems with P in ∼ 0.3 days and a relative deficit of periods between 0.5 and 2 days. Contact binaries are easier to discover by eclipses, and their relative abundance in the MSC can be a selection effect.
Most remarkably, the period distribution has no features around P in ∼ 10 days. Dynamical evolution of  misaligned triples involving tides (see section 3.7 and Naoz 2016) affects mostly inner periods of 10-100 days, shortening them to P in < 10 days. As a result, the number of subsystems in the 10-100 days range is reduced, and they should pile up at periods just below 10 days, as predicted by Fabrycky & Tremaine (2007). An apparent excess of inner subsystems with P in < 10 days noted by Tokovinin & Smekhov (2002) was likely caused by the selection (easy discovery of eclipsing subsystems), but modern, larger data do not confirm it. Absence of the tidal signature in the distribution of the inner periods suggests that the mechanism studied by Fabrycky & Tremaine could not be a dominant factor in the formation of close binaries, echoing the conclusions of Moe & Kratter (2018) who determine that only a minor fraction of close binaries coud result from the tidal evolution within hierarchies, while most of them are products of early migration. Our brief overview of statistical trends would be incomplete without considering the mass ratios. Figure 5 compares the mass ratios q in = M 2 /M 1 in the inner subsystems with the outer mass ratios q out = M 3 /(M 1 +M 2 ) for the MSC subset within 200 pc. Only the inner subsystems at the lowest hierarchy levels (simple pairs of stars without subsystems) are selected. Overall, there is no obvious correlation between the mass ratios, although the large number of inner subsystems with similar-mass components, q in ∼ 1 (twins), is notable. Among the 1783 inner systems in this plot, 483 (27%) belong to the secondary component of the outer pair (green symbols), while the remaining majority are found in the primary component (red symbols). However, subsystems in the secondary components are more difficult to discover, hence their observed fraction is only a lower limit.
A system where the masses of the tertiary component and the inner primary are equal, M 3 = M 1 , can be called outer twin, and in such case q out = 1/(1 + q in ). This relation is shown in Figure 5 by the black line. When both inner and outer mass ratios are close to one, the system can be called double twin. When the tertiary star is the most massive one, M 3 > M 1 , the subsystem belongs to the secondary component of the wide pair, and these points, mostly located to the right of the black line, are plotted in green. However, the distinction between pri-mary and secondary components of a wide pair does not always correlate with the mass: it can be modified for evolved stars, by the mass transfer, and when the tertiary component also contains a subsystem (i.e. in 2+2 quadruples). This latter situation, marked by squares, explains why some points appear on the wrong side of the dividing line. A slight concentration of points to this line, noted earlier (Tokovinin 2008), would imply a preference of outer twins with M 3 ≈ M 1 ; however, this effect does not appear to be statistically significant. The plot features a large number of hierarchies where all three stars have comparable masses (q in ≈ 1, q out ≈ 0.5). Remember however that the MSC data are burdened by the observational selection, and hierarchies with comparable masses are discovered more easily.

FORMATION OF MULTIPLE STARS
In this Section, physical processes leading to the formation of binary stars are briefly reviewed and extended to the formation scenarios of stellar hierarchies.

Physics of Star Formation
Stars form by hierarchical collapse of giant molecular clouds Vázquez-Semadeni et al. (see the review by 2019). It starts with slow accumulation of cold molecular gas on large spatial scales caused by the Galactic spiral structure and colliding flows in the interstellar medium. The collapse is highly inhomogeneous; it proceeds by creation of two-dimensional (sheets or slabs) and one-dimensional (filaments) structures. During collapse, the density of the cold gas increases and the Jeans mass decreases, leading to fragmentation in a cascade, from larger to smaller spatial scales ( Figure 6). The fragmentation cascade stops when the gas becomes opaque and heats adiabatically, increasing the Jeans mass; the smallest scale of the cascade creates stellar embrios (protostars). This co-called opacity limit to fragmentation sets both the initial mass of the protostars and their minimum separation, on the order of 10 au (Larson 1972). At low metallicity, the gas contains less dust, decreasing the opacity and the opacity limit to fragmentation and helping to form closer binaries. This explains the increase of the close-binary fraction at low metallicity found by Moe et al. (2019). Interestingly, wide binaries exhibit an opposite trend:   Gravitational physics dictates that collapse at small scales proceeds much faster than at large scales (Larson 2007), so when the smallest structures form, the larger structures still collapse. As a result, each level in a hierarchical collapse continues to accrete from the upper levels: protostars accrete from clumps, clumps from the filament, etc. Vázquez-Semadeni et al. (2019) also point out that star formation continues on a relatively long time scale (a few Myr) defined by the outermost structure. The final stellar masses result from the accretion of gas on time scales much longer than the free-fall time of their natal envelopes and may be unrelated to the initial masses of the small-scale clumps from which the protostars form. Continued accretion on to newly formed binaries shrinks their orbits Tokovinin & Moe 2020). The origins of the stellar mass distribution (the initial mass function, IMF, and the system mass function, SMF) and its relation to multiplicity are reviewed by Lee et al. (2020). The SMF can be explained as a consequence of the mass distribution of prestellar cores (if a given fraction of the core mass is converted into stars), or as a result of competitive accretion. Clark & Whitworth (2021) successfully model the observed SMF by assuming that turbulent fragmentation produces low-mass seeds at some rate and these seeds subsequently grow by competitive accretion. In their model, only 0.23 fraction of the gas mass is consumed by forming the seeds, the rest is accreted later.
The specific angular momentum of molecular gas exceeds by several orders of magnitude the specific angular momentum of a protostar rotating at breakup speed. The mass growth of protostars by accretion is possible only when the angular momentum is extracted from the infalling gas. Given that the angular momentum transport needs some moving mass, accretion cannot proceed without ejecting a fraction of gas needed to carry away the angular momentum excess. The inevitable relation between accretion and ejection is manifested by the jets and outflows from young stars.

Binary Star Formation
Hierarchical systems consist of nested binaries, so mechanisms of binary formation and early evolution are essential for understanding the origin of multiples. Large-scale simulation of collapse by Bate (2019) show that binary stars form by a variety of mechanisms acting alone or in combination. Evolution of newly formed binaries in a dense collapsing cluster changes their parameters further by dynamical interaction with gas and neighboring stars. Therefore, the statistics of the resulting binary population do not directly relate to the output of elementary formation mechanisms. The two basic channels of binary-star formation illustrated in Figure 7 are disk instability and core fragmentation.
Disk instability (DI). As already noted, star formation is impossible without shedding the excessive angular momentum of the infalling gas. If the angular momentum transport is slower than its influx, the accreted gas accumulates in a disk around the protostar. A massive disk can become unstable to fragmentation, forming one or several companions around the nascent central protostar , and references therein). The opacity limit to fragmentation dictates that the initial separations of binaries formed by the DI mechanism are larger than ∼10 au.
Disks sufficiently massive to become unstable are more likely to exist around massive stars ). However, accretion is known to be highly variable (episodic). A burst of accretion increases the disk mass temporarily and may lead to the disk fragmentation even when the average disk is too small to fragment. The opacity of the gas plays an important role here, and the DI is enhanced at low metallicity (Moe et al. 2019).
A binary formed by DI has initially a very small mass ratio. Continued accretion on to such a binary increases masses of both components, but the secondary component grows faster and the mass ratio always increases. At the same time, the orbit shrinks through interaction with the circumbinary disk (Heath & Nixon 2020). The accretion-driven binary migration is a complex process depending on several factors, such as orientation of the binary orbit relative to the angular momentum of the infalling gas, size of the inner circumstellar disks, gas temperature, etc. Statistical modeling of accreting binaries based on simplified prescriptions can reproduce the overall properties of the close-binary population (Tokovinin & Moe 2020). The final mass ratio of a DIformed binary correlates with the time of the companion's formation: binaries that formed early tend to have larger mass ratios and develop a sub-population of twins  with nearly identical masses. Conversely, companions formed by the end of the mass-assembly period have little chance to grow. The DI mechanism explains the large fraction of close binaries found among massive stars by the large amount of gas needed to form such stars, favoring binary formation and migration. Furthermore, strong migration can lead to mergers of binary components, which is more likely for massive binaries; early binary mergers can produce the most massive stars. Core fragmentation (CF) is probably the dominant binary formation mechanism, being a direct consequence of the hierarchical collapse. The last stage of the collapse is set by the opacity limit to fragmentation (∼10 au) which also defines the minimum separation of binaries that can be formed by CF. Numerical simulations of an isolated cloud collapse driven by internal turbulence show how several protostars usually are born from the individual over-densities (Offner et al. 2010;Lee et al. 2019). Given that the gravitational collapse creates twoand one-dimensional structures, the protostars can form in linear configurations along filaments. Pairs of neighboring protostars may end up in bound binary systems. Over-densities are usually produced by gas compression in converging flows, where the kinetic energy of the colliding flows is mutually canceled. Relative velocities of the clumps in a filament are therefore smaller than the typical large-scale gas velocities (Kuffmeier et al. 2019). Two protostars born in the same filament can be gravitationally bound from the outset if their relative velocity is less than the escape velocity (Tokovinin 2017b).
Simulations of Bate (2019) show that even originally unbound protostars in a cluster can end up in a binary, while the excessive kinetic energy of their encounter is dissipated by the gas. Dissipative capture is the dominant mechanism of binary formation in these simulations. It should also work in isolated collapsing filaments. Two neighboring protostars approach each other while falling to the common center of mass. At the same time, each protostar continues to attract and accrete gas. When the two protostars become sufficiently close to each other, their envelopes interact, the kinetic energy is dissipated, and they form a gravitationally bound pair with a common envelope, even if their initial relative velocity exceeded the escape velocity. The initial size of such binaries should be comparable to the size of their disks or envelopes, on the order of a few hundred au or less.
A binary pair formed by CF continues to evolve and migrate as it accretes more gas . This is similar to the evolution of a DI-formed binary, except that the initial spatial scales can be larger and the motions of the gas relative to the binary are more likely to be chaotic, as happens in nascent clusters. The end products of the CF and DI mechanisms can be very similar. In a dense cluster, a competing mechanism of binary evolution are dynamical interactions with other stars or binaries. They modify the orbits, can disrupt the binary, and often involve exchanges of the components (Bate 2019).
Binaries help to form stars by storing the angular momentum of the infalling gas in their orbits. Sterzik et al. (2003) noted that the specific angular momentum of cores is comparable to the typical angular momentum of binaries. According to their logic, low-mass binaries originate from small cores and have correspondingly small orbits, matching the observed trend in the binary separations vs. mass. If the size of binary orbits is indeed determined mostly by the core angular momentum, we can talk about rotationally-driven core fragmentation in the spirit of the early collapse simulations by Larson (1972). Fragmentation of disk-like or bar-like rotating structures also happens in the modern cluster simulations by Bate (2019), as evidenced by movies provided in the supplementary material to that paper.
The most efficient way to store the excessive angular momentum in a binary orbit is by forming two fragments of comparable masses on a near-circular orbit, because this configuration corresponds to the maximum specific (i.e. per unit mass) orbital angular momentum. As if by coincidence, these characteristics are typical for the lowest-mass binaries (Dupuy & Liu 2011). On the other hand, formation of companions to low-mass stars by DI is unlikely . Therefore, fragmentation of isolated low-mass cores might be the dominant mechanism of forming low-mass binaries. Rohde et al. (2021) simulated a large number of fragmenting cores of 1 M ⊙ total mass and found that the resulting binaries have a strong preference to equal masses, including many twins. The typical mass of the components, ∼0.3 M ⊙ , corresponds to M-type dwarfs (only a half of the initial core mass was converted into stars). In these simulations, episodic outflows are shown to be important, Figure 8.
Formation of a triple system by sequential disk fragmentation and migration according to the model of Tokovinin & Moe (2020). The lines show the fraction of mass in each star vs. fraction of total accreted mass. The squares and asterisks show the inner and outer semimajor axes, respectively, in relative units. The first binary (green) forms early, at 0.05 accreted fraction, and migrates from 30 au to 0.05 au, reaching the final inner masses of 0.80 and 0.75 M ⊙ (a twin). The tertiary component (blue) forms at 178 au separation (outside the plot limit) when 0.25 fraction of the mass was already accreted and migrates to 23 au, reaching the mass of 1.18 M ⊙ and becoming the primary component of the triple. In most cases, however, the tertiary forms later and does not outgrow the inner two stars.
decreasing the average final masses and yielding a more realistic multiplicity statistics. On the other hand, variations of the core radius, virial parameter, and turbulence spectrum had little effect on the resulting masses and multiplicity. These simulations also produced a substantial number of triples, some of which were dynamically unstable.
The existence of close (spectroscopic) low-mass binaries and the large fraction of the low-mass twins suggest that the CF-formed binaries also migrate. The toy model of accretion-driven binary migration proposed by Tokovinin & Moe (2020) is equally applicable to the CFformed binaries. When a CF-formed binary accretes substantial mass, the resulting close pair has the same properties as the DI-formed close binary, erasing the signature of its initial formation mechanism. One may even wonder whether the distinction between the DI and CF mechanisms of binary formation is as meaningful as it appears at the first sight.

Formation Scenarios of Hierarchical Systems
Hierarchical systems with three or more components can be formed by the combination of elementary mechanisms outlined above and further modified by the Nbody interactions. The sequence of events producing a hierarchy is called formation scenario. Potential scenarios are outlined below and summarized in Table 1. The proposed scenarios are hypothetical and qualitative. Examples of these scenarios can be found in the numerical simulations of star formation. In all cases, gas plays a critical role in the formation of hierarchies.

Sequential disk instability (DI+DI)
While a close binary formed by DI migrates inward, another companion can be formed by DI on its periphery (Tokovinin & Moe 2020). The DI+DI process can create compact triple systems and 3+1 quadruples. In this scenario, hierarchies are assembled inside-out, starting from the innermost and closest binary. The accreted mass is preferentially retained by the least massive outer component because it moves on a wider orbit around the center of mass, but when the components grow and become comparable, the accreted mass is distributed evenly, leading to the formation of twins with q ≈ 1. This consideration applies to both inner and outer subsystems. The DI+DI scenario makes a strong prediction regarding the outer mass ratio: it cannot exceed one. This limit is attained when most mass is accreted after formation of the triple. In such case, the resulting system is a double twin where both inner and outer mass ratios are close to one and the masses of the components are distributed approximately in the 1:1:2 proportion.
The outer, most massive component in a double twin outshines the fainter inner binary and reduces the chance of its discovery. Several double twins were discovered only recently using advanced observational techniques (Tokovinin 2018a). Double twins are still relatively rare compared to other known hierarchies (note the points near (1,1) in Figure 5). A case where the outer companion grows to become the most massive star in a triple system (but not yet a double twin) is illustrated in Figure 8. However, in the majority of simulated DI+DI triples the outer companion forms later and remains the least massive star of the system. The distributions of periods and mass ratios of simulated triples are plotted in Tokovinin & Moe (2020). Tobin et al. (2016) found a triple protostar where a low-mass tertiary component is apparently forming in the circumbinary disk.
Growth of the outer companion can be stopped or slowed down when yet another, fourth star forms at a larger separation, thus converting a nascent triple into a 3+1 planetary-type quadruple. Indeed, in the known planetary-type quadruples the tertiary and quaternary components are usually less massive than the central pair.
The DI+DI scenario implies a planar or quasi-planar orbital architecture if the angular momentum of the accreted gas maintains approximately constant direction during the mass assembly period. Otherwise, accretion of misaligned gas on to a triple or quadruple modifies the orientation of its orbits and destroys the coplanarity, partially or totally. Yet another consequence of the accretion is the inward migration of the orbits. A newly formed low-mass tertiary component overtakes the accretion flow, previously directed to the inner binary, and migrates faster. Later, if the masses become comparable, the gas is distributed more evenly and the inner binary also migrates. The inner and outer orbits, evolving jointly to closer separations, interact dynamically and can be trapped in a mean motion resonance (MMR), at least temporarily (Tremaine 2020). There is some analogy with the migration in multi-planet systems that also produces sometimes resonant or quasi-resonant architectures. On the other hand, if the tertiary component migrates faster than the inner binary, the system can become dynamically unstable and will be disrupted, ejecting some stars.  (CF+CF, DI+CF). A close binary can form by either core fragmentation or disk instability if it accretes enough gas and migrates. Then another component formed at a large distance can be captured, producing a triple system in a DI+CF or CF+CF scenario, as found in the simulations by Lee et al. (2019), Kuffmeier et al. (2019), and Rohde et al. (2021). A variant of this scenario is independent formation of two pairs by either DI or CF mechanisms and their encounter. If all gas is exhausted or removed by the time of the encounter, the dynamical interaction between two binaries can leave behind a triple system and an ejected star. If, on the other hand, the interaction is moderated by gas, both binaries can survive with modified orbits (Ryu et al. 2017).
Elongated filamentary structures produced by gravitational collapse tend to accumulate mass near their ends by means of the so-called edge collapse. Yuan et al. (2020) found that the structure of the 25-pc long filament S228 matches the edge collapse predictions. This mechanism can create two adjacent star clusters or groups of stars of comparable masses. Edge collapse is a consequence of the gravitational focusing that attracts gas to the opposite ends of elongated structures and drives the longitudinal gas motion. This mechanism should also work at smaller scales, producing wide binaries or multiples from elongated cores. In a series of papers, Bonnell & Bastien (1993) explored binary formation from elongated cores and compared predictions of their early simulations with observations. Sadavoy & Stahler (2017) found that the orientation of wide pairs of protostars in elongated cores correlates with the core axis at separations exceeding ∼500 au, while closer pairs are oriented randomly with respect to the core. These observations hint that protostars likely form near the ends of elongated structures and fall to the center of mass.
In the CF+CF scenario, the inner and outer subsystems are formed independently of each other, therefore their orbits can be oriented randomly and the component's masses can span a wide range. Typical triple systems with wide outer orbits qualitatively match these predictions. However, in the hierarchical collapse each structure accretes gas on a time scale longer than its freefall time, so a hierarchy formed by the CF+CF scenario continues to evolve and migrate. As in the DI+DI scenario, accretion of a substantial mass shrinks the orbits and equalizes the masses. Moreover, dissipative interaction with gas reduces the eccentricities and increases the orbit alignment. An initially wide and misaligned 2+2 quadruple can end up as a compact and aligned system of four equal-mass stars. This seems to be the only way to explain the origin of compact 2+2 quadruples with outer separations much less than the opacity limit for fragmentation.
One can envision a cascade of fragmentations where most of the angular momentum remains in the mutual orbit of two fragments which collapse and further fragment, forming a quadruple system of 2+2 hierarchy in the outside-in fragmentation sequence (Bodenheimer 1978). For efficient storage of the angular momentum in the fragments, their masses should be similar. Although such 2+2 architecture is typical of many multiple systems consisting of four stars of comparable masses, most modern hydrodynamical simulations do not favor the cascade hierarchical fragmentation scenario, featuring instead chaotic gas motions and sequential (rather than simultaneous) formation of protostars. In these simulations, hydrodynamic and magnetic transport of the angular momentum is sufficient to enable accretion and buildup of stellar masses. However, recent simulations of a solar-mass core fragmentation by Rohde et al. (2021) yield some triples with comparable (small) masses and aligned orbits resembling the expected products of a fragmentation cascade.

Collisions
Mutual interactions between nascent protostars play an important role in dense clusters, leading to captures, exchanges, disruptions, etc., as demonstrated in the simulations of Bate (2019). However, mutual interactions also can take place in a low-density star-forming regions because they are highly structured and the local density in the clumps can be high even when the average density is low. However, unlike cluster stars, members of small groups interact with each other for a short time, only during their first encounter (Delgado-Donate et al. 2003).
Interactions between gravitating point masses without gas ('dry') are discussed in the following subsection. In the star-forming environments, gas is dynamically important. Collisions at large spatial scales that compress the interstellar medium and trigger star formation (Vázquez-Semadeni et al. 2019) can be relevant at smaller scales as well. A collision between two starless cores creates a shock front and can lead to the collapse of each core into a binary, forming a 2+2 quadruple system in a single event. This mechanism has been proposed and simulated by Whitworth (2001).
An encounter of two accreting protostars can produce a multiple system in several ways. The encounter perturbs the gas envelopes and produces a burst of accretion on to the stars that can lead to the formation of inner subsystems by the DI mechanism (enhancement of accretion by a fly-by or collision is well known in the context of galaxy evolution). At the same time, the kinetic energy of the relative motion is dissipated and the two initially unbound protostars become a wide bound pair. The result can be a 2+2 quadruple, as in the collision between starless cores, or just a triple system. Simultaneous presence of subsystems in both components of wide binaries, noted above, suggests that the subsystems originated in a common event, which could be a collision of protostars surrounded by the gas envelopes or a collision of starless cores. While in the DI+DI or CF+CF scenarios the hierarchies are built from inside-out (inner pairs form first), the collision-formed hierarchies are born in one event.
Along with collisions between cores or protostars, one might envision a collision between a relatively wide binary and a gas cloud. An episode of enhanced accretion promotes buildup of circumstellar disks and can create subsystems around each binary component by the DI mechanism, converting the binary into a 2+2 quadruple. In this CF+DI2 (late disk instability) scenario, the hierarchy is assembled outside-in, starting from the outermost pair and adding the inner subsystems later. With a limited gas supply, the inner companions will not accrete much and the inner mass ratios will remain small. The angular momentum of the gas in the accretion burst is not necesarily aligned with the outer binary, so the two inner subsystems can be mutually aligned but not coplanar with the outer orbit. This hypothetic scenario matches properties of Castor-type quadruples, where the inner subsystems have small mass ratios, while the two outer primaries have comparable masses. A variant of this scenario where the gas cloud colliding with the wide binary also contains a protostar could produce a sextuple system.

Dynamical Interactions
To the first approximation, stellar hierarchies can be modeled by gravitationally interacting point masses (an N-body system). Their dynamics is well understood, despite the lack of analytic theory for systems of three or more gravitating points. The motion can be studied by means of the perturbation approach that represents the inner and outer systems by Keplerian two-body orbits with slowly changing (osculating) elements. This approximation works well for systems with a strong hierarchy, i.e. with a large ratio of periods or separations. Otherwise, the motion can be explored by direct numerical integration. A recent review of triple-star dynamics is published by Docobo et al. (2021), while the book by Valtonen & Karttunen (2006) is recommended for a deeper study of the three-body problem.
When the separations between three masses are com-parable, the triple system is non-hierarchical and dynamically unstable; a regular quasi-Keplerian motion is replaced by the chaotic 'interplay'. Close triple approaches or 'scrambles' lead to ejections of one star on a wide orbit, and during ejections the system temporarily appears as hierarchical, until the next interplay episode (Anosova 1986;Manwadkar et al. 2020). Eventually, one star (usually the least massive one) is ejected from the system, while the remaining binary shrinks. Assuming that the interplay erases memory of the initial conditions, Stone & Leigh (2019) derived the eccentricity distribution of binaries that remain after decay of unstable triples.
The distinction between non-hierarchical (unstable) and hierarchical (stable) triple systems has been studied by many authors. Several criteria of dynamical stability of triple systems based on numerical simulations have been proposed in the literature. For example, the criterion of Mardling & Aarseth (2001) for coplanar orbits can be recast as P out /P in > 4.7(1 − e out ) −1.8 (1 + e out ) 0.6 (1 + q out ) 0.1 , (1) where e out is the eccentricity of the outer orbit and q out is the ratio of the distant-companion mass to the combined mass of the inner binary. The limiting period ratio of 4.7 corresponds to the solid line in Figure 3. There exist multiple systems with period ratios close to this limit (e.g. LHS 10170 and ζ Aqr presented below), hence with strong dynamical interaction between the inner and outer orbits.
Dynamically unstable multiple systems with comparable separations, so-called trapezia, are short-lived, making their discovery in mature stellar populations unlikely. The best chance to find trapezia is by studying the premain sequence (PMS) stars; indeed, several interesting young trapezia candidates are known. Dynamical decay of young unstable multiples leads to the ejection of components (single stars or tight binaries) with a high speed, on the order of the orbital speed at the moment of close interaction. Runaway stars that move away from young star-formation regions witness the dynamical decay of young multiples, while their velocities provide an order of magnitude estimate of the size of those systems at the moment of disruption. Existence of hierarchies with moderate period ratios is an indirect evidence that some unstable young systems have decayed, while the stable and marginally-stable ones survived.
Triple systems that experienced a chaotic dynamical evolution bear characteristic imprints of this process: their orbits are usually misaligned and have large eccentricities, while the period ratios are moderate (not too far from the stability limit). These features are inferred from the simulations of cluster decay (Sterzik & Tokovinin 2002) and from the numerical scattering experiments (Antognini & Thompson 2016). The eccentricity distribution is thermal, f (e) = 2e, or even steeper than thermal (Stone & Leigh 2019).
Stars typically form in small groups with comparable separations. These small-N clusters should evolve dynamically, ejecting some members. When such cluster is immersed in a massive gas cocoon, the ejected star can be pulled back by the additional gravity of the gas and return, instead of being lost. Reipurth & Mikkola (2012) suggested such ejections as a mechanism for forming very wide binaries and called it unfolding. In this scenario, each wide binary should contain an inner close pair (unless that pair migrated strongly and merged), explaining the observed correlation between wide multiples and triples. The orbit of the unfolded tertiary is necessarily very eccentric because its angular momentum is comparable to the angular momentum of the initial unstable triple which, according to these authors, can have a size of a few hundred au. If the wide companion repeatedly returns to the central gas cloud, its orbit might be circularized and the triple will no longer be wide. The scenario proposed by Reipurth & Mikkola is elegant, but some predictions of the unfolding mechanism do not match observations. In populations younger than a few Myr, wide binaries are already abundant, whereas unfolding postulates their delayed formation. The outer eccentricities of wide triples are not always large, signaling that unfolding is unlikely to be the dominant channel of wide binary (or triple) formation; there are other, less exotic ways (Tokovinin 2017a).
When the distance between binary components becomes small, comparable to the stellar radii, the stars can no longer be treated as point masses. Tidal forces cause precession of the inner orbit and its circularization. A relatively wide triple with large mutual inclination, 39 • < Φ < 141 • , experiences Lidov-Kozai (LK) cycles that modulate both Φ and e in (see Naoz 2016, for a review). The increased inner eccentricity may 'switch on' tidal interaction in the inner orbit near its periastron. As a result, the LK cycles break up and the inner orbit shrinks and circularizes. This mechanism of close-binary formation is called Kozai cycles with tidal friction, KCTF (Eggleton & Kisseleva-Eggleton 2006;Eggleton 2006). It can be viewed as a kind of migration where the angular momentum of the inner subsystem is transferred to the orbit of the tertiary companion. Formation of close binaries by the KCTF mechanism matches their preference to be inner subsystems in multiples. This is particularly true at periods shorter than a few days: such binaries cannot form very early because the radii of PMS stars are large and the binary would have merged early on. The strong tendency of close binaries to be members of hierarchical systems (Tokovinin et al. 2006) finds its explanation in the KCTF mechanism, at least partially. Fabrycky & Tremaine (2007) explored statistics of close binaries formed by this channel, although Moe & Kratter (2018) argue that it cannot be the dominant mechanism of close-binary formation.
Dynamical interactions in mature hierarchical systems are astrophysically important and can produce some exotic objects (Hamers 2020;Toonen et al. 2020). However, they are outside the scope of this review.

FAMILIES OF HIERARCHICAL SYSTEMS
In this Section, I explore the relation between hypothetical formation scenarios of hierarchies and the properties of real systems, trying to develop their 'genetic' classification. The knowledge of formation processes is still qualitative, their predictions are tentative, and a given system can be associated with several formation scenarios. Hopefully, this attempt will stimulate further theoretical and observational work in this area. Figure 3 reproduces the P in , P out diagram and marks the location of groups discussed below. However, the classification is based not only on the periods, but on other parameters (mutual orbit orientation, eccentricities, mass ratios) as well, so different groups overlap in this diagram. The upper left corner is occupied by the hierarchies with short inner and long outer periods, P out /P in > 100. These hierarchies can be called typical. Owing to the large period ratios, interaction between their inner and outer orbits is negligible. Hierarchies with smaller period ratios are more interesting, and the following discussion is focused mostly on these systems. To help visualize the forthcoming text, Figure 9 shows schematically various families of triple and quadruple systems. Table 2 gives the main characteristics of several selected hierarchical systems discussed in the following subsections. The information is extracted from the MSC. The systems are identifed by the WDS codes and common names. Then the total number of known components N and the hierarchy type H are given. The last column gives the orbital periods in the bracket notation, like P out (P in1 ; P in2 ) for a 2+2 quadruple.

Planar Systems
There is a general trend of increasing orbital alignment with decreasing outer separation. Visual triples with outer separations less than ∼50 au tend to have aligned orbits: their average mutual inclination Φ is about 20 • (Tokovinin 2017b). Moreover, in the corotating (hence likely aligned) systems the inner eccentricities are, on average, smaller than in the counter-rotating systems. Planar architecture and moderate eccentricities suggest that these systems were possibly formed by the DI+DI scenario. The inner binary forms first, grows in mass and evolves to closer separations, while the outer companion(s) form later. Further accretion of substantial mass by such a triple should produce a double twin with both inner and outer mass ratios close to one. Several such visual triples were identified in Tokovinin (2018a) and called 'dancing twins'. Moderate period ratios in these systems imply strong dynamical interaction between the inner and outer orbits.
One of the most remarkable dancing twins is the lowmass triple system LHS 1070 (00247−2653). The main star A has a mass of 0.12 M ⊙ , and the inner subsystem B,C is composed of two stars of 0.075 M ⊙ each, on the borderline between hydrogen-burning stars and brown dwarfs. Both inner and outer mass ratios are close to one. The inner period is 18.2 yr and the outer period (still preliminary) is 99 yr according to Xia et al. (2019). Both inner and outer orbits are nearly circular and coplanar (Φ = 1.7 • ). The period ratio of 5.44 is close to the limit of dynamical stability; however, numerical integration performed by Xia et al. confirms stability on a ∼1 Gyr time scale. Tokovinin (2018a) found earlier the period ratio of 4.5, suggestive of a 9:2 MMR. A longer time coverage is needed to better constrain the outer period and thus to confirm or refute the MMR hypothesis. Interaction between the orbits causes measurable deviations from the Keplerian motion in the inner pair of LHS 1070. The structure of LHS 1070 resembles low-mass triples formed by fragmentation of isolated cores in the simulations of Rohde et al. (2021).
Double twins are difficult to discover because the bright tertiary outshines the faint inner pair. This means than   only a handful of them are studied so far. Newly discovered double twins still lack the coverage of the outer orbits. Most known planar triples with outer separations below ∼50 au are not double twins,nhowever. Their mutual inclinations are measured when both inner and outer orbits are resolved, either directly or through the astrometric wobble. A special case where the mutual inclination is deduced from dynamical modeling of compact triples is discussed below in section 4.2.
As noted above, the growth of the tertiary companion is stopped when a more distant fourth companion forms and overtakes the accreted gas flow. In such case, the result is a quasi-planar 'planetary' hierarchy of 3+1 type. An example of such system is HD 91962 (10370−0850) studied by Tokovinin et al. (2015). It is called planetary because three small companions of K, M spectral types revolve around the central most massive star of G0V type and also because their orbits have moderate eccentricities (0.30, 0.125, 0.135, counting outside-in). The mutual inclination between the outer and middle orbits is 11 • , and the less certain mutual inclination between the innermost and middle orbits, determined later by the astrometric wobble, is 32 • ± 12 • . The mutual inclinations and the eccentricities in HD 91962 are larger than in the solar system, but smaller than in typical multiples. The ratio of the middle and inner periods, 18.97 ± 0.06, is close to an integer number, but a 1:19 MMR is very weak, hence unlikely.
The 3+1 quadruple HIP 78842 (16057−3252) has the inner triple with periods of 10.5 and 131 yr (period ratio 12.6), again with nearly coplanar and circular orbits and with a mutual inclination of 12 • (Tokovinin 2018a). The inner pair Ba,Bb is a twin with masses of 0.75 M ⊙ each, but the mass of the main star A is only 0.96 M ⊙ , so this is not a double twin. The reason might be the existence of a fourth companion C at 9.3 ′′ separation (the estimated period of AB,C is ∼4 kyr) with a mass of 0.64 M ⊙ . The parameters of the outer orbit cannot be determined owing to its long period. However, the small eccentricity of the intermediate 131-yr orbit, e = 0.03, implies the absence of the LK cycles, so the outer mutual inclination should be less than 39 • . Compared to HD 91962, this planetary-type quadruple is much wider, and all its subsystems are resolved.
On the other hand, a compact 3+1 quadruple system composed of three very similar K7V dwarfs and a lowermass outer companion, HIP 41431 (08270+2157), was discovered and studied by Borkovits et al. (2019). The orbital periods are 2.9 days, 59 days, and 3.9 yr. The mutual inclination between the two inner orbits is 2.2 • , and the outer orbit is inclined by ∼21 • . The inner pair is eclipsing. Accurate eclipse timing reveals strong dynamical interaction between all three orbits. The intermediate 59-day orbit precesses under the influence of the outer 4yr companion. The system is not resolved visually, the mutual inclinations are determined by dynamical modeling.
Summarizing, the family of planar hierarchies is distinguished by their approximately aligned orbits (Φ < 40 • ), moderate period ratios P out /P in < 50, and moderate eccentricities (e < 0.5). The outer separations are also modest, typically < 50 au. These features suggest formation by the DI+DI scenario. Resonances predicted by this scenario are not yet proven to exist in the real systems. Although Quirrenbach et al. (2019) discovered that the giant star ν Oph is orbited by two brown dwarfs in a MMR (periods 503 and 1385 days, ratio 1:6), the small mass ratios are more typical for planetary systems. The planar family includes double twins and 3+1 planetary quadruples, but the majority of its members are more common triples. A subset of compact planar hierarchies with short outer periods is discussed in the following section.

Compact Hierarchies
Looking at Figure 3, we note a relatively small number of points with P out < 10 3 days in the lower left corner of the diagram. Rareness of such hierarchies is even more evident in the volume-limited sample of solartype stars (Tokovinin 2014). I call these systems compact; they belong to the planar family, as the orbits are normally aligned. An example of a compact planar system is HIP 41431 (08270+2157) presented above (Borkovits et al. 2019).
Close tertiary companions to spectroscopic binaries are readily detectable by variation of the systemic velocity. Despite this, such cases are rare owing to the intrinsically low fraction of compact hierarchies. Among several known examples of compact (but not eclipsing) spectroscopic triples, I can mention HIP 7601 (01379−8259, periods 1.75 yr and 19.4 days) and HIP 15988 (2366−0034, periods 271 and 5.9 days).
Tertiary companions to close binaries can be discovered by the astrometric acceleration they produce. One of the most sensitive probes is the so-called proper motion anomaly (PMA) -the difference between the long-term proper motion deduced from the positions measured by the Gaia and Hipparcos space missions, and the accurate short-term motion measured by Gaia (Brandt 2018;Kervella et al. 2019). Examination of the PMA of close binaries in the 67-pc sample of solar-type stars revealed a few dozen new triples that were added to the MSC. However, the periods of the astrometric tertiary companions are unknown, and only a fraction of those triples could be compact.
The best way to detect and study compact hierarchies is based on the analysis of eclipses in the inner subsystems. A large number of such hierarchies were discovered by eclipse timing variations in the sample of Kepler eclipsing binaries by Borkovits et al. (2016). They are indeed close to coplanarity, with a few exceptions. The Kepler compact hierarchies are mostly beyond 200 pc from the Sun, and for this reason they are not plotted in Figure 3. This technique has led to the discovery of ultra-planar compact hierarchies. For example, Borkovits et al. (2020) found a triple system TIC 209409435 with periods of 121.9 and 5.7 days where the mutual inclination is constrained to be Φ < 0.25 • . Several other compact ultra-planar systems were discovered recently by this method. One of the shortest outer periods, 33.9 days, is found for KIC 5897826 (KOI-126, 19499+4107); the inner period is 1.8 days (Borkovits et al. 2016). The young triply eclipsing system HD 144548 (16073−2204) in the Uppper Scorpius association has very similar periods of 33.9 and 1.63 days (Alonso et al. 2015). Borkovits et al. (2020) give a list of 17 compact triples contaning eclipsing binaries; their outer periods range from 33 to 1100 days.
The most compact among known 2+2 quadruples is VW LMi (11029+3035, HD 95660) with an outer period of only 355 days (Pribulla et al. 2020). The inner periods are 7.93 and 0.48 days. All four stars have similar masses of ∼1 M ⊙ and are tightly packed in an outer orbit of only 1.6 au size. The coplanarity of orbits in this unresolved system is inferred from the dynamical analysis. The doubly eclipsing system V994 Her (18278+2442) is similar to VW LMi, but its outer period is longer, 2.9 yr; the inner periods are 2.08 and 1.42 days, and all four masses are comparable, from 2 to 3 M ⊙ . The compact 2+2 quadruple with the second-longest outer period of 432 days is TIC 454140642 (04191+0054), a doubly-eclipsing star discovered by Kostov et al. (2021). Its inner periods are 13.6 and 10.4 days, and the masses of all four stars are remarkably similar, from 1.1 to 1.2 M ⊙ . The system is not young (estimated age 1.9 Gyr).
Modern large-scale photometric surveys such as OGLE revealed a vast number of eclipsing binaries. Some of those turned out to be doubly eclipsing, being in fact 2+2 quadruples or higher-order hierarchies (Zasche et al. 2019). Double eclipses are more likely when the subsystems are mutually aligned, which is expected in compact hierarchies like VW LMi. Indeed, in some doublyeclipsing quadruples it was possible to detect motion in the outer orbit by the anti-correlated variations of the eclipse time. For example, OGLE BLG-ECL-145467 (17521−2920) has the outer period of 4.2 yr and the inner periods of 4.9 and 3.9 days. One of the most intriguing findings of this work is the hint on mutual resonances between the periods of the inner subsystems: the 3:2 ratio of inner periods is frequent and the 1:2 ratio seems to be rare. This can be explained by the interaction of both subsystems with their outer orbit, which is possible only in compact configurations. Yet another necessary condition for such resonances is migration of the orbits (Tremaine 2020).
Proposing a formation scenario of compact 2+2 quadruples of VW LMi type is a challenge. Some of their properties (compactness, coplanarity, similar masses) point to the accretion-driven migration. However, the DI+DI scenario for compact triples, where companions form and migrate sequentially, does not work for 2+2 quadruples. It seems that relatively wide 2+2 quadruples formed by one of the mechanisms outlined above (CF+CF, CF+DI2, collision) evolved to their present compact configurations through accretion-driven migration. This scenario predicts that all compact 2+2 quadruples must have large inner mass ratios, as observed in VW LMi, TIC 454140642, V994 Her, and OGLE BLG-ECL-145467. Coplanarity is also predicted because the gas accreted by the inner subsystems is aligned with the outer orbit.

Castor-type Quadruples
Castor (α Gem, 07346+3153) is a system containing 6 or 7 stars. Its central binary A,B consists of two similar bright early-A type stars and it has been measured since 1778; however, the orbital period of 460 yr is not yet fully covered. Both components of this visual binary contain spectroscopic subsystems with periods of 9.2 and 2.9 days and minimum masses of 0.2 and 0.36 M ⊙ , respectively. Small inner mass ratios distinguish Castor from the majority of 2+2 quadruples where all four com-ponents typically have similar masses.
The architecture of the sextuple system 88 Tau (04357+1010) is very similar to that of Castor, but its inner quadruple has a period of 18 yr and its inner spectroscopic subsystems have periods of 7.9 and 3.6 days and relatively large mass ratios (spectral types from G2V to A6V). The interferometric study by Lane et al. (2007) revealed that the inner subsystems in this quadruple are not aligned with the 18-yr outer orbit and the mutual inclinations Φ are 143 • and 82 • . In this work, the inner subsystems were not resolved directly, and the orientation of their orbits was determined from the wobble in the motion of the 18-yr pair.
In hierarchies assembled inside-out, the inner subsystems typically have large mass ratios resulting from the accretion-driven migration. The architecture of Castor is quite distinct and suggests that it could form in reverse order by the CF+DI2 scenario, where an existing binary A,B experienced an episode of late accretion burst that formed the inner subsystems by disk instability. The accretion episode could be caused by an encounter with another protostar surrounded by its own disk that also became unstable and formed a close pair that remained bound to the inner binary A,B, making it a sextuple system. Although the proposed scenario looks exotic, it can explain the unusual architecture of Castor-type hierarchies.
Recently, another sextuple system with an architecture resembling Castor has been discovered by Powell et al. (2021). This is the triply-eclipsing star TIC 168789840 (04141−3155). It contains three eclipsing subsystems with periods of 1.3, 1.5, and 8 days; all six sets of eclipses (3 primary and 3 secondary) are observed. The two closest binaries are paired in an orbit with a period of a few years, while the 8-day pair belongs to the outermost component on a wide orbit with a period of ∼2 kyr. The edge-on orientation of the three inner orbits could be a mere coincidence. More likely, however, the inner orbits in this triply-eclipsing hierarchy are alined (the same probabilistic argument was advanced for the doubly-eclipsing systems). Another particularity of TIC 168789840 consists in the remarkable similarity of the component's masses and inner mass ratios. The primary components of all eclipsing subsystems have masses about 1.4 M ⊙ , and their secondaries are all about 0.6 M ⊙ . It looks as though all three inner subsystems were made at the same factory according to the same 'blueprint'. A potential formation scenario of this system is CF+DI2, where a young triple composed of similar-mass stars encounters a gas cloud, and the accretion burst forms secondary subsystems with similar properties around each component of the original triple.

Quadruples of 2+2 Hierarchy
Quadruple systems of 2+2 hierarchy (two close binaries on a wide orbit around each other) are rather common. Their fraction among solar-type stars is 4%, compared to the 1% fraction of 3+1 quadruples (Tokovinin 2014). Compact 2+2 quadruples like VW LMi and the Castortype quadruples discussed above belong to this group, which contains hierarchies with a wide range of periods and mass ratios.
The classical visual quadruple system ǫ Lyr (18443+3940, HR 7051-7054) consists of four simi-lar stars of spectral types from F1V to A4V, arranged in two pairs with periods of 1800 and 700 yr (these orbits are poorly constrained) at a large projected separation of 11.5×10 4 au. The similarity of all four masses and of the two inner periods is the characteristic of this family, called ǫ Lyr type for brevity (Tokovinin 2008). Other, more compact members of this family are known, e.g. the visual quadruple FIN 332 (18455+0530, HR 7048), where all four components are similar (spectral types from A0V to A1V), the outer projected separation is 432 au, and the two inner orbits have periods of 40 and 48 yr; both inner orbits also have large eccentricities of 0.82 and 0.84 (Tokovinin 2020b). Yet another quite typical example is HIP 43947 (08571−2951), where the two pairs of similar solar-type stars, at 100 au projected separation, have separations of 3-4 au each and estimated periods of ∼100 yr.
Maximum periods of inner subsystems in a 2+2 quadruple are limited by the dynamical stability criterion, hence some correlation between the two inner periods is expected. However, the statistical analysis of such quadruples in (Tokovinin 2008) shows that the similarity of the inner periods is stronger than required for the stability alone, hinting that the inner pairs have not formed independently of each other. Likewise, the two inner mass ratios in 2+2 quadruples appear to be correlated.
The ǫ Lyr type 2+2 quadruples resemble products of cascade hierarchical fragmentation driven by the conservation of the angular momentum in the orbits of the subsystems, envisioned by Bodenheimer (1978). The relatively wide outer and inner separations (above the opacity limit for fragmentation) and similar masses of the components match the predictions of this scenario. Nevertheless, the approximate coplanarity of orbits implicit for this mechanism does not hold. Although the sample of resolved visual 2+2 quadruples with known sense of revolution is rather small, the numbers of apparently co-rotating and counter-rotating inner pairs are similar, demonstrating the lack of mutual alignment between the inner subsystems. However, the outer separations in this sample are larger than ∼1000 au, and such wide triples also have misaligned orbits. In the more compact quadruple FIN 332 the two inner orbits could be mutually aligned, but the outer pair rotates in the opposite sense, excluding alignment of the whole system. The ǫ Lyr system is wide, unlike the compact 2+2 quadruples presented above, but both groups have one common feature, namely the similarity of masses. However, 2+2 quadruples with unequal masses are also quite common.
A relevant example is HIP 12548 (02415−7128), a classical visual binary with an outer period of 100 yr and a semimajor axis of 0.6 ′′ . Its two solar-type components A and B, resolved visually or interferometrically, are never separated by the seeinglimited spectroscopy. Existence of a spectroscopic subsystem was suspected from the occasional doubling of the spectral lines. Regular monitoring revealed that each component is in fact a spectroscopic binary. Their periods are relatively long, ∼2010 and 110 days, and the inner mass ratios are small, ∼0.3 (Tokovinin 2021, in preparation). Such quadruples are very difficult to discover. Other visual binaries where each component has a spectroscopic subsystem are known, e.g. HIP 41171 (08240−1548) with an outer period of ∼900 yr and the inner periods of 25.4 and ∼960 days. The masses of all components of HIP 41171 are comparable and in some orbital phases the spectrum contains four distinct systems of spectral lines, although most of the time the lines are blended.
Summarizing, 2+2 quadruples have a wide range of parameters, from wide configurations of ǫ Lyr type consisting of similar-mass stars to the most compact known quadruple VW LMi, also with similar masses. Most 2+2 quadruples are situated in the middle of this range, and their masses can be quite diverse. Such quadruples are difficult to discover, especially when their secondary components are faint and do not have signatures in the spectrum. One notable exception are doubly-eclipsing systems. The diversity of 2+2 quadruples suggests that they could be formed by a variety of mechanisms.

Misaligned Hierarchies
The trend to mutual orbit alignment weakens with increasing outer separation. However, misaligned hierarchies exist even at small separations. Several such systems are well documented.
A typical example of a misaligned triple is ζ Aqr (22288−0001, HR 8559+8558) (Tokovinin 2016). This is a bright visual binary known since 1777. A wobble in its motion on the 540-yr outer orbit was detected in 1955, and the inner subsystem Aa,Ab was directly resolved in 2009; its 26-yr orbit is now well defined. The most remarkable feature of this triple is the large inner eccentricity, e in = 0.87, and the large mutual inclination of 140 • : the inner and outer pairs are definitely counter-rotating. The period ratio of 20.8±0.6 is small and, considering the outer eccentricity of 0.42, this triple is located near the limit of dynamical stability (eq. 1). The visual components A and B have similar masses of 1.4 M ⊙ (an outer twin), but the inner secondary is substantially less massive, 0.6 M ⊙ .
Although the modest period ratio and the outer semimajor axis of 100 au roughly match parameters of planar visual triples, other characteristics of ζ Aqr are quite distinct. The large inner eccentricity and the counterrotation resemble products of dynamical interactions. Usually, the least massive star is ejected from a dynamically unstable triple, but we may envision a reverse process, namely a capture. A low-mass star can encounter the existing binary, interact with it dynamically, and become captured. The mechanics of point masses is timereversible and allows capture (opposite of ejection), given the appropriate initial conditions. A capture is probably more likely in an encounter of two binaries because in this case one star can be ejected, removing the excess of energy and angular momentum that is needed for a capture (in a triple-star capture, the energy and momentum are absorbed by the outer binary). Equal masses of the main components A and B speak for their common origin, while the low-mass inner companion could be formed independently and captured later.
This example is by no means unique. Another very similar system is 02460−0457 (HD 17251). This is a 1 ′′ long-period visual binary composed of two solar-type stars (masses 1.4 and 1.0 M ⊙ ). A wobble has been suspected in its motion, and the faint inner companion has been eventually directly detected in 2016. Its orbit, still

kyr
A,B Figure 10. Hierarchical structure of two wide nearby quintuple systems κ Tuc and ξ Sco. Periods of the subsystems (green circles) and masses of the stars (pink circles) are indicated.
preliminary, has a period of 38 yr, an eccentricity of 0.6, and the rotation sense opposite to the outer binary; the mass of the inner companion is only 0.44 M ⊙ . The resemblance to ζ Aqr is quite obvious. In yet another nearby triple 20396+0458 (HIP 101955), the mutual inclination is 65 • , so the orbits are closer to being perpendicular rather than coplanar. Not surprisingly, the inner eccentricity is substantial, 0.69, and the LK cycles are certainly going on. The periods are 38.7 and 2.51 yr (ratio 15.4), and the masses of all three stars are comparable, between 0.6 and 0.7 M ⊙ .
The statistical study of mutual orbit orientation in visual triples (Tokovinin 2017b) hints at a decreased alignment with increasing stellar mass. Massive stars form in dense environments and are more prone to dynamical interactions with their cluster neighbors and with their own companions. Moreover, assembly of massive stars implies strong accretion which favors companion formation by DI and their fast inward migration. Migrating outer companions de-stabilize the inner subsystems, leading to ejections (runaway stars), and leave behind hierarchies with misaligned and eccentric orbits. The most compact known triple, λ Tau (04007+1229, HR 1239) has a B3V primary component; its periods are 33.07 and 3.953 days (ratio 8.37, near the limit of dynamical stability). The inner pair is eclipsing, and the inclination of the 33-day pair is estimated at 65 • , suggesting misaligned orbits despite the short outer period. A better documented example of a misaligned massive triple is σ Ori (05387−0236, HR 1931), an O9.5V member of the Orion Trapezium cluster. Schaefer et al. (2016) resolved the inner 143-day spectroscopic subsystem with an eccentric (e = 0.77) orbit and determined that its inclination to the outer visual orbit of 157 yr period is either 130 • or 114 • , so the subsystems counter-rotate.

High-Multiplicity Systems
Hierarchies containing 5, 6, or 7 stars can be collectively called large-N systems. Architecture of one such system, 65 UMa, is illustrated in Figure 2. Castor (α Gem, 07346+3153) contains 6 or 7 components. Its central quadruple AB, discussed above, is paired with the twin eclipsing binary C = YY Gem at projected separation of 1100 au. The period of this pair Ca,Cb is 0.81 days. Wolf et al. (2018) discovered eclipse time variation with a 54 yr period that could be caused by a substellar companion of 0.05 M ⊙ associated with C. If this companion is real, Castor is a septuple system. It has only three hierarchical levels, allowing a maximum of 2 3 = 8 components. In this system, almost all available slots are filled, which is atypical.
A large-N hierarchy can be decomposed into simpler constituents. For example, the Castor system is a combination of the 2+2 quadruple and a triple. The large number of components is normally associated with the large range of separations: the outer systems in large-N hierarchies are wide, and the inner subsystems are close, as in Castor.
Two nearby (within 30 pc) quintuples that defy the general trend and contain only relatively wide subsystems, namely κ Tuc (01158−6853) and ξ Sco (16044−1122), are discussed by Tokovinin (2020a) and illustrated in Figure 10. Their outer periods are ∼300 kyr (projected separations ∼8 kau), and all inner periods are longer than 20 yr. Both hierarchies are ∼2 Gyr old and, naturally, are not associated with any young kinematic group. The wide outer subsystems indicate an absence of strong dynamical interactions with neighbors, suggesting formation in a low-density environment like Taurus. The age of these systems speaks for their stability, ruling out strong internal dynamical interactions. In harmony with the wide outer separations, the orbits of the inner visual subsystems appear to be oriented randomly (some pairs revolve in opposite directions). Masses of the stars in these hierarchies are similar (from 0.8 to 1.5 M ⊙ ) and do not match a random selection from the IMF. The formation scenario of such systems proposed in (Tokovinin 2020a) is a filament fragmentation in relative isolation (see above). The similarity of masses is possibly explained by the prolonged accretion of gas on to the systems; however, for some unknown reason the accretion has not shrunk the inner orbits.
Given that the multiplicity is a strong function of stellar mass, one might expect that large-N hierarchies are rare among low-mass stars; at least, such is the prediction of the IMM statistical model (section 2.4). Yet, low-mass quintuples are found in the immediate vicinity of the Sun, casting doubt on their rareness. GJ 2069 (08316+1924, HIP 41824, parallax 60 mas) contains five similar stars of spectral types from M4V to M3.5V arranged in a 3+2 architecture: two visual binaries at 10.4 ′′ from each other, with one of those containing a 2.8-day spectroscopic and eclipsing subsystem CU Cnc. Even more interesting is the quintuple GJ 644 (16555−0820, HIP 81817, Wolf 630) located at 6 pc from the Sun. Its innermost compact triple composed of M3V dwarfs (periods 1.7 yr and 3 days) is accompanied by an M4V dwarf at 21 ′′ and another M7V dwarf at 232 ′′ , making it a rare 4-tier hierarchy.

Young Hierarchies
In the standard paradigm of low-mass star formation, the envelope is accreted in less than 10 5 yr, so the chances of observing multiple-star formation directly are small. On the other hand, in the hierarchical collapse framework the accretion and mass assembly continue for 1-10 Myr and is still ongoing in the nearby star-forming regions (SFRs) like Taurus-Auriga and Orion. Formation and early evolution of hierarchical systems is therefore directly observable. Study of young hierarchies helps to develop and test the formation scenarios.
The prototype of pre-main sequence (PMS) stars, T Tau (04220+1932), can also be considered as a prototype of typical triple systems. Its heavily obscured southern component S, at 100 au projected separation from the northern star N, is itself a binary with a period of 27 yr and masses of 2.1 and 0.5 M ⊙ . The binary Sa,Sb is accreting from the circumbinary disk and drives an outflow; star N is also accreting. The kinematics of the gas and stars is complex, indicating that this system is misaligned and 'in disarray', using the words of Kasper et al. (2020). Most of the mass is concentrated in the component S hosting the inner subsystem, suggestive of the inside-out formation scenario by the DI+CF or CF+CF mechanisms. The fact that star N also accretes speaks against its dynamical ejection; more likely, it formed in the same cloud, approached S, and was captured. Continued accretion and dynamical evolution will shrink both inner and outer orbits, and in the end T Tau might become a typical triple system with an outer period on the order of a few centuries, its primary component hosting a close subsystem Sa,Sb.
Interesting insights on the multiple-star formation are provided by the PMS triple GW Ori (05291+1152) located in the λ Ori SFR with an age of ∼1 Myr. The inner binary consists of stars with masses of 2.5 and 1.4 M ⊙ on a nearly circular orbit with a period of 241.6 days. The 1.4 M ⊙ tertiary component C has a period of 11.5 yr, and the eccentricity of its orbit is 0.38. The orbits are well-constrained, their mutual inclination is 14 • . This system is actively accreting from the surrounding gas which has a remarkable structure. The dust emission comes mostly from two coplanar outer rings with the radii of 180 and 350 au and from the tilted inner ring of 43 au radius. Interferometric observations in the mm range and high-contrast imaging in the IR revealed the 3D structure of these rings in unprecedented detail ). All rings rotate in the same sense, but the inner ring is inclined to the outer orbit of the triple by 39 • . Hydrodynamic simulations can reproduce the system's geometry and show that the inner ring was torn from the disk by its dynamic interaction with the triple; for the same reason the outer disk is warped. Like T Tau, this object proves that mass assembly can last for ∼1 Myr and the angular momentum of the infalling gas can be misaligned relative to the stellar system.
When several stars form in close proximity, their dynamical interactions are inevitable. Although the overall stellar density in the Taurus-Auriga SFR is low, the association is highly structured and contains compact dense groups of stars. According to Ménard et al. (2020), the triple system UX Tau (04301+1814) presents evidence of a recent fly-by of another star. Its main components A and B are separated by 5.7 ′′ (825 au); B is a 0.1 ′′ binary, A is surrounded by a disk. There is another star C at 2.7 ′′ from A. Comparable separations of AB and AC imply dynamical instability. However, Ménard et al. show that star C is likely a fly-by. It interacted recently with the disk of A, leaving a typical signature in the form of spirals, and captured some gas into its own small disk; a tenuous gas bridge still connects A and C. The authors cite another examples of fly-bys of PMS stars evidenced by the residual gas structures such as bridges, tidal tails, and spirals. A nascent trapezium system has been discussed by Reipurth & Friberg (2021). This is IRAS 05417+0907 in the λ Ori SFR associated with the Herbig-Haro flow HH 175. This object contains six sources at comparable separations on the order of a few thousand au. One of those stars is immersed in a dense envelope and drives the outflow. The authors suggest that this is a binary formed by dynamical interaction in the trapezium and that the associated accretion burst created the outflow; its dynamical age is 6 000 yr. An eccentric binary continues to accrete gas, preferentially near periastron, and evolves to closer separation. Some Herbig-Haro jets (e.g. HH 111) have quasi-periodic knots that could be produced by accretion bursts driven by the inner evolving binary of ∼10 au separation (Reipurth 2000). Apart from the N-body dynamics, binaries driving the jets can result from gasassisted captures, as happens in the cluster simulations.
GG Tau (04325+1732) is a famous PMS hierarchy. Its components A and B, separated by 10 ′′ (1400 au), are close binaries with circumbinary disks, and AB is surrounded by a circum-quadruple disk, earning this object the nickname 'ring world'. Moreover, the component Ab is also a tight 0.03 ′′ pair, making this system a quintuple. The masses of these stars range from 0.2 to 0.65 M ⊙ . A recent paper by Keppler et al. (2020) explores the structure of the disk around Aa,Ab (cavity, filaments, shadows, accretion streams), trying to infer orientations of the unresolved circumstellar disks.
The 2+2 PMS quadruples are represented by the system HD 98800 (11221−2447), also known as TWA 4 because it belongs to the TWA association (age ∼10 Myr). The outer system A,B is a classical visual binary with a period of 200-300 yr and a moderate eccentricity of 0.4; the orbit of A,B is highly inclined. Both stars are also spectroscopic binaries with periods of 262 and 315 days, respectively. Moreover, the double-lined subsystem Ba,Bb was resolved using Keck interferometer, allowing Boden et al. (2005) to determine the mutual inclination and the masses (0.7 and 0.6 M ⊙ ). B is surrounded by a dusty and gaseous circumbinary disk responsible for the large IR excess, while A has no disk. The disk around B, directly imaged by ALMA (Kennedy et al. 2019), is a ring of 3.5 au radius and 2 au width, oriented almost 'face-on' (inclination 154 • ). Hence, this disk is perpendicular to the orbits of Ba,Bb and A,B. The subsystem Aa,Ab also has been recently resolved interferometrically (Zúñiga-Fernández et al. 2021).
A nascent system of three protostars L1448 IRS3B in the Perseus SFR (age less than 150 kyr) has been discovered and studied by Tobin et al. (2016). The inner binary consists of two solar-mass stars A and B at projected separation of 61 au. It is immersed in a massive (0.3 M ⊙ ) gaseous disk with a spiral structure, and another star C is forming in this disk at projected separation of 183 au. Star C has a low mass, but looks brighter than A and B owing to the faster accretion. The authors present this system as a classical case of the DI+DI formation scenario; the massive disk surrounding A and B is gravitationally unstable in the region where C is forming.

SUMMARY AND DISCUSSION
Many observational and theoretical works mentioned in this review are quite recent, evidencing that architecture of stellar hierarchies is presently a very active research topic. The interest is stimulated by the synergy with exoplanets, where hierarchies help to understand the origin and evolution of the angular momen-tum of stars and their systems and inform us on the dynamics of gas from which planets will form. New observational techniques such as interferometry and accurate photometry and astrometry from space also stimulate this research. Our concepts related to formation of stars and their systems evolve in response. These shifting paradigms are summarized below, followed by an outlook of observational and theoretical progress expected in the near-term future.
The main conclusion from this review is the diversity of hierarchical systems and the inferred diversity of their formation scenarios. The three main ingredients of this process are fragmentation, accretion, and dynamics. The scenarios proposed here are based on general considerations (e.g. hierarchical collapse and angular momentum extraction) and recent simulations, but they remain tentative.

Evolution of Concepts 1. Star formation is a process, not an event.
Early works on star formation considered collapse of isolated clouds happening on a relatively short free-fall time scale of ∼ 10 4 yr (Larson 1972). Now we know that assembly of stellar masses, as well as formation and evolution of multiple systems, extends over 1-10 Myr; this time scale is defined by the size of the largest collapsing structures, molecular clouds and clusters (Vázquez-Semadeni et al. 2019). Star formation proceeding during this time is now likened to a 'conveyor belt' (Krumholz & McKee 2020), and this notion helps to model the IMF/SMF (Clark & Whitworth 2021). The architecture of stellar systems is defined by their early evolution as much as by the elementary mechanisms of their formation; in fact, systems formed by different mechanisms such as DI and CF may end up in similar configurations. Complex evolution of stellar systems in a young collapsing cluster is beautifully illustrated by the simulations of Bate (2019).
2. The critical role of gas in the formation and dynamics of stellar systems is incontestable. Accretiondriven migration forms close binaries during the massassembly period (Tokovinin & Moe 2020). Dissipation of kinetic energy in the gas is also relevant for the formation of wide systems by disk-assisted capture. In the past, multiple systems were often treated simply as gravitating point masses. Only at close separations, the point-mass dynamics was modified by considering the finite stellar size and tidal forces (Naoz 2016). At large separations, interactions with neighboring stars destroy wide pairs by the so-called 'dynamical processing' (Duchêne & Kraus 2013). However, the N-body dynamics alone fails to match the architecture of real hierarchical systems, as demonstrated, e.g., by the simulations of Sterzik & Tokovinin (2002). Accreting stars and their systems evolve jointly and self-consistently (Bate 2019). Early works where a 'dynamical processing' of binary population in a cluster was studied by assuming some ab initio binaries with universal properties slowly lose their relevance. This said, dynamical evolution of mature hierarchies (after gas dispersal) leads to several astrophysically important outcomes (Hamers 2020;Toonen et al. 2020), while dynamical interactions during or shortly after formation of hierarchies certainly take place and leave their imprints such as misaligned hierarchies or very eccentric orbits.
3. Non-universality of star formation in general, and multiplicity in particular, is not yet fully appreciated. The old concept of a universal IMF was traditionally extended in the past to assume a universal multiplicity. High multiplicity in the SFRs like Taurus, compared to the field, has been attributed to the dynamical processing in clusters. Now we know that the metallicity of the gas affects its fragmentation scale and leads to a strong inverse dependence of the close-binary fraction on metallicity (Moe et al. 2019). Conversely, the fraction of wide binaries positively correlates with metallicity (Hwang et al. 2021). At intermediate separations of 300-3000 au, the fraction of binaries in loose associations is twice lager than in the open clusters (Deacon & Kraus 2020). This difference highlights the dependence of multiplicity on the density of the formation environment (it would be difficult to explain the difference by dynamical processing, which is effective only at larger separations). The notions of binary-rich and binary-poor formation environments help us to understand certain aspects of hierarchies mentioned above.

Observational Perspective
Comparing predictions of formation theories with the observed statistics of stellar systems has always been the main way of testing the theory and gaining new insights. However, the population of stellar systems in the field results from a mixture of various formation channels. Moreover, statistical properties of binary and multiple systems depend substantially on their formation environments, as emerges from the recent studies. This consideration somewhat undermines the usefulness of statistics derived from the nearby field population.
On the other hand, a reasonable completeness over the full range of periods can be reached only for the nearby field. Specifically, the multiplicity statistics of the nearby low-mass stars and brown dwarfs can be revealing because the DI formation channel is unlikely and the accretion-induced evolution of the orbits may be less important. Thus, the architecture of low-mass hierarchies is more directly related to the CF-based formation scenarios. Two nearby quintuples with stellar masses of ∼0.4 M ⊙ are mentioned in section 4.6. The triple LHS 1070 and the young quadruple with substellar masses studied by Bowler & Hillenbrand (2015) are extreme examples of such systems. Discovery and study of more low-mass hierarchies are needed.
The origin of stellar hierarchies is also elucidated by the detailed study of selected 'benchmark' systems, linking their properties to the formation theories. This approach, most popular in the study of nascent stars (section 4.7), can be usefully extended to mature stellar hierarchies because attempts to explain their architecture might be challenging and informative at the same time (e.g. Powell et al. 2021). Here I explored a hybrid strategy by dividing known hierarchies into groups according to some relevant parameters and associating these groups with different formation scenarios. This classification attempt can be rightly criticized for being subjective and ambiguous. For example, Castor is a large-N system, it contains a 2+2 inner quadruple, and this quadruple is attributed to the special 'Castor family' presumably formed by the CF+DI2 scenario; λ Tau is compact, but, at the same time, likely misaligned.
Accumulation of new data on stellar hierarchies is currently progressing at an increased pace. The classical, object-by-object studies of hierarchies are nowadays complemented or replaced by large-scale surveys from space, e.g. Gaia (Gaia Collaboration et al. 2018) and TESS (Ricker et al. 2014), that detect new hierarchical systems in large numbers. The ground-based photometric and spectroscopic surveys also make a substantial contribution (e.g. Zasche et al. 2019). Well-defined coverage and uniform detection limits of surveys favor their use for statistical studies. However, the sky surveys are usually not designed for the study of hierarchies and bring only partial information that helps to discover such systems but is often insufficient for their study. For example, the cadence of the Gaia RV data will allow determination of spectroscopic orbits in a restricted range of periods and only for single-lined systems, because splitting blended lines requires a larger spectral resolution and dedicated methods. Similarly, the Gaia mission duration will restrict the periods of future Gaia astrometric orbits to less than a few years.
Combination of surveys with dedicated follow-up observations of selected objects is nowadays the mainstream of observational progress on multiples. By screening large samples, surveys help to identify rare systems that either extend the explored parameters space (e.g. compact and ultra-planar systems found in Kepler and TESS) or offer prospects of particularly detailed study. Follow-up observations of these selected systems reveal their structure in detail (e.g. Powell et al. 2021). Highresolution studies of cold gas by ALMA add a new dimension to the classical field of multiple stars, especially promising for the PMS hierarchies (see section 4.7).
The uniform coverage and the large volume of modern surveys can partially compensate for the insufficient information on individual systems by using statistical approach, where the statistics of the underlying population are inferred by modeling the observed parameters (e.g. astrometric accelerations in the future Gaia data releases) and solving the inverse problem. For example, thousands of eclipsing binaries discovered in photometric surveys, coupled with modeling of the observational selection, provide solid information on the underlying population of close binaries. Statistical exploitation of large surveys in terms of stellar multiplicity is expected to augment in the near future.
Many thousands of new binaries discovered by Gaia and their accurate astrometry can contribute substantially to the study of relative orbit orientation by the indirect method of counting co-and counter-rotating systems; one could probe mutual inclination as a function of mass and age. Relative motions in wide pairs also contain statistical infromation on their eccentricities, allowing to test the unfolding mechanism of Reipurth & Mikkola (2012). Intriguingly, Shatsky (2001) found that wide pairs containing subsystems tend to have less eccentric orbits than simple binaries, in direct contradiction with the predictions of unfolding.

Theoretical Perspective
The way forward from the current tentative classification of hierarchies based on the biased MSC sample to a more robust and informative connection between architecture of hierarchies and their formation mechanisms is, as usual, offered by combination of theoretical advances with new observations. I am not an expert in the star formation theory and can only hope for more quantitative and predictive studies of different formation channels. The inherently chaotic nature of star formation and the diversity of initial conditions and relevant processes are natural obstacles in this regard. Although all simulations are informative, the greatest return is expected by simulating many random cases followed by statistical analysis of the results. In this way, predictions of the formation scenarios can be sharpened and quantified, allowing better comparison with observations. Establishing the relative role of different formation channels and quantitative predictions of their outcomes is one of the major outstanding issues.
Formation of close binaries and their relation to higherorder multiples is still an open question. Although the crude modeling by Tokovinin & Moe (2020) and the analysis by Moe & Kratter (2018) highlight accretion as the main agent of orbital migration, many details are still obscure: Is migration in triples always associated with orbit alignment or not? Can simulations reproduce the formation of compact 2+2 hierarchies?
Even if in the future the goal of building a predictive theory of multiple-star formation will prove elusive (as it is now), every step in this direction will bring new valuable insights on the formation of stars and planets.
I appreciate comments on the draft of this review provided by N. Shatsky and S. Rappaport. This work used the SIMBAD service operated by Centre des Données Stellaires (Strasbourg, France), bibliographic references from the Astrophysics Data System maintained by SAO/NASA, and the Washington Double Star Catalog maintained at USNO.