11institutetext: Niels Bohr Institute, University of Copenhagen, Øster
Voldgade 5-7, DK-1350, Copenhagen K, Denmark
11email<EMAIL_ADDRESS>
# The contribution of binary star formation on core-fragmentation scales on
protostellar multiplicity

Rajika L. Kuruwita 11 Troels Haugbølle 11

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. Observations of young multiple star systems find a bimodal
distribution in companion frequency and separation. The origin of these peaks
have often been attributed to binary formation via core and disc-
fragmentation. However, theory and simulations suggests that young stellar
systems that form via core-fragmentation undergo significant orbital
evolution.

Aims. Using simulations of star formation in giant molecular clouds we
investigate the influence of environment on multiple star formation pathways
and the contribution of core-fragmentation is on the formation of close
($<100\,\mathrm{AU}$) binaries.

Methods. Simulations are run with the adaptive mesh refinement code RAMSES
with sufficient resolution to resolve core-fragmentation beyond
$400\,\mathrm{AU}$ and dynamical evolution down to $16\,\mathrm{AU}$, but
without the possibility of resolving disc-fragmentation. The evolution of the
resulting stellar systems is followed over millions of years.

Results. We find that star formation in lower gas density environments is more
clustered, but despite this, the fractions of systems that form via dynamical
capture and core-fragmentation are broadly consistent at $\sim$40% and
$\sim$60% respectively. In all gas density environments, we find the typical
scale at which systems form via core-fragmentation is
$10^{3-3.5}\,\mathrm{AU}$. After formation, we find that systems that form via
core-fragmentation have slightly lower inspiral rates ($\sim
10^{-1.75}\,\mathrm{AU/yr}$ measured over first $10000\,\mathrm{yr}$) compared
to dynamical capture ($\sim 10^{-1.25}\,\mathrm{AU/yr}$). We then compared the
simulation with conditions most similar to the Perseus star forming region to
determine whether the bimodal distribution observed by Tobin et al. (2016b)
can be replicated. We find that it can be replicated, but it is sensitive to
the evolutionary state of the simulation.

Conclusions. Our results indicate that a significant number of binary star
systems with separations $<100\,\mathrm{AU}$ can be produced via non-disc-
fragmentation pathways due to efficient inspiral, suggesting disc-
fragmentation is not the dominant formation pathway for low-mass close
binaries in nature.

###### Key Words.:

Star Formation – Binary stars; Simulations – MHD

##  1 Introduction

Around half of all stars exist in binary or multiple star systems (Moe & Di
Stefano 2017), and many are born with a companion (Chen et al. 2013). Models
of star and planet formation must account for this multiplicity.

The main proposed pathways for binary star formation are:

  1. 1.

Disc-fragmentation ($a\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}100\,\mathrm{AU}$): can
occur in massive, gravitationally unstable (i.e. Toomre Q parameter $<1$)
protostellar discs, if the gravitationally unstable regions are able to
efficiently cool (Gammie 2001). This has been modelled in simulations
(Takaishi et al. 2021), and with high angular resolution imaging, it is
possible to see what could be the outcome of disc-fragmentation
observationally (Tobin et al. 2016a). Bate (2019) investigated the dependency
of metallicity on fragmentation and find in their models that $15-20\%$ of
systems formed via this pathway.

  2. 2.

Core-fragmentation ($100\,\mathrm{AU}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
3.0pt\hbox{$\sim$}}\hss}\hbox{\raise 2.0pt\hbox{$<$}}}}a\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 3.0pt\hbox{$\sim$}}\hss}\hbox{\raise
2.0pt\hbox{$<$}}}}10000\,\mathrm{AU}$): occurs in turbulent protostellar
cores. Simulations of core-fragmentation into binaries have often displayed
significant orbital migration from their initial separation (Ostriker 1999;
Wurster et al. 2018; Lee et al. 2019; Kuruwita et al. 2020; Saiki & Machida
2020), and separations can shrink due to the ejection of tertiary companions
(Armitage & Clarke 1997) down to $<100\,\mathrm{AU}$. Observations of wide
binaries with misaligned discs have also been attributed to formation via
core-fragmentation (Lee et al. 2017).

  3. 3.

Dynamical Capture: In dense stellar environments it is possible for stars that
are marginally unbound to become bound or vice versa (Parker & Meyer 2014).
There is observational evidence to support that primordial wide binaries
formed in clusters can disintegrate (Elliott & Bayo 2016). There is also
evidence suggesting a third of multiple stars did not form together and
instead became bound via dynamical interactions (Murillo et al. 2016).

Observations of the companion frequency distribution of protostars have found
a bimodal distribution (Tobin et al. 2016b, 2022). For the Perseus molecular
cloud, Tobin et al. (2016b) finds a peak at approximately $75\,\mathrm{AU}$
and $3000\,\mathrm{AU}$. The authors suggested that the peak at
$75\,\mathrm{AU}$ is caused by disc-fragmentation, while the peak at
$3000\,\mathrm{AU}$ is likely caused by fragmentation of protostellar cores.
They hypothesise that core-fragmentation is a dominant mechanism for system
with separations $>200\,\mathrm{AU}$. However, the authors also hypothesise
that the lack of binaries with separations $>1000\,\mathrm{AU}$ may be due to
rapid in-spiralling to lower separations.

Surveys of circumstellar disc sizes find the average radius to be $\leq$75 au
(Cox et al. 2017; Ansdell et al. 2018) for low-mass young stars and the very
large and massive discs that could fragment to form stars are uncommon. While
there is evidence that disc-fragmentation can form close binaries, it is
unlikely to be the sole source of binaries with separations of $\sim$10s of
au.

In this work we investigate the contributions of core-fragmentation and
dynamical capture to the origin of the observed bimodal separation
distribution. We run multiple simulations of star formation from giant
molecular clouds with different initial masses, to investigate how the
environment may influence formation pathways for multiple star systems, and
orbital evolution of these systems.

In Section 2 we cover the RAMSES code and the simulation set up. In Section 3
we look at the formation and evolution of the multiple star systems that form,
including a comparison with observations. In Section 4 we discuss the overall
results of this work and the implications. In Section 5 we discuss some of the
limitations of this paper and suggestions for future work.

$M_{gas}$ |  $T_{end}$ ($T_{5}$) |  SFE |  $M_{\star}$ ($M_{\star,5}$) |  $N_{\star}$ ($N_{\star,5}$)
---|---|---|---|---
M⊙ |  Myr |  % |  M⊙ |
1500 |  26.8 (26.8) |  5 |  75 (75) |  86 (86)
3000 |  25.0 (23.8) |  13 |  393 (150) |  413 (191)
3750 |  23.7 (23.5) |  7 |  252 (187) |  411 (313)
4500 |  23.5 (23.2) |  8 |  341 (225) |  506 (394)
6000 |  23.0 (22.8) |  8 |  461 (300) |  895 (624)
12000 |  22.5 (22.4) |  8 |  923 (600) |  2474 (1795)
Table 1: Columns show the initial gas mass ($M_{gas}$), simulation duration
($T_{end}$), final star formation efficiency (SFE), final total mass in sink
particles ($M_{\star}$) and final number of sink particles produced
($N_{\star}$). The values in parenthesis are the corresponding column value at
SFE=0.05.

##  2 Method

###  2.1 Simulation set-up

We use a locally developed version of the publicly available
magnetohydrodynamics (MHD) Adaptive mesh refinement (AMR) code Ramses
(Teyssier 2002), which is described in Haugbølle et al. (2018). The simulation
setup is a $4\,\mathrm{pc}^{3}$ box with an isothermal equation of state and
solenoidal turbulent driving resulting in a typical velocity dispersion of 2
km s-1, with periodic boundary conditions. Six simulations were run with
initial gas masses of $1500,3000,3750,4500,6000$ and
$12000\,\mathrm{M}_{\odot}$. The box contains a magnetic field of
$7.2\mu\mathrm{G}$, initially aligned with z-axis, corresponding to an
alfvenic Mach number of 5. The simulations ran until a star formation
efficiency (SFE) of at least 0.05 was reached. The SFE is defined as the
fraction of the initial gas mass that is accreted onto sink particles. The
refinement strategy is based on the overdensity. We refine the root grid when
the density reaches a threshold density such that the Jeans length is resolved
by 14.4 cells, and then refine the grid every time the density increases by a
factor of four, keeping the minimum number of cells per Jeans length constant
on all levels, except the highest level of refinement. Sink particles are
allowed to form on the highest level of refinement from gas that has reached a
density corresponding to resolving the Jeans length with at least 2 cells. The
maximum resolution is $50\,\mathrm{AU}$, and the density threshold for sink
formation is $1.7\,10^{9}$ cm-3. Details of each simulation are summarised in
Table 1.

####  2.1.1 Sink particle model

A sink particle model is implemented to capture star formation. A summary of
the sink particle formation parameters is below. For a detailed description of
the sink particle model, please refer to Haugbølle et al. (2018). In order for
a sink particle for form within our simulations, the collapsing gas must meet
the following criteria:

Figure 1: _Top_ : Total mass accreted in sink particles, _Middle_ : number of
sink particles, _Bottom_ : Star formation efficiency over time. The circles
indicated the value of a quantity when the simulation has SFE=0.05.

Figure 2: Projections of the density over the entire computational domain from
the formation of the first sink (_Left column_) to SFE=0.05 (_Right column_),
for each simulation. The sink particles are annotated with blue particles. The
colour bars are logarithmically centred on the mean gas density of the initial
gas distribution.

Figure 3: Illustrates the three formation pathways that a system with a new
sink ID is classified into: _Top_ : If a sink is formed and is gravitationally
bound to another star at birth, this is bound core-fragmentation, _Middle_ :
If a sink is formed and is gravitationally unbound at birth, but later becomes
bound to the sink it was most bound to at birth, this is unbound core-
fragmentation, and _Bottom_ : If a sink is formed and is gravitationally
unbound at birth, and later becomes bound to a sink that is not the sink it
was most bound to at birth, this is dynamical capture.

Figure 4: Separation vs SFE of all systems formed in each simulation. The
markers indicate the formation pathway of the system if it contains a new sink
id. _Blue boxes_ : Bound core-fragmentation, _purple triangles_ : unbound
core-fragmentation, _red circles_ : dynamical capture.

  1. 1.

The cell is on the highest level of refinement.

  2. 2.

The density within a cell must be above the threshold number density of
$\rho_{sink}=1.7\times 10^{9}$cm${}^{-}3$. A typical value for $\rho_{sink}$
in the runs presented here is $10^{5}$-$10^{6}$ times the average density.

  3. 3.

There must be a gravitational potential minimum in the cell where the sink
particle will form.

  4. 4.

The velocity field is converging at the cell.

  5. 5.

No other sink particle exist within an exclusion radius. For our simulations
this is set to 8 cells or $400\,\mathrm{AU}$.

This sink particle prescription is similar to others used in different MHD
codes (Bate et al. 1995; Krumholz et al. 2004; Federrath et al. 2010; Gong &
Ostriker 2013).

When a sink particle forms, it initially has no mass, but mass is accreted
from cells within the accretion radius ($r_{acc}$). In our simulations
$r_{acc}=$4 cells. The accretion rate, $\dot{m}$, in the time step is
dependant on the gas density, distance of the gas to the sink particle, and
relative speed; see Haugbølle et al. (2018) for details.

Protostellar stellar jets and outflows carry away a portion of the mass and
momentum that is accreted to the inner disc regions (Tomisaka 2000, 2002;
Federrath et al. 2014; Kuruwita et al. 2017). However, the launching region of
these high velocity outflows are not resolved in our simulations. Therefore to
account for mass loss via outflows, an accretion efficiency,
$\epsilon_{\textrm{acc}}$, is defined, such that the actual accreted mass is
given by $\epsilon_{acc}\dot{m}\rm{dt}$, where $\rm{dt}$ is the time step. In
our simulations $\epsilon_{\textrm{acc}}=0.5$.

Sink particle positions and velocities are updated using a leap-frog
integrator. The smoothing length for the gravitational potential of the sinks
is set to be $0.33\Delta x/\approx 17\,\mathrm{AU}$.

With the selected simulations parameters and sink particle model, our
simulations resolve core-fragmentation but not disc-fragmentation. This choice
is necessary to limit the already very high computational cost, makes it
reasonable to use an isothermal equation of state. Furthermore, it is a
perfect setup to study the formation and evolution of multiple star systems
produced exclusively via core-fragmentation and dynamical capture.

##  3 Results and discussion

The evolution of the total accreted mass, number of stars, and SFE of all
simulations are shown in Figure 1, and values are summarised in Table 1.
Figure 1 shows that simulations with higher initial gas mass have a higher
star formation rate (Padoan & Nordlund 2011; Padoan et al. 2012; Federrath &
Klessen 2012) and also a higher rate of stars formed (Haugbølle et al. 2018)
in accordance with theory. The high mass simulations also begin to produce
sink particles earlier than the lower mass simulations, with the first sinks
in the $12000$ M⊙ simulation forming around $22$ Myr, while the $1500$ M⊙
simulation begins forming its first sinks around $23$ Myr, as shown in Table
1.

Projections of the density over the entire computational domain for each
simulation is shown in Figure 2. Each row is a simulation and each column
shows a projection for SFE=0%, 1.7%, 3.3%, and 5% respectively. The colour
bars are centred logarithmically on the initial mean gas density of each
simulation (defined as $M_{gas}/4$pc3), and have the same range. This was
chosen to visualise degrees of gas concentration between simulations. The sink
particles are annotated by blue dots. From the projections we see that most
sink particles reside in regions of high column density, which, in particular
at early times, is close to their region of birth.

###  3.1 Identifying multiple star systems and their formation pathways

As the sink particles accrete and interact, multiple star systems form and
disintegrate. Multiple star systems were identified by:

  1. 1.

Finding the sink particle pairs that are below a separation threshold and
gravitationally bound (where the $E_{potential}+E_{kinetic}=E_{total}<0$). For
our work, the maximum separation is set to $10\,000\,\mathrm{AU}$.

  2. 2.

Sort pairs from lowest to highest separation.

  3. 3.

A new ‘particle’ is saved with the position being the centre of mass of the
bound sink particle pair and mass that is the total mass of the particle pair.

This process is carried out iteratively, until all bound systems with up to
six components are identified. Six components is chosen as the maximum number
of stars in a system, because higher multiples systems tend to be short-lived
with very weakly bound oter companions in a highly dynamical environment, such
as close to the hub-like clusters in the simulations. Sextuples are also the
largest systems identified in the observations by Tobin et al. (2016b). When
comparing with observations we change the approach for identifying multiple
systems, since three dimensional distance and energies are not available for
observations. This is discussed in Section 3.5.

When iterating, to avoid double counting sink particles that are already in a
bound systems, only systems are used to calculate pairs. These are either an
actual sink particle if it is single, or one of the systems that was created
to collapse a bound pair.

As highlighted in 2.1.1 with given the numerical resolution in our models, we
can only capture core-fragmentation. However, some sink particles can form
isolated and unbound, but later become bound to another star. We use the birth
conditions of a sink particle found in a binary/multiple to determine the
pathway in which the system formed. We chronologically check bound systems,
and if a new sink particle ID is found, we determine the system’s formation
pathway using the following criteria:

  1. 1.

Bound core-fragmentation: If a sink is form bound to another sink or system of
sinks, then it was formed via bound core-fragmentation.

  2. 2.

Unbound core-fragmentation: If a sink is formed unbound, but the sink or
systems of sinks it was most bound to at birth is in the current multiple
system, then it was formed via unbound core-fragmentation.

  3. 3.

Dynamical capture: If a sink is formed unbound, and the sink or system of
sinks it was most bound to at birth is _not_ in the current multiple system,
then it was formed via dynamical capture.

These formation pathways are illustrated in Figure 3. For systems with more
than two stars, the formation pathways is determined by the newest sink
particle ID ie. the most recently formed star. This classification only rely
on the energetics of the sink particles, and not that of the gas. Therefore,
some of the systems found to be formed by unbound core-fragmentation may in
reality be bound, when making a detailed account of the gas distribution. We
have chosen not to do this, since it would entail relating different volume
elements to different sink systems, which is not easily done in a unique
manner, and beyond the scope of this article.

In Figure 4, we plot the separations of all multiple systems for each
simulation. For binaries, this is the true separation, but for higher order
multiples, this may be the separation between stars and the centre of mass of
an inner system, or two centres of masses, depending on the system
configuration. If a new sink particle ID is found in a system, we also
annotate its formation pathways with a marker at the initial separation of the
system. In Figure 4, the blue boxes indicate systems that formed via bound
core-fragmentation, the purple triangles indicate systems that formed via
unbound core-fragmentation, and the red circles indicate systems that formed
via dynamical capture. We see that all formation pathways occur throughout the
simulations. We also see in each panel that many systems experience rapid
orbital shrinkage, where the final separations of the systems is more than an
order of magnitude smaller than the initial separation.

Figure 5: Fractions of systems formed via the pathways (illustrated in Figure
3) for each simulation: bound core-fragmentation (blue), unbound core-
fragmentation (purple), and dynamical capture (red).

In Figure 5 we plot the fractions of systems that form via the different
pathways for each simulation. The fractions are broadly consistent over all
simulations. The fraction of systems formed via bound core-fragmentation
varies between $35-45\%$, unbound core-fragmentation varies between $10-25\%$,
and dynamical capture varies between $35-50\%$.

When looking at only the systems that are born unbound, in the lowest mass
simulations, more than a third of these systems follow the unbound core-
fragmentation pathway. This drops to $20\%$ in the highest mass simulations.
The inverse correlation between whether a star that is born unbound and later
becomes bound to the star it was most bound to, or another star, may be due to
the increasing stellar density with gas mass. In higher stellar density
environments the likelihood of dynamical capture would increase.

Overall, the fractions of systems formed via core-fragmentation (bound and
unbound) is slightly higher in the lower mass simulations. This may be because
the lower mass simulations require a larger over density in order for
protostellar cores to collapse, and this results in more systems forming via
the core-fragmentation pathways. This implies that star formation in lower
mass GMC is more clustered than in higher mass one.

###  3.2 Measuring clustering with two point correlation function

In this section we measure three-dimensional two point correlation functions
(TPCF) to gain an insight into how clustered the sink particles distributions
are in each simulation. A TPCF quantifies clustering by measuring all possible
separations between points (for $N$ stars, the number of separations is
$N(N-1)/2$)) and binning the separations into logarithmic separation bins
($r$). This is calculated for the observed data points ($P_{\mathrm{DD}}(r)$)
and for randomly generated points that are uniformly distributed
($P_{\mathrm{RR}}(r)$). These two distributions are normalised using the
number of data points and randomly generated stars to obtain:

$\mathrm{DD}(r)=\frac{P_{\mathrm{DD}}(r)}{N_{\mathrm{D}}(N_{\mathrm{D}}-1)},\mathrm{RR}(r)=\frac{P_{\mathrm{RR}}(r)}{N_{\mathrm{R}}(N_{\mathrm{R}}-1)},$
(1)

where $N_{\mathrm{D}}$ and $N_{\mathrm{R}}$ is the number of data points and
randomly generated points respectively. The basic TPCF (Peebles 1980) is
derived by dividing the normalised distribution of the data (DD$(r)$) by the
normalised distribution of the random points (RR$(r)$), i.e.:

$1+\omega(r)=\frac{\mathrm{DD}(r)}{\mathrm{RR}(r)}$ (2)

Figure 6: TPCFs of all simulations at SFE=0.05

If the data is clustered, the resulting TPCF will have values $\gg 1$ in
smaller separations and typically has a power law function as a function of
separation with a negative slope.

In order to calculate a TPCF an underlying distribution must be derived, for
comparison. Given the periodic boundary conditions of our simulation, the
largest separation two sinks can have is $\sqrt{3(L/2)^{2}}\sim$3.5 pc which
is equivalent to $\sim$7$\times 10^{5}\,$au. With this largest separation and
the smoothing length of 16.6 au, 10 separation bins spaced logarithmically
between $10-10^{6}\,$au were selected to calculate the TPCF. In order to
derive the underlying distribution, $5\times 10^{5}$ positions were randomly
generated within the simulation domain, and separations were calculated. The
large number of randomly generated positions was necessary in order to
populate the smallest separation bin. 100 instances of $5\times 10^{5}$
randomly generated positions was performed, and the final underlying
distribution ($P_{\mathrm{RR}}(r)$) was taken to be the mean of these 100
instances, and the error is the standard deviation in a separation bin over
the 100 instances.

For the simulations, the distribution ($P_{\mathrm{DD}}(r)$) is derived by
calculating all separations between stars and binning the separations. The
error of this distribution is taken to be the Poisson noise, because of the
low number of stars per separation in the simulation, i.e.
$\sigma_{P_{\mathrm{DD}}(r)}=\sqrt{N}$.

The TPCFs for all simulations at SFE $=0.05$ is shown in Figure 6. The error
bars for the TCPF is derived by summing the relative errors of DD$(r)$ and
RR$(r)$. After calculating the TPCF a power-law is fit to the function, with a
steeper slope indicating stronger clustering. At SFE $=0.05$, the lower mass
simulations show stronger clustering, however, this is only at one snapshot in
time.

To understand how the clustering evolves over the course of the simulations,
TPCFs were calculated for every time step over each simulations, and power-law
fits were derived. In Figure 7, we plot the median gradient of the power-law
for each simulation smoothed over a window of SFE $0.05\%$, and the shaded
region shows the standard deviation. Initially the gradient is not well
constrained due to the low number of sink particles. After SFE $\approx 0.02$,
the $M_{gas}=1500\,$M⊙ simulation consistently has stronger clustering over
the course of the simulations, while the higher mass simulations typically
have shallower gradients. This is consistent with the hypothesis that the
lower mass simulations produces more multiple star systems via core-
fragmentation pathways due to stronger clustering.

Observationally, the two-point correlation functions in the nearby Taurus,
$\rho$ Oph, and Orion regions are found to be a power law with exponents in
the range $-2<\beta<-1$ (Gomezet al. 1993; Larson 1995; Gomez & Lada 1998),
which compares well to our results, given that exponent should be increased by
one when projecting from three to two dimensions.

Figure 7: Gradient of two point correlation function (TPCF) over each
simulation. The gradient was smoothed over a window of SFE$\pm 5\times
10^{-4}$. The solid line shows the median of this smoothing and the shaded
regions show the $1\sigma$ variation of the smoothing.

###  3.3 System formation scales

Figure 8: Histograms of the initial separations of systems with new sink IDs
for each simulation, for the three different formation pathways: _Blue_ :
bound core-fragmentation, _purple_ : unbound core-fragmentation, and _red_ :
Dynamical capture. The separations were binned in log(Separation), and the
histogram for unbound core-fragmentation is stacked on top of the bound core-
fragmentation histogram. The black line is the fitted Gaussian to the core-
fragmentation distribution.

The star forming environment may affect the initial separation of systems
formed via core-fragmentation pathways, which might give insight into a
fragmentation scale. In Figure 8 we plot the histogram of the initial
separations of the systems formed via the three different formation pathways
for each simulation pathway. For the core-fragmentation pathways, the
histograms are stacked in the left column, while the initial separation of
dynamical capture is shown on the right column of Figure 8.

We fit Gaussian distributions to the sum of the histograms for systems formed
via bound and unbound core-fragmentation to characterise a typical bound core-
fragmentation scale. The means and widths of these distributions are plotted
against GMC mass in Figure 9. If we assume that the cores on average can be
described as critical mass Bonnor-Ebert spheres, the external pressure on the
cores will increase with increasing cloud mass, leading to the core
fragmentation scale depending inversely proportional to the square root of the
cloud mass, indicated by the dashed line in the figure. The typical
fragmentation scale over all simulations appears to be around
$10^{3\mbox{--}3.5}\,\mathrm{AU}$ which is consistent with observations that
find a peak of around $\sim 3000\,\mathrm{AU}$ in the separation distribution
of protostars (Pokhrel et al. 2018).

When looking at the initial separations of systems formed via dynamical
capture (shown on the right column of Figure 8), we see that for all GMC
masses the histogram grows with larger initial separations. For simulation
$12000\,\mathrm{M}_{\odot}$ the initial separation histogram plateaus beyond
$1000\,\mathrm{AU}$. This plateau is likely due to the higher stellar density
which creates a situation where it is more likely for stars to capture other
stars at smaller separations first before capturing stars at larger
separations.

Figure 9: Mean fragmentation scale as a function of the GMC mass. The curve
shows the peak of Gaussian distributions fitted to the (stacked) histogram of
initial separations of systems formed via core-fragmentation pathways plotted
in Figure 8. The errorbars is the width of the Gaussian distributions, and not
the error on the peaks. The dashed line show the theoretical prediction from
assuming that cores on average can be described as critical Bonnor-Ebert
spheres.

###  3.4 Orbital evolution of the different pathways

Simulations of binary star formation via protostellar core-fragmentation often
show that the binaries undergo significant orbital migration. There are two
classes of mechanisms that can efficiently provide migration. Dynamical
friction (Chandrasekhar 1943; Ostriker 1999), slow-down due to accretion of
material with lower angular momentum, and the initial high-eccentricity
configuration of newly-formed cores is a combination of reasons why naturally
binaries with separations below $\sim 1000\,\mathrm{AU}$ suffer rapid
inspiral. Prestellar cores form in general along filaments, which are
kinematically cold in the direction transverse to the filament. Therefore,
binaries formed through core fragmentation will initially have an angular
momentum, much lower than what is required for circular motion and launch on
highly eccentric orbits. This low angular momentum leads to rapid inspiral of
the young binary system. Systems formed via dynamical capture were likely to
be marginally unbound and became bound. When these systems first form, their
angular momentum is high, and thus, these systems are not expected to undergo
fast orbital evolution. On the other hand, dynamical capture provides for the
possibility of three-body interaction, which can lead to interaction and
hardening of the binary systems (Reipurth 2000).

To investigate inspiral, we measure the rate of change of the semi-major axis
in the $1000$ and $10\,000\,\mathrm{yr}$ after the first periastron of a
system. The semi-major axis is calculated at all time steps using:

$a=\frac{G\,M_{tot}}{2\epsilon_{orb}}\,,$ (3)

where $M_{tot}$ is the total mass of the system, and $\epsilon_{orb}$ is the
mass specific orbital energy of the system. $\epsilon_{orb}$ is calculated by
summing the specific kinetic and potential energies of the system.

The median inspiral rate for systems of different formation pathways is
plotted against simulation gas mass in Figure 10. The error bars indicate the
standard deviation of the measured inspiral rates. For both baselines, the
inspiral rates are not affected by the density if the star forming
environment. The inspiral rates are generally higher ($10^{-1}$au yr-1) in the
$1000\,\mathrm{yr}$ baseline compared to the $10\,000\,\mathrm{yr}$ baseline
($10^{-1.5}$au yr-1) for all pathways. This variation is not surprising
because we expect the inspiral rate to be the greatest soon after the
formation of the system, and then decrease at later stages.

For the $1000\,\mathrm{yr}$ baseline we do not see any significant variation
between formation pathways. For the $10\,000\,\mathrm{yr}$ baseline, we see
that the median inspiral rate of systems formed via dynamical capture is
consistently higher ($\sim 10^{-1.25}$AU yr-1) than bound core-fragmentation
($\sim 10^{-1.75}$AU yr-1).

The resulting inspiral rates measured from the $10\,000\,\mathrm{yr}$ baseline
is contrary to the previously stated hypothesis. This may suggest that systems
that form via dynamical capture undergo an interaction which hardens the
system quickly, and prolong the inspiral process, while bound core-
fragmentation systems reach a steady orbit sooner.

Overall, we see that many systems experience significant orbital evolution
within the early stages of their lifetime.

Figure 10: Shows the median inspiral rate of systems formed via the different
formation pathways (see Figure 3) measured over the $1000\,\mathrm{yr}$
(_Top_) and $10\,000\,\mathrm{yr}$ (_Bottom_) after the first periastron of
the system vs initial gas mass of each simulation. The error bars show the
standard deviation of the measured inspiral rates.)

###  3.5 Comparison with observations

We now aim to compare the multiplicity statistics produced from these
simulations with observations of protostellar multiplicity. In particular we
compare our results with that of Tobin et al. (2016b) of the Perseus star
forming region. Tobin et al. (2016b) observed protostars in this region and
measured how the companion frequency (CF) evolved with separation of the
components. The companion frequency (CF) as defined by Reipurth & Zinnecker
(1993) is the average number of companions a star has and is given by:

$CF=\frac{B+2T+3Q+...}{S+B+T+Q+...}.$ (4)

where $B,T,Q,...$ are the number of single, binary, triple, quadruple, and
higher order system respectively.

Tobin et al. (2016b) observed a bimodal distribution in the CF versus
separation, with peaks at $\sim$75$\,\mathrm{AU}$ and
$\sim$3000$\,\mathrm{AU}$ for Class 0 and I objects. The peak at
$75\,\mathrm{AU}$ was attributed to binaries formed via disc-fragmentation
while the peak at $3000\,\mathrm{AU}$ was attributed to systems formed via
core-fragmentation. However, we have shown in the previous sections that
multiple systems that formed on core-fragmentation scales often experience
significant orbital evolution, down to separations $\leq 75\,\mathrm{AU}$. We
aim to determine what the contribution of core-fragmentation is to the
observed bimodal distribution.

###  3.6 The Perseus star forming region and finding the best simulation for
comparison

The Perseus star forming region consists of multiple star forming clusters.
Arce et al. (2010) present a comprehensive analysis of the mass, volume,
turbulence, and star forming efficiency of various star forming regions in
Perseus. From Arce et al. (2010) the estimated current mass, including gas
mass, mass in outflows, and young stellar objects, in the L1448, NGC1333,
B1-Ridge, B1, IC348, and B5 star forming regions is approximately
$3220\,\mathrm{M}_{\odot}$. The estimated total volume of all star forming
regions is $53.5\,\mathrm{pc}^{3}$, which is slightly smaller than the
$4^{3}=64\,\mathrm{pc}^{3}$ volume of the simulations. Based on these mass and
volume estimates, the $M_{gas}=3000\,\mathrm{M}_{\odot}$ or
$3750\,\mathrm{M}_{\odot}$ simulations are best for the comparison.

To further refine which simulation is best for comparison, we calculate which
simulation produces a similar number of sink particles with luminosities that
would be observable by the observations of Tobin et al. (2016b).

The bimodal distribution found by Tobin et al. (2016b) is only prominent in
Class 0 and I objects, and it is not seen in Class II objects. Therefore, for
this comparison, we aim to target sink particles that would be classified as
Class 0/I. Accretion in protostars is generally higher in the early
protostellar stages than in the later stages, so an accretion limit of
$10^{-7}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ is applied to select sinks
that are likely to be in the Class0/I stages. Observationally, Class I objects
have observed accretion rates from $\sim$10-9 to
$10^{-6}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, however Class II are not
observed to have accretion rates above
$10^{-7}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ (Fiorellino et al. 2021). With
the selected accretion limit we may miss some Class I objects with lot
accretion rate, but avoid selecting Class II objects.

The number of Class 0/I objects observed by Tobin et al. (2016b) was $92$
stars, with luminosities between $0.1$ and $55\,\mathrm{L}_{\odot}$. In Tobin
et al. (2022), the upper limit of the Perseus observations is stated to be
approximated $120\,\mathrm{L}_{\odot}$. In order to find the best simulation
for comparison with these observations, we calculated the number of sink
particles in a simulation that would be observable based on both the highest
luminosity object ($55\,\mathrm{L}_{\odot}$), and the theoretical maximum
luminosity ($120\,\mathrm{L}_{\odot}$). We assume that the luminosity is
dominated by the accretion luminosity. To calculate the accretion luminosity
we use:

$L_{acc}=f_{acc}\frac{GM\dot{M}}{R_{\star}}$ (5)

Figure 11: Shows the total number of sink particles (black) and the number of
visible sink particles for maximum luminosity of $120\,\mathrm{L}_{\odot}$
(blue), and $55\,\mathrm{L}_{\odot}$ (orange) as calculated by Equation 5 over
SFE for each simulation. The horizontal dashed line indicates the number of
Class 0/I objects observed by Tobin et al. (2016b) (i.e. 92 stars).

where $\dot{M}$ is the mass accretion rate for the sink particle, $R_{\star}$
is the radius of the protostar, and $f_{acc}$ is the fraction of potential
energy from accretion that is converted into radiation. For these calculations
we take $R_{\star}=2\mathrm{R}_{\odot}$ for all protostars and $f_{acc}=0.5$
(1.e. 50% of potential energy is converted to radiation).

The evolution of the number of visible stars over time is shown in Figure 11.
We see that for all simulations, before SFE $\sim 0.005$ essentially all sink
particles are visible. However at later times, the total number of sink
particles and visible sink particles diverges. This is because older sink
particles may have stopped accreting or have very low accretion rates, but
newer sinks are forming which have higher accretion rates. In all simulations
a steady state in the number of visible sinks is established after SFE $\sim
0.01$. Based on the steady state visible star numbers and the number of
objects observed by Tobin et al. (2016b) (annotated by the horizontal dashed
line in Figure 11) the $M_{gas}=3750\,\mathrm{M}_{\odot}$ simulation produces
the closest number of visible stars to the observations. Therefore, for
further comparisons with observations, this simulation will be used.

####  3.6.1 Calculating characteristic SFE of Perseus

Observations capture the star formation evolution at a particular time, while
simulations evolve over time. To make an appropriate comparison with
observation, we need to find the best time of the simulations to compare with
observations. To do this, we estimate a characteristic star formation
efficiency of the observed star forming regions in Perseus.

Arce et al. (2010) provide right ascension and declination boundaries for the
L1448, NGC1333, B1-Ridge, B1, IC348, and B5 star forming regions and
observational estimates for the star formation efficiency in these regions.
Based on the boundaries of these regions, we find how many of the Class 0 and
I objects from Tobin et al. (2016b) are in each region. We find that 12 stars
are in L1448, 39 stars are in NGC1333, 4 stars are in B1-Ridge, 10 stars are
in B1, 13 are in IC348, 2 are in B5, and 12 are unclassified. We re-classified
the unclassified stars based on which star forming region they are closest to
on the sky, and this adds 11 stars to B1-Ridge and 1 star to B1.

We then calculate a weighted average of the star forming efficiency based on
the number of stars in each region and the observed SFE. Based on Arce et al.
(2010), the SFE of each region is 1.5% for L1448, 4.9% for NGC1333, 2.4% for
B1-Ridge, 2.1% for B1, 9% for IC348 and 0.4% for B5. From this, the
characteristic SFE we retrieve is $\sim 4.2\%$.

Figure 12: The median companion frequency vs. separation at star formation
efficiency of 4.2%, integrated over a SFE $\pm 0.1\%$ window, for the
$M_{gas}=3750\,\mathrm{M}_{\odot}$ simulation. The solid histograms in each
panel shows the resulting CF vs. Separation histogram for the different
settings summarised in Table 2. For settings with the luminosity limits, the
blue and orange show the results with an upper limit of $55$ and
$120\,\mathrm{L}_{\odot}$ respectively. The black dashed histogram shows the
observed distribution found by Tobin et al. (2022) for all class 0/I objects,
and the solid black line is the histogram derived from pairs with boundness
likelihood $>0.68$. The error-bars are calculated using the binomial
statistics described by equation 3 in Tobin et al. (2022).

###  3.7 Data processing to make CF vs Separation histogram.

In order to accurately carrying out a comparison with observations, we must
process our data the same way the observations were processed. While the true
sink particle separations and separations between centres of mass for sub-
systems was used to find bound systems previously, to be consistent with Tobin
et al. (2022), we also calculate the separations between mid-points of
components and sub-systems. Please refer to Figure 1 of Tobin et al. (2022)
for details on measuring separations in a multiple star system.

For some of the comparisons with observations, we save a projected midpoint-
separation, e.g., save the true separations projected onto the xy-plane. We
also vary whether accretion and luminosity limits are applied to be consistent
with Tobin et al. (2016b), and whether to count unbound pairs (pairs with
projected separation $<10,000\,\mathrm{AU}$). The various analysis settings
are summarised in Table 2. After all the multiple systems are found based on
the settings and limitations applied, we then process systems to determine how
the companion frequency evolves with projected separation.

Name | Projection | Bound |  $\dot{M}$ and $L_{acc}$
---|---|---|---
| into 2D | Check | Limits
3D-Full | False | True | False
3D-Limits | False | True | True
2D-Bound | True | True | True
2D-Unbound | True | False | True
Table 2: Settings used for CF vs. separation histogram. _Name_ : of set up,
_Projection into 2D_ : projects midpoint separation onto the xy-plane, _Bound
Check_ : only counts gravitationally bound pairs, _Limits_ : only counts
‘visible’ stars based on accretion and luminosity limits (see Section 3.6).

We create 12 separation bins logarithmically spaced from
$10^{1}-10^{4}\,\mathrm{AU}$, in the same fashion of Tobin et al. (2016a) and
Tobin et al. (2022). Within each bin, if a multiple system has components with
separations smaller than the lower bound of the bin, it is considered a single
star, and if there are components with projected separations larger than the
upper bound then that system is separated into smaller systems.

It is not appropriate to compare the observations with a single time step in
the simulations because accretion rates and luminosity can vary between time
steps. Therefore, the CF versus separation histogram was integrated over an
SFE $\pm 0.1\%$ window, and the median value in each bin is used to produce
the resulting histogram.

For comparison with observations, we retrieved the observations for Perseus
from Tobin et al. (2022). We removed objects that were not Class 0/I and
processed the data to create our own histogram. We use the entire data set to
produce the dashed histogram shown in Figure 12. To approximate the correction
for unbound pairs observed by Tobin et al. (2022), we recalculate the
histogram only for separation with boundness likelihoods $>0.68$ (Refer to
Table 4 of Tobin et al. (2022)). This produces the solid line histogram shown
in Figure 12. The error-bars are calculated using the binomial statistics
described by equation 3 in Tobin et al. (2022).

The results for SFE=4.2% is shown in Figure 12. For the 3D-Full settings, the
true CF distribution of all bound system appears to be relatively uniform at
CF=0.1. However, the CF drops significantly in the last separation bin,
indicating that there are fewer true bound systems at these large separations.
When observational limits are imposed in 3D-Limits, we observe an overall
reduction in CF in all bins. We also see a bimodal distribution appear with
peaks at $\sim 150\,\mathrm{AU}$ and $\sim 2000\,\mathrm{AU}$. When projected
onto 2D in 2D-Bound, the peak at $\sim 150\,\mathrm{AU}$ is still seen, but
larger peak is smeared out. Projecting separations into 2D will affect larger
separations that smaller separations, therefore it is expected that CF
histogram would be more affect at larger separations. When unbound pairs are
counted in the 2D-Unbound analysis we see that the CF increases significantly
above projected separations of $>100\,\mathrm{AU}$.

####  3.7.1 YSO density and contamination from chance alignments

Tobin et al. (2022) correct their observations for possible chance alignments.
The likelihood of a chance alignment is dependant on the observed young
stellar object (YSO) density. Based on our 2D-Bound and 2D-Unbound analysis,
Tobin et al. (2022) appear to have to sufficiently filtered out chance
alignments which may contaminate their results.

Our simulation produced significantly more unbound pairs than Tobin et al.
(2016a), and we hypothesise that the chosen simulation has a significantly
higher stellar density that Perseus. To investigate this we calculate the YSO
density distribution for our simulation at different stages of evolution.

The YSO density is calculated using the same method described by Tobin et al.
(2022). The YSO density around the sink particle is calculated using:

$\Sigma_{\mathrm{YSO}}=\frac{10}{\pi r_{11}^{2}}$ (6)

where $r_{11}$ is the separation to the eleventh nearest visible neighbour.

We measure the YSO density at SFE = 0.5%, 1%, 2%, 3%, 4%, and 5%. The
resulting cumulative distributions are plotted in Figure 13 against the
observed YSO density of Perseus for the Class 0/I objects. At SFE $>0.5\%$,
the YSO density measured in the simulations is larger than Perseus. In the
simulations, we observe that the maximum YSO density increases up to SFE=2%
and then decreases. This may reflect that initial burst of star formation in
clusters, and then the later dispersal of these clusters.

The higher YSO density in the simulations contributes to the significantly
higher CF measured with the 2D-Unbound set-up in Figure 12. The median YSO
density for both Perseus and our simulation at SFE=4% is approximately
100pc-2, however the simulation has an extended tail at higher densities. The
highest YSO density measured among the Perseus objects is $10^{2.8}$pc-2,
while the highest YSO density in the simulation at SFE=4% is $10^{3.8}$pc-2,
an order of magnitude greater. Therefore around some objects we may expected
up to 10 times more unbound pairs, contributing to the significant excess in
companion frequency seen in Figure 12.

Despite the significant increase in CF when counting unbound pairs, a peak in
CF is still visible in this set-up at $\sim$200$\,\mathrm{AU}$.

Figure 13: YSO density CDF for our simulations at different SFE (blue to pink
gradient). The measured YSO density for Tobin et al. (2016b, 2022) is shown by
the dashed black line.

####  3.7.2 Finding best fit using two-sample KS test

While we find a bimodal distribution at the calculated characteristic SFE, we
aim to determine if there is another time in the simulation where a stronger
bimodal distribution is found. We carry out a two-sample Kolmogorov–Smirnov
(KS) test (Massey 1951) because it is a non-parametric comparison of the shape
of two distributions. While we show the result of the bimodal distribution
found in Perseus, we are not concerned by the values of the resulting CF
versus separation histogram, but more specifically the shape.

In order to calculate the KS statistic, we:

  1. 1.

Order all the separations within systems observed in Perseus, and then
calculate a cumulative distribution function (CDF). This distribution is
normalised by dividing the CDF by the sum of all values (such that the bounds
are 0 and 1)

  2. 2.

For each time in the simulation, we used the 2D-Bound analysis to find the
visible stars. All separations in all bound systems are sorted to produce a
normalised CDF for the simulated systems.

  3. 3.

Calculate the KS statistic, which is the largest vertical separation between
the two normalised CDFs. We used the scipy package, which has the
scipy.stats.ks_2samp function to calculate the values.

We calculate the KS statistic between SFE=0.01 and 0.05, we also calculate a
corresponding critical value, below which, the two distribution are
significantly similar. The critical value is found using:

$KS_{crit}=\sqrt{-ln\left(\frac{\alpha}{2}\right)\times\left(\frac{1+\frac{m}{n}}{2m}\right)}$
(7)

where $m$ and $n$ are the number of values in each of the two samples (i.e.
the number of separations, which is 39 for Perseus), and $\alpha$ is derived
from the confidence level. We use $\alpha=0.99$, therefore if the KS statistic
is below the calculated critical value, the two samples are consistent to a
confidence level of 99%. We calculate KS statistics for the number for visible
stars using limits in 2D-Bound analysis, for both luminosity limits.

We smoothed the KS statistic and critical value over a window of SFE$=0.01\%$,
and the results are shown in Figure 14. The shaded regions show the 16th and
84th percentile of the integrated values. For both luminosity limits, we see
that bimodal distributions appear at various stages in the simulation. This is
not surprising because we expect many systems to evolve from large to small
separations, meaning that at some points in the simulation, the valley
observed between $\sim 300-1000\,\mathrm{AU}$ is filled in.

Based on the KS test, the time when we produce a distribution that is most
consistent with the Perseus observation is at approximately SFE$=4.1\%$ for
both $L_{max}=55$ and $120\,\mathrm{L}_{\odot}$.

Figure 14: The smoothed KS statistic measured over SFE using the 2D-Bound
settings (See Table 2) for upper luminosity limits of
$120\,\mathrm{L}_{\odot}$ (_top_) and $55\,\mathrm{L}_{\odot}$ (_bottom_). The
smoothing window is SFE$=0.01\%$. The median value is plotted with an opaque
line, and the shaded regions show the 16th and 84th percentile of the
integrated values. The critical KS statistic (calcuated with Equation 7 is
shown by the black line. Figure 15: Same as Figure 12, but for SFE=$4.1\%$

####  3.7.3 Best fit to observation

Based on the time of best fit found by the two sample KS test, we reproduce
Figure 12 for SFE=$4.1\%$. This is shown in Figure 15. Although it is only
slightly earlier than the CF versus separation histogram shown in Figure 12,
the bimodal distribution is more pronounced when luminosity limits are applied
in 3D-Limits, and 2D-Bound. While the values of the CF may vary between the
simulations and observations, the positions of the peaks and valley are in
good agreement.

From this we show that the observed bimodal distribution can be produced by
multiple star formation via non-disc-fragmentation pathways alone. This is not
to say that disc-fragmentation does not occur, but that a significant number
of binaries with separations $<100\,\mathrm{AU}$ can form on core-
fragmentation scales and migrate. The CF values from the simulation are
generally higher than observations, and in order to reduce the overall CF
value, more single stars are needed.

If our simulations resolved disc-fragmentation we expected a higher CF at
lower separations and this may shift the inner peak to lower separations.

##  4 Discussion

###  4.1 Star formation environment on formation pathways

The classical formation pathways of binary star formation considered mostly
disc and core-fragmentation, and dynamical capture was thought to not occur
often (Tohline 2002). However, simulations of clustered star formation and
observations of young stars may suggest that flybys may occur frequently.
Flybys by unbound stars has been proposed as the trigger of the excitation of
spiral arms observed in some circumstellar discs Pérez et al. (2016); Cuello
et al. (2019), and potentially trigger accretion bursts Borchert et al.
(2022). Pfalzner & Govind (2021) used N-body simulations of the stellar
dynamics in young clusters and determined that stellar flybys are probably
more common that initially expected in low mass clusters, and 10%-15% of discs
should show evidence of this. Pfalzner & Govind (2021) suggest that while the
overall stellar density is lower in low mass star forming regions, the
distribution of stars is more clustered, which can aid stellar interactions in
the central regions. In our work, we also confirm that star formation is more
clustered in lower density environments.

With the evidence that stellar flybys are not uncommon, it is not unreasonable
to suggest that some flyby events may turn into dynamical capture events. Our
work finds that for stars that are born unbound but later become bound, in
higher density environments, they are more like to follow the dynamical
capture pathways. Murillo et al. (2016) using SED fitting, inclination
effects, and outflows, measure how coeval young multiple star systems in
Perseus were. They find that approximately a third of the objects were non-
coeval, suggesting these systems did not necessarily form together. This
fraction is consistent with our $M_{gas}=3750\,\mathrm{M}_{\odot}$ simulation
which also found around a third of systems followed the dynamical capture
formation pathway (see Figure 5). Our simulations also suggested that the
fraction of systems that form via dynamical capture increases in higher
density environments.

While we find that the star forming environment may influence the formation
pathways of multiples, it does not seem to significantly affect the subsequent
orbital evolution. It is difficult to observe evidence of orbital evolution in
wide binaries ($a>1000\,\mathrm{AU}$) due to the long orbital periods for
individual systems. However, observations of OB associations find that the
velocity dispersion increases with age (Ramírez-Tannus et al. 2021). This
relationship between velocity dispersion and age is suggested to be caused by
the hardening of binaries, supporting orbital migration forming close
binaries. However, this velocity dispersion evolution is observed on the order
of $\,\mathrm{Myr}$, and but simulations of binary star formation from core-
fragmentation frequently produce binaries that experience significant orbital
evolution very early (Ostriker 1999; Wurster et al. 2018; Saiki & Machida
2020; Lee et al. 2019; Kuruwita et al. 2020).

The multiple star systems formed in the simulations of this paper also show
significant early orbital evolution (see Figure 10). Curiously, we find that
systems that formed via dynamical capture experienced faster inspiral rates,
when measured over the first $10\,000\,\mathrm{yr}$ of evolution since first
periastron. This was a surprise because we hypothesised that these systems
would inspiral slowly. However, many of these systems may have experienced an
interaction which quickly hardens previously unbound pairs.

###  4.2 Reproducing observations

In Section 3.5 we investigated whether our simulations reproduced the bimodal
CF versus separation distributions observed in protostars in Perseus (Tobin et
al. 2016b). This bimodal distribution has also been observed in Orion (Tobin
et al. 2022), and Encalada et al. (2021) appear to resolve the inner peak in
Ophiuchus. The valley between the two peaks has also been observed in Class I
objects in various star forming regions across the sky by Connelley et al.
(2008). Because this bimodal distribution has been observed in multiple
regions it is expected to be a long lived feature.

The origin of the two peaks have been attributed to disc-fragmentation
producing close binaries, and core-fragmentation producing wide binaries.
However, as shown throughout this paper, binaries that are formed on core-
fragmentation scales often experience significant orbital evolution early in
their formation. We calculated CF versus separation distribution of multiple
stars systems, looking at the true distribution, as well as applying
observational limits. We successfully produce bimodal distributions at the
calculated characteristic SFE of the observed objects in Perseus (4.2%), and
later find a better match to the Tobin et al. (2016b) observations at SFE =
4.1%.

As stated previously, this observed bimodal distribution may be long lived and
it is not clear if the evolutionary state of a star forming cloud would affect
when this distribution is observed. Multiple star systems form throughout all
of the simulations as shown in Figure 4, therefore we expect to constantly
observe systems forming at large separations ($>1000\,\mathrm{AU}$) and in-
spiralling to smaller separations. The derived best-fit is close to the
derived characteristic SFE, and it is not clear if this is coincidental.

While the best fit between simulations and observations occurs near when the
SFE of the simulation matched observations, when we look at the calculated KS
statistic against the critical value in Figure 14, we see that significant
fits appear throughout the evolution of the simulations. Looking at the
evolution of the KS statistic, bimodal distributions with similar shapes to
that observed in Perseus appear and disappear at different stages. There is an
extended period where significant fits are found between SFE $\sim$3.9-4.3%,
which spans $\sim$52$\,\mathrm{kyr}$ in simulation time. This bimodal
distribution will be visible for approximately 10% of the simulation run time.
This makes it clear that the feature in our simulation is not spurious, but is
a robust feature.

Overall, we can reproduce the observed bimodal distribution in protostellar CF
using multiple star systems formed on core-fragmentation scales only due to
orbital migration. This is not to say that disc-fragmentation does not occur,
but that a significant fraction of closer binaries (¡100AU) can form via non-
disc-fragmentation pathways.

##  5 Limitations and caveats

###  5.1 Numerical resolution

On the highest level of refinement, the resolution of our simulations is
$50\,\mathrm{AU}$, which means typical circumstellar disc sizes
($\sim$75$\,\mathrm{AU}$; Cox et al. (2017); Ansdell et al. (2018)) are only
resolved over a couple of cells. However, this work only investigates star
formation from molecular cloud fragmentation, therefore, resolving discs is
not crucial to this work.

Disc are an important part of the mechanism that determines how mass much is
accreted onto stars, and what is ejected via outflows, and this is discussed
in the next section.

The sink particle motion is calculated using an leap-frog integrator, and the
smoothing length is a third of a cell length on the highest resolution, i.e.
$16.6\,\mathrm{AU}$. The sink particle motion is calculated using the
gravitational potential from both other sink particles and the gas potential.
The integrator can accurately calculate the N-bobdy particles until the
separation is near the softening length. However, the gravitational potential
contribution from the gas can be less resolved as the two sink particle
approach separations near the softening length.

Overall, the numerical resolution is limited near and sink particles, and may
not accurately simulate orbital evolution when the separations are comparable
to the softening length. The work in this paper has primarily focused on
resolving the bimodal distribution observed in protostars, with the inner peak
at projected separations of 75$\,\mathrm{AU}$. With our simulation set up, we
are able to sufficiently resolve this inner peak.

###  5.2 Non-ideal MHD effects and resolving outflows

Our simulations compute the ideal-magnetohydrodynamic equations. For molecular
clouds, the typical fractional ionisation (abundance of electron) is
$\sim$10${}^{-6}-10^{-8}$. With this ionisation fraction, ambipolar diffusion
is present to help dissipate magnetic field flux to allow protostellar cores
to collapse. While we do not explicitly include ambipolar diffusion in the
computation of the MHD, Hennebelle & Inutsuka (2019) argue that the numerical
diffusion naturally present in hydrodynamic simulation is sufficient to
reproduce the effect of ambipolar diffusion. Other non-ideal effects are not
dominant in the molecular cloud regime, however, (Wurster et al. 2019)
investigated full non-ideal MHD with clustered star formation with smoothed
particle hydrodynamic simulations and find there is less magnetic breaking
leading to discs forming more easily. With our resolution, we do not resolve
discs, therefore, we do not expect the inclusion of other non-ideal effects to
affect our results.

Magnetic fields are responsible for launching jets and outflows which return
mass and momentum to the surrounding protostellar environment. With the
resolution of $50\,\mathrm{AU}$, jets are not self-consistently produced while
weak outflows can be produced with larger discs. Not resolving outflows will
not affect the motion of the sink particles significantly, but outflows
regulate what fraction of mass if accreted onto star and what is lost via
outflows. As described in 2.1.1, we assume a mass accretion fraction of
$50\%$. The fraction of mass that is accreted is not well constrained and
various disc wind models give accretion fractions from 40-90% (Seifried et al.
2012; Fendt & Sheikhnezami 2013). Observations of T-Tauri stars find mass
accretion fractions of 50-99% (Nisini et al. 2018), and observations of
Herbig-Haro objects find fractions of $\sim$90% (Ellerbroek et al. 2013).

The accretion fraction used in our simulations may be on the lower end of
typical accretion fractions but it is still consistent with models and
observation.

###  5.3 Radiation feedback

An isothermal equation of state is used in our simulations, therefore, the
internal energy remains constant throughout the simulations. This equation of
state is satisfactory for our simulations because with our resolution we do
not enter the regime where adiabatic heating of protostellar cores occurs
($\sim$3.8$\times$10${}^{-13}\,\mathrm{g}\,\mathrm{cm}^{-3}$; Masunaga &
Inutsuka (2000)). This is because the density threshold for sink particle
formation does not exceed $1.3\times 10^{-14}\,\mathrm{g}\,\mathrm{cm}^{-3}$.

While an isothermal equation of state is justified in this work for the
resolution at which we simulation hydrodynamics, the sink particles are used
as proxies for stars, which would be able to produce radiation feedback.
Radiation pressure from star formation can inject energy back into the star
forming environment, however, radiative MHD simulations carried out by Rosen &
Krumholz (2020) find that magnetic fields are dominant over radiation pressure
feedback, even in massive stars. Observations also suggest that gas
temperature may not have a strong impart on fragmentation and rather, mass and
density are key factors in fragmentation Murillo et al. (2018).

Radiation feedback from the protostars is also found to suppress disc-
fragmentation (Offner 2011). However, in this work we are only concerned with
multiple star formation on core-fragmentation scales.

##  6 Summary and Conclusion

We find that systems that form on core fragmentation scales undergo
significant orbital evolution. The density of the star forming environment may
influence the formation pathways, with star formation in low density
environments being more clustered leading to a higher fraction of systems
forming via core-fragmentation pathways rather than dynamical capture.

We are also able to reproduce the bimodal distribution observed in Perseus
(Tobin et al. 2016b). This is possible due to orbital migration from large
separations to $<100\,\mathrm{AU}$, suggesting that disc-fragmentation is not
a dominant formation pathways for binaries of close separations.

## Acknowledgements

The research leading to these results has received funding from the
Independent Research Fund Denmark through grant No. DFF 8021-00350B (TH,RLK).
This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No. 847523 ‘INTERACTIONS’. The astrophysics HPC facility at the
University of Copenhagen, supported by research grants from the Carlsberg,
Novo, and Villum foundations, was used for carrying out the simulations, the
analysis, and long-term storage of the results. yt (Turk et al. 2011) was used
to help visualise and analyse these simulations.

## References

  * Ansdell et al. (2018) Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21
  * Arce et al. (2010) Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170
  * Armitage & Clarke (1997) Armitage, P. J. & Clarke, C. J. 1997, MNRAS, 285, 540
  * Bate (2019) Bate, M. R. 2019, MNRAS, 484, 2341
  * Bate et al. (1995) Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362
  * Borchert et al. (2022) Borchert, E. M. A., Price, D. J., Pinte, C., & Cuello, N. 2022, Monthly Notices of the Royal Astronomical Society: Letters, 510, L37
  * Chandrasekhar (1943) Chandrasekhar, S. 1943, ApJ, 97, 255, aDS Bibcode: 1943ApJ….97..255C
  * Chen et al. (2013) Chen, X., Arce, H. G., Zhang, Q., et al. 2013, ApJ, 768, 110
  * Connelley et al. (2008) Connelley, M. S., Reipurth, B., & Tokunaga, A. T. 2008, AJ, 135, 2526
  * Cox et al. (2017) Cox, E. G., Harris, R. J., Looney, L. W., et al. 2017, ApJ, 851, 83
  * Cuello et al. (2019) Cuello, N., Montesinos, M., Stammler, S. M., Louvet, F., & Cuadra, J. 2019, A&A, 622, A43
  * Ellerbroek et al. (2013) Ellerbroek, L. E., Podio, L., Kaper, L., et al. 2013, A&A, 551, A5
  * Elliott & Bayo (2016) Elliott, P. & Bayo, A. 2016, MNRAS, 459, 4499
  * Encalada et al. (2021) Encalada, F. J., Looney, L. W., Tobin, J. J., et al. 2021, ApJ, 913, 149
  * Federrath et al. (2010) Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
  * Federrath & Klessen (2012) Federrath, C. & Klessen, R. S. 2012, ApJ, 761, 156
  * Federrath et al. (2014) Federrath, C., Schrön, M., Banerjee, R., & Klessen, R. S. 2014, ApJ, 790, 128
  * Fendt & Sheikhnezami (2013) Fendt, C. & Sheikhnezami, S. 2013, ApJ, 774, 12
  * Fiorellino et al. (2021) Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43, publisher: EDP Sciences
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174
  * Gong & Ostriker (2013) Gong, H. & Ostriker, E. C. 2013, ApJSS, 204, 8
  * Haugbølle et al. (2018) Haugbølle, T., Padoan, P., & Nordlund, Å. 2018, ApJ, 854, 35
  * Hennebelle & Inutsuka (2019) Hennebelle, P. & Inutsuka, S.-i. 2019, Frontiers in Astronomy and Space Sciences, 6
  * Krumholz et al. (2004) Krumholz, M. R., McKee, C. F., & Klein, R. I. 2004, ApJ, 611, 399
  * Kuruwita et al. (2020) Kuruwita, R. L., Federrath, C., & Haugbølle, T. 2020, A&A, 641, A59
  * Kuruwita et al. (2017) Kuruwita, R. L., Federrath, C., & Ireland, M. 2017, MNRAS, 470, 1626
  * Lee et al. (2019) Lee, A. T., Offner, S. S. R., Kratter, K. M., Smullen, R. A., & Li, P. S. 2019, ApJ, 887, 232
  * Lee et al. (2017) Lee, J.-E., Lee, S., Dunham, M. M., et al. 2017, Nature Astronomy, 1, 0172
  * Massey (1951) Massey, F. J. 1951, Journal of the American Statistical Association, 46, 68
  * Masunaga & Inutsuka (2000) Masunaga, H. & Inutsuka, S.-i. 2000, ApJ, 531, 350
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJSS, 230, 15
  * Murillo et al. (2016) Murillo, N. M., van Dishoeck, E. F., Tobin, J. J., & Fedele, D. 2016, A&A, 592, A56
  * Murillo et al. (2018) Murillo, N. M., van Dishoeck, E. F., Tobin, J. J., Mottram, J. C., & Karska, A. 2018, A&A, 620, A30
  * Nisini et al. (2018) Nisini, B., Antoniucci, S., Alcalá, J. M., et al. 2018, A&A, 609, A87
  * Offner (2011) Offner, S. S. R. 2011, Proceedings of the International Astronomical Union, 6, 231
  * Ostriker (1999) Ostriker, E. C. 1999, ApJ, 513, 252, aDS Bibcode: 1999ApJ…513..252O
  * Padoan et al. (2012) Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJL, 759, L27, publisher: American Astronomical Society
  * Padoan & Nordlund (2011) Padoan, P. & Nordlund, Å. 2011, ApJ, 730, 40, publisher: American Astronomical Society
  * Parker & Meyer (2014) Parker, R. J. & Meyer, M. R. 2014, MNRAS, 442, 3722
  * Peebles (1980) Peebles, P. J. E. 1980, The large-scale structure of the universe
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519
  * Pfalzner & Govind (2021) Pfalzner, S. & Govind, A. 2021, ApJ, 921, 90
  * Pokhrel et al. (2018) Pokhrel, R., Myers, P. C., Dunham, M. M., et al. 2018, The Astrophysical Journal, 853, 5
  * Ramírez-Tannus et al. (2021) Ramírez-Tannus, M. C., Backs, F., de Koter, A., et al. 2021, A&A, 645, L10
  * Reipurth (2000) Reipurth, B. 2000, The Astronomical Journal, 120, 3177
  * Reipurth & Zinnecker (1993) Reipurth, B. & Zinnecker, H. 1993, A&A, 278, 81
  * Rosen & Krumholz (2020) Rosen, A. L. & Krumholz, M. R. 2020, AJ, 160, 78
  * Saiki & Machida (2020) Saiki, Y. & Machida, M. N. 2020, ApJL, 897, L22
  * Seifried et al. (2012) Seifried, D., Pudritz, R. E., Banerjee, R., Duffin, D., & Klessen, R. S. 2012, MNRAS, 422, 347
  * Takaishi et al. (2021) Takaishi, D., Tsukamoto, Y., & Suto, Y. 2021, PASJ, 73, L25
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Tobin et al. (2016a) Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016a, Nature, 538, 483
  * Tobin et al. (2016b) Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016b, ApJ, 818, 73
  * Tobin et al. (2022) Tobin, J. J., Offner, S. S. R., Kratter, K. M., et al. 2022, ApJ, 925, 39
  * Tohline (2002) Tohline, J. E. 2002, ARA&A, 40, 349
  * Tomisaka (2000) Tomisaka, K. 2000, ApJ, 528, L41
  * Tomisaka (2002) Tomisaka, K. 2002, ApJ, 575, 306
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
  * Wurster et al. (2018) Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 475, 1859
  * Wurster et al. (2019) Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719


$\mathscr{G}_{\widetilde{A}}=({\rm V},\mathscr{E}_{\widetilde{A}})$. There are
two possibilities (i) $\\{1,2\\}\notin\mathscr{E}_{\widetilde{A}}$; (ii)
$\\{1,2\\}\in\mathscr{E}_{\widetilde{A}}$.

  * •

Proof under Case (i): if $\\{1,2\\}\notin\mathscr{E}_{\widetilde{A}}$, then

$\psi\big{(}[\widetilde{A},F_{12}]\big{)}:=\mathscr{G}=({\rm V},\mathscr{E})$

has neither self-loops nor multiple edges, and $|\mathscr{E}|\geq 1$. With the
notation $\widehat{A}=[\widetilde{A},F_{12}]$, we have
$\mathscr{G}=\psi(\widehat{A})$. Applying Lemma 14, there exists a finite
sequence of node pairs $(i_{1},j_{1}),\dots,(i_{w},j_{w})$ for some integer
$w\geq 1$, such that the cardinality of the edge set of graph
$\mathcal{H}_{i_{w}j_{w}}\big{(}\cdots\mathcal{H}_{i_{2}j_{2}}\big{(}\mathcal{H}_{i_{1}j_{1}}(\mathscr{G})\big{)}\big{)}$
is equal to one. Let
$\mathcal{H}_{i_{w}j_{w}}\big{(}\cdots\mathcal{H}_{i_{2}j_{2}}\big{(}\mathcal{H}_{i_{1}j_{1}}(\mathscr{G})\big{)}\big{)}=\big{(}{\rm
V},\big{\\{}\\{i^{*},j^{*}\\}\big{\\}}\big{)}$. We conclude that

$\displaystyle\psi\big{(}[F_{i_{w}j_{w}},\dots,[F_{i_{2}j_{2}},[\widehat{A},F_{i_{1}j_{1}}]]]\big{)}=\big{(}{\rm
V},\big{\\{}\\{i^{*},j^{*}\\}\big{\\}}\big{)}.$

Hence

$[F_{i_{w}j_{w}},\dots,[F_{i_{2}j_{2}},[\widehat{A},F_{i_{1}j_{1}}]]]=a^{*}F_{i^{*}j^{*}},$

where $a^{*}$ is the coefficient generated during the operation of the Lie
brackets. Since $F_{i_{s}j_{s}}\in{\rm S}^{z}$ for $s=1,\dots,w$, we have
$F_{i^{*}j^{*}}\in\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}$. In addition, with the definition of
$\widetilde{A}$, we have $F_{i^{*}j^{*}}\notin{\rm S}^{z}$. Therefore,
$\mathscr{G}_{\\{F_{i^{*}j^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}}$ has multiple edges. Proposition 1 now yields
$\big{\\{}\\{F_{i^{*}j^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$. This implies that
$\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$, and thus
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$.

  * •

Proof of Case (ii): in $\mathscr{G}_{\widetilde{A}}$, the nodes adjacent to
node $1$ except node $2$ are denoted as $i_{1},\dots,i_{m}$, and, the nodes
adjacent to node $2$ except node $1$ are denoted as $j_{1},\dots,j_{n}$, where
$m,n\geq 0$. If $m=0$, let $i\in{\rm V}$ and $i\neq 1\neq 2$. Then it is
evident that graph $\psi\big{(}[\widetilde{A},F_{1i}]\big{)}$ has neither
self-loops nor multiple edges, and the cardinality of its edge set is greater
than or equal to two. This is due to the fact that
$\\{1,i\\}\notin\mathscr{E}_{\widetilde{A}}$. If $m\geq 1$, then the edge set
of $\psi\big{(}[\widetilde{A},F_{12}]\big{)}$ contains edge $\\{2,i_{1}\\}$.
Let $i_{*}\in{\rm V}$ and $i_{*}\neq i_{1}\neq 1\neq 2$. We see at once that
edge $\\{i_{1},i_{*}\\}$ is not in the edge set of
$\psi\big{(}[\widetilde{A},F_{12}]\big{)}$. Therefore,
$\psi\big{(}[[\widetilde{A},F_{12}],F_{i_{1}i_{*}}]\big{)}$ has neither self-
loops nor multiple edges, and the cardinality of its edge set is greater than
or equal to one.

In summarizing, there must exist a matrix
$\widehat{A}\in\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}$ such that $\psi(\widehat{A})$ has neither self-loops
nor multiple edges, and the cardinality of its edge set is greater than or
equal to one. Analysis similar to that in the proof of Case (i) shows that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$.

The proof above gives more, namely
$\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$ if
$\widetilde{A}\neq 0$.

We next prove that if
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
S}$ has a cycle with an odd number of Red edges, then
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$. We continue to use the
definition of $\widetilde{A}$ in the above. Now $\widetilde{A}\neq 0$, which
is due to the fact that $\mathscr{G}_{\mathrm{S}^{z}}$ dose not have a cycle
with an odd number of Red edges, while
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
S}$ has. It follows easily that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm
LA}=\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$.

(Necessity) Assume that
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
S}$ has neither self-loops nor a cycle with an odd number of Red edges.
Proposition 1 shows that the Lie algebra generated by ${\rm S}$ is not equal
to $\operatorname{\mathfrak{su}}(n)$. In addition,
$\mathscr{G}_{\mathrm{S}^{z}}$ is complete, and it has neither a self-loop nor
a cycle with an odd number of Red edges. We must have
$\mathscr{E}_{A}\subseteq\mathscr{E}_{\mathrm{S}^{z}}$, because
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{{\rm
S}^{z}}$ has neither a self-loop nor a cycle with an odd number of Red edges.
This implies that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\\{{\rm S}\\}_{\rm
LA}\neq\operatorname{\mathfrak{su}}(n)$. Therefore, if the Lie algebra
generated by
$\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathrm{S}$
is equal to $\operatorname{\mathfrak{su}}(n)$, then
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
S}$ must have a self-loop or a cycle with an odd number of Red edges. We have
now completed the proof of Lemma 15. ∎

We are now in a position to present the proof of Theorem 4. Let
$\mathscr{G}_{\rm contr}^{{}^{\beta}}$ be connected. Since $\mathscr{G}_{\rm
drift}^{{}^{\alpha}}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
contr}^{{}^{\beta}}$ has a self-loop or a cycle with an odd number of Red
edges, there must exist a matrix $A\in\Sigma_{\rm
r}(\operatorname{\mathfrak{su}}(n)_{\alpha})$ such that
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
contr}^{{}^{\beta}}$ has a self-loop or a cycle with an odd number of Red
edges. Using Lemma 15, we see that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\operatorname{\mathfrak{so}}(n)_{\beta}\big{\\}}_{\rm
LA}=\operatorname{\mathfrak{su}}(n)$. From this we conclude that the system
(1) is structurally controllable on the Lie group ${\rm SU}(n)$ with respect
to the pair of zero patterns $\Sigma_{\rm zero}:=\big{(}\Sigma_{\rm
r}(\operatorname{\mathfrak{su}}(n)_{\alpha}),\Sigma_{\rm
f}(\operatorname{\mathfrak{su}}(n)_{\beta})\big{)}$.

### F. Proof of Theorem 5

The statement in Theorem 5 will be proved once we prove the lemma below.

######  Lemma 16.

Consider a subset
$\mathrm{S}\subseteq\mathpzc{B}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathpzc{C}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathpzc{D}$
with the associated multigraph
$\mathscr{G}_{\mathrm{S}}=(\mathrm{V},\mathscr{E}_{\mathrm{S}})$ and a matrix
$A\in\operatorname{\mathfrak{su}}(n)$ with the associated multigraph
$\mathscr{G}_{A}=(\mathrm{V},\mathscr{E}_{A})$. The Lie algebra generated by
$\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathrm{S}$
is equal to $\operatorname{\mathfrak{su}}(n)$ if the following conditions
hold:

$(i)$ Each connected component of $\mathscr{G}_{\mathrm{S}}$ contains at least
three nodes;

$(ii)$ $\mathscr{G}_{A}$ has no multiple edges and the union graph
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\mathrm{S}}$
is connected;

$(iii)$
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\mathrm{S}}$
has a self-loop or a cycle with an odd number of Red edges.

###### Proof.

If $\mathscr{G}_{\rm S}$ is connected, then by Lemma 15, the Lie algebra
generated by
$\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathrm{S}$
is equal to $\operatorname{\mathfrak{su}}(n)$. Now assume that
$\mathscr{G}_{\rm S}$ has $m$ connected components with $m\geq 2$. Let
$\mathscr{G}_{i}=(\mathrm{V}_{i},\mathscr{E}_{i})$ denote the $i$-th connected
component of $\mathscr{G}_{\rm S}$ for $i=1,\dots,m$. We continue to use the
definition of ${\rm S}^{k}$ in the proof of Proposition 1. Write
$\overline{{\rm S}}=\bigcup_{k=1}^{\infty}{\rm S}^{k}$. There must exist a
nonnegative integer $z$ such that ${\rm S}^{z}={\rm S}^{z+1}=\dots$, which
then implies $\overline{{\rm S}}={\rm S}^{z}$. By the definition of
$\overline{{\rm S}}$, we see that ${\rm S}^{z}$ spans $\\{{\rm S}\\}_{\rm
LA}$. Since $\mathscr{G}_{\rm S}$ has $m$ connected components,
$\mathscr{G}_{{\rm S}^{z}}$ has $m$ connected components too. In addition, for
any $u<v\in\mathrm{V}_{i}$, $B_{uv}\in{\rm S}^{z}$ or $C_{uv}\in{\rm S}^{z}$.

We continue to use the definition of $\widetilde{A}$ in the proof of Lemma 15.
As the union graph
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{{\rm
S}^{z}}$ is connected while $\mathscr{G}_{{\rm S}^{z}}$ is not, we always have
$\mathscr{E}_{\widetilde{A}}\neq\emptyset$, i.e., $\widetilde{A}\neq 0$. This
allows us to further conclude that

$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm
LA}=\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}.$ (11)

The remainder of the proof is divided into two steps.

Step 1. We first prove the statement when $\mathscr{G}_{\rm S}$ contains only
two connected components, i.e., $m=2$. Because
$\mathscr{G}_{A}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{\rm
S}$ is connected,
$\mathscr{G}_{\widetilde{A}}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{{\rm
S}^{z}}$ is also connected by the definition of $\mathscr{G}_{\widetilde{A}}$
and $\mathscr{G}_{{\rm S}^{z}}$. In addition, there must exist an edge
$\\{u,v\\}\in\mathscr{G}_{\widetilde{A}}$ such that $u\in{\rm V}_{1},v\in{\rm
V}_{2}$. Let $v_{11}$ be the node in $\mathrm{V}_{1}$ with $\deg(v_{11})=k>0$,
and $v_{i_{1}j_{1}},\dots,v_{i_{k}j_{k}}$ be its neighbors. For convenience,
we let $v_{i_{1}j_{1}},\dots,v_{i_{r}j_{r}}\in\mathrm{V}_{2}$, $1\leq r\leq
k$. Consider the node $v_{12}\in\mathrm{V}_{1}$. We must have
$F_{v_{11}v_{12}}\in{\rm S}^{z}$, where $F=B$ or $F=C$.

Let $v_{i_{*}j_{*}}\neq v_{i_{1}j_{1}}\in\mathrm{V}_{2}$. There holds
$F_{v_{i_{1}j_{1}}v_{i_{*}j_{*}}}\in{\rm S}^{z}$. Lemma 10 now shows that

$\psi\big{(}[[\widetilde{A},F_{v_{11}v_{12}}],F_{v_{i_{1}j_{1}}v_{i_{*}j_{*}}}]\big{)}:=\mathscr{G}=\big{(}{\rm
V}_{1}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
V}_{2},\mathscr{E}\big{)}$

is a bi-graph, with $|\mathscr{E}|\geq 1$. Write
$[[\widetilde{A},F_{v_{11}v_{12}}],F_{v_{i_{1}j_{1}}v_{i_{*}j_{*}}}]=\widehat{A}$.
It is immediate that $\mathscr{G}=\psi(\widehat{A})$. From Lemma 13, there is
an integer $z\geq 1$ and a finite sequence of node pairs either
$v_{s},u_{s}\in\mathrm{V}_{1}$ or $v_{s},u_{s}\in\mathrm{V}_{2}$ for
$s=1,2,\dots,z$, such that
$\mathcal{H}_{v_{z}u_{z}}\big{(}\cdots(\mathcal{H}_{v_{1}u_{1}}(\mathscr{G}))\big{)}=\big{(}\mathrm{V}_{1}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathrm{V}_{2},\mathscr{E}_{v_{z}u_{z}}\big{)}$
is also a bi-graph with $|\mathscr{E}_{v_{z}u_{z}}|=1$. Let
$\mathscr{E}_{v_{z}u_{z}}=\big{\\{}\\{v^{*},u^{*}\\}\big{\\}}$, where
$v^{*}\in\mathrm{V}_{1},u^{*}\in\mathrm{V}_{2}$. Note that
$\\{v_{s},u_{s}\\}\in\mathscr{E}_{{\rm S}^{z}}$, we have
$F_{v_{s}u_{s}}\in{\rm S}^{z}$ for $s=1,\dots,z$.

Based on Definition 10 and Lemma 10, it can be concluded

$\displaystyle\psi\big{(}[F_{v_{z}u_{z}},\dots,[F_{v_{2}u_{2}},[\widehat{A},F_{v_{1}u_{1}}]]]\big{)}=\mathcal{H}_{v_{z}u_{z}}\big{(}\cdots(\mathcal{H}_{v_{1}u_{1}}(\mathscr{G}))\big{)}$

without considering the color of the edges. This yields

$[F_{v_{z}u_{z}},\dots,[F_{v_{2}u_{2}},[\widehat{A},F_{v_{1}u_{1}}]]]=a^{*}F_{v^{*}u^{*}},$

where $a^{*}$ is the coefficient generated during the operation of the Lie
brackets. Therefore,
$F_{v^{*}u^{*}}\in\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}$, and, together with the connectivity of
$\mathscr{G}_{\\{F_{v^{*}u^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}}$, Lemma 15 now leads to

$\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm
LA}=\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\\{F_{v^{*}u^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n).$

By (11) it is obvious that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n).$

Step 2. In this step, we proceed to establish the result for the general case
by induction on the number of connected components of $\mathscr{G}_{\rm S}$.

Induction Hypothesis. If graph $\mathscr{G}_{\rm S}$ contains $m\geq 2$
connected components, then
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$. We will prove it for
$m+1$.

To do this, consider $\mathscr{G}_{\widetilde{A}}$. Let $v_{11}$ be the node
in $\mathrm{V}_{1}$ with $\deg(v_{11})=k>0$, and
$v_{i_{1}j_{1}},\dots,v_{i_{k}j_{k}}$ denote the nodes that are adjacent to
$v_{11}$. There is no loss of generality in assuming
$v_{i_{1}j_{1}},\dots,v_{i_{r}j_{r}}\in\mathrm{V}_{2}$, with $1\leq r\leq k$.
Analysis similar to that in the step 1 shows that

$\psi\big{(}[[\widetilde{A},F_{v_{11}v_{12}}],F_{v_{i_{1}j_{1}}v_{i_{*}j_{*}}}]\big{)}=\psi(\widehat{A})=\mathscr{G}_{a}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\mathscr{G}_{b},$

where $\mathscr{G}_{a}=({\rm
V}_{1}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
V}_{2},\mathscr{E}_{a})$ is a bi-graph with $|\mathscr{E}_{a}|\geq 1$,
$\mathscr{G}_{b}=(\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}_{i=3}^{m+1}{\rm
V}_{i},\mathscr{E}_{b})$ is a empty graph. Again we apply Lemma 13 to
$\mathscr{G}_{a}$ and obtain
$F_{v^{*}u^{*}}\in\big{\\{}\\{\widetilde{A}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}^{z}\big{\\}}_{\rm LA}$ with
$v^{*}\in\mathrm{V}_{1},u^{*}\in\mathrm{V}_{2}$. It follows immediately that
$F_{v^{*}u^{*}}\in\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}$. Now
$\mathscr{G}_{\\{F_{v^{*}u^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}}$ has $m$ connected components, and by our induction hypothesis,
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}\\{F_{v^{*}u^{*}}\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$. Thus we arrive at the
conclusion that
$\big{\\{}\\{A\\}\mathbin{\mathchoice{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\textstyle\bigcup$}}}{\vbox{\hbox{$\scriptstyle\bigcup$}}}{\vbox{\hbox{$\scriptscriptstyle\bigcup$}}}}{\rm
S}\big{\\}}_{\rm LA}=\operatorname{\mathfrak{su}}(n)$, and the proof is
complete. ∎

## References

  * [1] X. Wang, B. Li, J.-S. Li, I. R. Petersen, and G. Shi, “Controllability over Graphs for Bilinear Systems over Lie Groups,” in 59th IEEE CDC, December 2020.
  * [2] X. Wang, B. Li, J.-S. Li, I. R. Petersen, and G. Shi, “Controllability and Accessibility on Graphs for Bilinear Systems over Lie Groups,” arXiv:2007.11929.
  * [3] A. Jadbabaie, J. Lin, and A. S. Morse, “Coordination of groups of mobile autonomous agents using nearest neighbor rules,” IEEE Trans. Autom. Control, vol. 48, no. 6, pp. 988–1001, 2003.
  * [4] R. Olfati-Saber and R. Murray, “Consensus problems in the networks of agents with switching topology and time delays,” IEEE Trans. Automat. Control, vol.49, pp. 1520-1533, 2004.
  * [5] S. Martinez, J. Cortes, and F. Bullo, “Motion coordination with distributed information,” IEEE Control Systems Magazine, 27, pp. 75–88, 2007.
  * [6] A. Nedic, A. Ozdaglar, and P. A. Parrilo, “Constrained consensus and optimization in multi-agent networks,” IEEE Trans. Autom. Control, vol. 55, no. 4, pp. 922–938, Apr. 2010.
  * [7] M. Mesbahi and M. Egerstedt, Graph Theoretic Methods in Multiagent Networks. Princeton University Press, 2010.
  * [8] J. L. Gross and J. Yellen,  Graph Theory and Its Applications. Chapman $\&$ Hall/CRC, 2006.
  * [9] D. L. Elliott, Bilinear Control Systems: Matrices in Action. Dordrecht: Springer, 2009.
  * [10] D. L. Elliott and T. J. Tarn, “Controllability and observability for bilinear systems,” in  Proceedings of the SIAM National Meeting, 1971.
  * [11] R. W. Brockett, “Lie theory and control systems defined on spheres,” SIAM Journal on Applied Mathematics, vol. 25, no. 2, pp. 213–225, 1973.
  * [12] M. Aoki, “Variable structure systems with application to economics and biology,” in  Some Examples of Dynamic Bilinear Models in Economics. New York: Springer, 1975, pp. 163–169.
  * [13] C. Altafini, “Controllability of quantum mechanical systems by root space decomposition of su(N),” Journal of Mathematical Physics, vol. 43, no. 5, pp. 2051–2062, 2002.
  * [14] F. Albertini and D. D’Alessandro, “Notions of controllability for bilinear multilevel quantum systems,” IEEE Transactions on Automatic Control, vol. 48, no. 8, pp. 1399–1403, 2003.
  * [15] H. J. Sussmann and V. Jurdjevic, “Controllability of nonlinear systems,” Journal of Differential Equations, pp. 95–116, 1972.
  * [16] V. Jurdjevic and H. Sussmann, “Control systems on Lie groups,” Journal of Differential Equations, vol. 12, no. 2, pp. 313–329, 1972.
  * [17] W. Boothby, “A transitivity problemom control theory,” Journal of Differential Equations, vol. 17, pp. 296–307, 1975.
  * [18] G. Dirr and U. Helmke, “Lie theory for quantum control,” GAMM-Mitt, vol. 31, no. 1, pp. 59–93, 2008.
  * [19] A. Y. Khapalov and R. R. Mohler, “Reachable sets and controllability of bilinear timeinvarian systems: a qualitative approach,” IEEE Transactions on Automatic Control, vol. 41, no. 9, pp. 1342–1346, 1996.
  * [20] W. Zhang and J. S. Li, “Control of ensemble systems on special orthogonal groups,” Proc. 55th IEEE CDC, Las Vegas, USA, pp. 2529–2534, 2016.
  * [21] J. S. Li, W. Zhang, and L. Wang, “Computing controllability of systems on $\operatorname{SO}(n)$ over graphs,” Proc. 56th IEEE CDC, Melbourne, Australia, pp. 5511–5516, 2017.
  * [22] W. Zhang and J. S. Li, “Analyzing Controllability of Bilinear Systems on Symmetric Groups: Mapping Lie Brackets to Permutations,” IEEE Transactions on Automatic Control, DOI: 10.1109/TAC.2019.2963164.
  * [23] F. Albertini and D. D’Alessandro, “The Lie algebra structure and controllability of spin systems,” Linear Algebra and Its Applications, pp. 213–235, 2002.
  * [24] R. W. Brockett and N. Khaneja, “On the stochastic control of quantum ensembles,” in System Theory: Modeling, Analysis, and Control. Kluver Academic Publisher, Boston, pp. 75–96, 2000.
  * [25] T. Boukhobza and F. Hamelin, “Observability analysis for structured bilinear systems: A graph-theoretic approach,” Automatica, vol. 43, no. 11, pp. 1968–1974, 2007.
  * [26] Z. J. Ji and H. S. Yu, “A new perspective to graphical characterization of multiagent controllability,” IEEE transactions on cybernetics, vol. 47, no. 6, pp. 1471–1483, 2017.
  * [27] X. D. Chen, M. A. Belabbas, and T. Basar, “Controllability of formations over directed time-varying graphs,” IEEE Transactions on Control of Network Systems, vol. 4, no. 3, pp. 407–416, 2017.
  * [28] C. Aguilar and B. Gharesifard, “Graph controllability classes for the Laplacian leader-follower dynamics,” IEEE Transactions on Automatic Control, vol. 60, no. 6, pp. 1–13, 2015.
  * [29] A. Tsopelakos, M. A. Belabbas, and B. Gharesifard, “Classification of the structurally controllable zero-patterns for driftless bilinear control systems,” IEEE Transactions on Control of Network Systems, in press.
  * [30] A. Rahmani, M. Ji, M. Mesbahi, and M. Egerstedt, “Controllability of multi-agent systems from a graph-theoretic perspective,” SIAM Journal on Control and Optimization, vol. 48, no. 1, pp. 162–186, 2009.
  * [31] Z. Ji, Z. Wang, and H. Lin, “Interconnection topologies for multi-agent coordination under leader-follower framework,” Automatica, vol. 45, pp. 2857–2863, 2009.
  * [32] G. Parlangeli and G. Notarstefano, “On the reachability and observability of path and cycle graphs,” IEEE Transactions on Automatic Control, vol. 57, no. 3, pp. 743–748, 2012.
  * [33] B. Gharesifard, “Stabilization of bilinear sparse matrix control systems using periodic inputs,” Automatica, vol. 77, pp. 239-245, 2017.
  * [34] X. Chen, M. A. Belabbas, and T. Basar, “Controllability of formations over directed time-varying graphs,” IEEE Transactions on Control of Network Systems, vol. 4, no. 3, pp. 407–416, 2015.
  * [35] F. Liu and A. S. Morse, “A graphical characterization of structurally controllable linear systems with dependent parameters,” IEEE Transactions on Automatic Control, vol. 64, no. 11, pp. 4484–4495, 2019.
  * [36] J. Jia, H. L. Trentelman, W. Baar and M. K. Camlibel, “Strong structural controllability of systems on colored graphs,” IEEE Transactions on Automatic Control, vol. 65, no. 10, pp. 3977–3990, 2020.
  * [37] J. Jia, H. J. van Waarde, H. L. Trentelman and M. K. Camlibel, “A unifying framework for strong structural controllability,” IEEE Transactions on Automatic Control, vol. 66, no. 1, pp. 391–398, 2021.
  * [38] H. Qi, B. Mu, I. R. Petersen, and G. Shi, “Measurement-Induced Boolean Dynamics and Controllability for Quantum Networks,” Automatica, 114: 108816, 2020.