GMC Collisions as Triggers of Star Formation. III. Density and Magnetically Regulated Star Formation
Benjamin Wu Jonathan C. Tan Duncan Christie Fumitaka Nakamura Sven Van Loo David Collins







National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan
Department of Physics, University of Florida, Gainesville, FL 32611, USA
ben.wu@nao.ac.jp

Department of Physics, University of Florida, Gainesville, FL 32611, USA
Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan

School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK

Department of Physics, Florida State University, Tallahassee, FL 32306-4350, USA




We study giant molecular cloud (GMC) collisions and their ability to
trigger star cluster formation. We further develop our three
dimensional magnetized, turbulent, colliding GMC simulations by
implementing star formation sub-grid models. Two such models are
explored: (1) “Density-Regulated,” i.e., fixed efficiency per
free-fall time above a set density threshold; (2) “Magnetically-
Regulated,” i.e., fixed efficiency per free-fall time in regions that
are magnetically supercritical. Variations of parameters associated
with these models are also explored. In the non-colliding simulations,
the overall level of star formation is sensitive to model parameter
choices that relate to effective density thresholds. In the GMC
collision simulations, the final star formation rates and efficiencies
are relatively independent of these parameters. Between non-colliding
and colliding cases, we compare the morphologies of the resulting star
clusters, properties of star-forming gas, time evolution of the star
formation rate (SFR), spatial clustering of the stars, and resulting
kinematics of the stars in comparison to the natal gas. We find that
typical collisions, by creating larger amounts of dense gas, trigger
earlier and enhanced star formation, resulting in 10 times higher SFRs
and efficiencies. The star clusters formed from GMC collisions show
greater spatial sub-structure and more disturbed kinematics.









§ INTRODUCTION




Most stars are thought to form in clusters within giant molecular
clouds (GMCs). GMCs have typical hydrogen number densities of n_
  H=100  cm^-3, diameters of ∼tens of parsecs, masses
of up to 10^6 M_⊙, and average temperatures of ∼
10-30  K. Dense clumps within GMCs, potentially traced as,
e.g., Infrared Dark Clouds (IRDCs), are recognized as being the likely
precursors to star clusters
<cit.>. 
IRDCs have such high mass surface densities (Σ≳ 0.1 
  g cm^-2) that they are dark at mid-IR (∼10 μ m)
and even far-IR (∼70 μ m) <cit.>.  
Their low temperatures (10-20K) <cit.>,
high volume densities (n_ H>10^5  cm^-3), relatively
compact sizes (∼few pc), and masses (∼
10^2-10^5 M_⊙) indicate that they have the potential to be
the precursors to most of the observed mass range of star clusters
known in the Galaxy. The initial and early stages of star cluster
formation can also be traced by dust continuum emission 
<cit.> and by samples based on emission of dense gas
tracers <cit.>. Surveys of young embedded stars also
can probe the structure <cit.>, age distribution
<cit.> and kinematics of young clusters 
<cit.>.


Currently, the dominant processes that induce the collapse and
fragmentation of GMCs into star-forming clumps are
poorly understood. 
Various theoretical models include 
regulation by turbulence <cit.>,
regulation by magnetic fields <cit.>,
triggering by stellar feedback <cit.>, triggering by converging atomic
flows <cit.>, and triggering via
converging molecular flows, i.e., GMC-GMC collisions
<cit.>.


Semi-analytic models <cit.> and numerical simulations
<cit.> of global
galactic disks have shown that GMCs collide relatively frequently due
to the approximately 2D geometry of a thin disk and interaction rates
driven by differential rotation of galactic orbits. The average
timescale between GMC collisions was found to be about 20% of a local
orbital period within a flat rotation curve disk
<cit.> <cit.>.
A growing number of numerical studies have also
shown that collisions between molecular clouds can provide conditions
favorable for massive star and star cluster formation <cit.>.
We note that in general comparison of the results between the
  simulations of different groups is complicated by the use of
  different initial conditions, different numerical methods and
  different included physics.

Our approach here is to systematically build-up realism for our
GMC collision simulations by including additional physics step by
step that allows an understanding of the relative importance of
different input assumptions. <cit.> 
and <cit.>
developed a numerical study of GMC-GMC collisions, focusing on
understanding the physical mechanisms as well as using them to predict
observational diagnostics. Comparing magnetized, supersonically
turbulent GMCs in colliding and non-colliding cases over a wide
parameter space and investigating a varied array of potential
observational signatures, they found that a number of indicators
suggest similarities between the colliding scenarios to observed GMCs
and IRDCs. Further, dynamical virial analysis suggested that dense
^13CO-defined structures created through GMC collisions were more
likely to collapse and form massive star clusters when compared with
more quiescently evolving structures.

The next stage in our work is the crucial transition from collapsing
clumps into star clusters. Properties of the stars that form, along
with their dynamical evolution shortly thereafter, may provide insight
into the dominant star formation mechanisms.  The goal of this study
is to answer the question: do realistic models of GMC collisions
create star clusters that closely match the properties of observed
young star-forming regions?  We approach this question by further
building upon our previous numerical framework of GMCs through the
development of star formation sub-grid models, one of which is a novel
magnetically-regulated model. We combine our existing gas-focused
observational diagnostic methods with additional information from the
population of star particles. Thus, we hope to provide insight to the
star formation process by analyzing the evolution of IRDC-type
structures into young star clusters.

Section <ref> describes our numerical setup and the
various star formation models.  We then present our results in
Section <ref>, which include: 
gas and star cluster morphologies (<ref>), 
properties of star-forming gas (<ref>),
global star formation rates (SFRs) and efficiencies 
(<ref>), spatial clustering (<ref>), 
and star particle kinematics (<ref>). 
In Section <ref> we discuss our conclusions.




§ NUMERICAL MODEL






 §.§ Initial Conditions


We further develop the numerical framework described in Paper II and
introduce two star formation routines.  Our GMCs are identical to
those initialized in Paper II, which are motivated by observed GMC
properties. The clouds are self-gravitating, supersonically turbulent,
and magnetized. They are initialized with identical densities and
offset by an impact parameter. The clouds are embedded in an ambient
medium of ten times lower density (i.e., an atomic cold neutral
medium, CNM), which for the colliding case, is converging along with
the GMCs. The initial simulation properties are summarized in
Table <ref>.

The simulation domain is (128  pc)^3 and contains two
neighboring GMCs. The GMCs are initially uniform spheres, with
Hydrogen number densities of n_ H,GMC = 100  cm^-3 and
radii R_ GMC = 20.0 pc. This gives each GMC a mass M_ GMC
= 9.3 × 10^4 M_⊙.  The ambient gas represents the atomic
cold neutral medium (CNM) and has a density of n_ H,0 = 10 
  cm^-3. The centers of the GMCs are offset by 2R_ GMC in
the collision axis (x), 0 in the y-axis, and b=0.5R_ GMC in
the z-axis.

To approximate the density and velocity structures observed in GMCs,
our clouds are initialized with a supersonic turbulent velocity field
which is random, purely solenoidal, and follows the v_k^2∝
k^-4 relation, where k=π/d is the wavenumber for an eddy
diameter d. Conventionally, the “k-mode” is normalized to the
simulation box length. The gas within the GMC is initialized with Mach
number ℳ_s≡σ/ c_s= 23 (for T=15 K
conditions), of order virial. We set our fiducial k-modes to be
{2,...,20}, where each mode within this range is
excited.  This is representative of the large-scale turbulent
velocities (small k) spanning from the GMC diameters down to a small
enough minimum scale (large k), which is numerically resolved, but
expected to cascade to smaller scales. We do not drive turbulence,
instead letting it decay within a few dynamical times.  Note also that
turbulence is initialized only within the initial volume of the GMCs
while we leave the ambient medium non-turbulent. Note also the GMC
collision will also drive turbulence in the clouds in that case.

A large-scale uniform magnetic field of strength 10 μ G is
initialized throughout the box at an angle θ=60^∘ with
respect to the collision (x-) axis. This choice of |B| is
motivated by the Zeeman measurements of typical GMC field strengths,
summarized by <cit.>.

In the fiducial colliding case, the bulk flows (including both the
ambient gas and the GMCs) have a relative velocity of v_
  rel=10  km s^-1. In the non-colliding case, there is no
bulk velocity flow.

The simulations are run for 5  Myr to investigate the onset of
star formation. Note that this is 1  Myr longer than the
simulations described in Paper II, which focused on gas properties of
the pre-star-forming clump. Note also that the freefall time given the
initial uniform density GMCs is t_
  ff=(3π/[32Gρ])^1/2≃ 4.35 Myr. However, the values of
t_ ff for the denser substructures created by turbulence and by
the collision are much shorter.
Star formation is expected to occur in both non-colliding and 
colliding cases, with the detailed properties of resulting star 
clusters acting as the key point of our investigation.


llllclc
Initial Simulation Properties 
4
0pt

                       
GMC          Ambient 


n_ H            (cm^-3)         100                       10      

R                    (pc)              20                         ...    

M                    (M_⊙)             9.3×10^4          ...    

T                    (K)                       15                        150     

t_ ff           (Myr)                     4.35                       ...    

c_ s            (km/s)                    0.23                      0.72    

v_A                (km/s)                    1.84                      5.83    

v_ vir          (km/s)                    4.9                        ...    

σ               (km/s)                    5.2                        ...    

ℳ_ s            ...             23                         ...    
 
ℳ_A               ...              2.82                       ...    
 
k-mode               (k_1,k_2)           (2,20)                     ...    

v_ bulk         (km/s)                    ±5                    ±5  

B                    (μ G)           10                        10      

λa          ...          4.3                       1.5     

βb            ...          0.015                     0.015   

anormalized mass-to-flux ratio: λ = (M/Φ)/(1/2π G^1/2)
bthermal-to-magnetic pressure ratio: β=8π c_s^2 ρ_0 / B^2










 §.§ Numerical Code



Our models are run using
[http://enzo-project.org (v2.4)], a
magnetohydrodynamics (MHD) adaptive mesh refinement (AMR) code
<cit.>. We use the Dedner-MHD method, which solves the
solves the MHD equations using the Harten-Lax-van Leer with
Discontinuities (HLLD) method and a piecewise linear reconstruction
method (PLM). The time is evolved using the MUSCL 2nd-order
Runge-Kutta method. The ∇·B=0 solenoidal constraint
of the magnetic field is maintained via a hyperbolic divergence
cleaning method <cit.>.

The simulation domain is realized with a top level root grid of
128^3 with 3 additional levels of AMR. Our models thus have an
effective resolution of 1024^3, with a minimum grid cell size of
0.125 pc. We refine solely on the local Jeans length, setting a
necessary requirement of resolving by 8 cells, i.e., higher than the 4
cells typically used to avoid artificial fragmentation
<cit.>). Our higher resolution leads to larger
volumes of the GMCs being better resolved and thus generally better
resolution of, e.g., shocks <cit.>.
We note that the Jeans criterion assumes purely thermal support. The
gas in our simulations also has some magnetic support, so its
effective “magneto-Jeans length” will be significantly larger than
the thermal Jeans length in directions perpendicular to the magnetic
field.

Due to the relatively high bulk velocities and 
potentially strong magnetic fields, we require the use of the
“dual energy formalism”
<cit.> , which separately solves the internal
energy equation as well as the total energy equation, ensuring
  accurate calculation of pressures and temperatures in these
  conditions. If the ratio of thermal to total energy is less than
0.001, then the temperature is calculated from the internal
pressure. Otherwise, the total energy is used.

Additionally, we employ the “Alfvén limiter” (described in Paper
II) to avoid exceedingly small timesteps set by Alfvén waves. This
acts by choosing a maximum Alfvén velocity, v_A,
    max=B/√(4πρ_ min)=1×10^7  cm
   s^-1, and setting a density floor that is determined by the
magnetic field. This predominantly affects only small pockets of very
low-density gas with which we are less interested, and thus the
dynamical results are deemed unaffected by this limiter.




 §.§ Thermal Processes



We assume a constant mean particle mass (μ=2.33 m_ H)
throughout the simulation domain for simplicity, as our focus is on
the dense molecular gas of GMCs.  We also choose a constant adiabatic
index γ=5/3. Note that this essentially ignores certain
excitation modes of H_2 that may be relevant (i.e., shocks),
but it is still the most appropriate single-valued choice of γ,
given our focus on the dynamics of cold H_2. Also, we assume
n_ He=0.1n_ H, giving a mass per H of
2.34×10^-24  g.

The PDR-based heating and cooling functions developed in Paper I are
again used in these simulations. The assumptions are: (1) FUV
radiation field of G_0=4 (i.e., appropriate conditions for the
inner Galaxy, e.g., at Galactocentric distances of ∼ 4 kpc) and
(2) a background cosmic ray ionization rate of ζ = 1.0 ×
10^-16  s^-1.  The heating/cooling functions trace the
atomic to molecular transition and recreate a multi-phase ISM. They
span density and temperature ranges of 10^-3≤ n_ H/
  cm^-3≤ 10^6 (extended to 10^10  cm^-3 via
extrapolation) and 2.7 ≤ T/ K≤ 10^7, respectively.

We use the  external chemistry and cooling
library[https://grackle.readthedocs.org/]
<cit.> to incorporate our heating/cooling
functions in tabular form into , modifying the energy
equation.

In order to avoid numerical instabilities related to the
heating/cooling processes, we limit the timestep on each AMR level to
a factor of 0.2 the minimum cooling time. Additionally, we set a hard
floor for the minimum cooling timestep of 625  yrs.






 §.§ Star Formation



We utilize the particle machinery of  to model star
formation.  Specifically, star particles (i.e., collisionless, point
particles with mass m_⋆) form within a simulation cell if
certain local criteria are met.  Two star formation routines are
developed: (1) density-regulated star formation; (2)
magnetically-regulated star formation.



  §.§.§ Density-regulated star formation



Our first star formation routine is a “density-regulated” model,
based on that of <cit.> <cit.>.
Stars are formed within a cell only if they have been
refined to the finest level of resolution and the density exceeds a
particular threshold value, n_ H,sf. The fiducial star
formation density threshold is chosen to be n_ H,sf =
10^6  cm^-3, which is set partly based on observed densities
of pre-stellar cores <cit.>.
We will consider variation of this parameter by a factor of two to
higher and lower values. The temperature in the cell is also required
to be <3000 K, to avoid star formation in dense, shock-heated
regions, but we will see that this constraint is not of practical
concern for the simulations presented here. Note that there is no
requirement for gravitational boundedness of gas in the cell. Nor is
there a requirement for net convergence of gas flow to the cell. These
choices are motivated by the fact such conditions are not well
resolved on the local scales associated with an individual cell. In
addition, we expect that processes such as turbulence and diffusion of
magnetic flux that are occurring on sub-grid scales (or scales near the
grid scale that are not well resolved) will regulate star formation,
e.g., creating local conditions that are gravitationally unstable,
perhaps via converging flows. With these points in mind, this star
formation sub-grid model using a density threshold is thus designed to
be as simple as possible, enabling us to gain a clear understanding of
how the results depend on its input parameters.

In cells meeting the above conditions, star particles are then
produced so that the SFR is, on average, equal to that expected if
there is a fixed star formation efficiency per local free-fall time,
ϵ_ ff, where the local free-fall time, t_ ff, is
expressed as 

    t_ ff    =    (3π/32Gρ)^1/2
        =     4.4×10^4 n_ H,6^-1/2  yr,

i.e., the value for collapse of a uniform density sphere, where
n_ H,6≡ n_ H/10^6  cm^-3. We adopt a fiducial
choice of ϵ_ ff=0.02, motivated by observations of GMCs,
their star-forming clumps and stellar populations in embedded
clusters, which suggest fairly low and density-independent values of
ϵ_ ff <cit.>. 
Thus the SFR is

    ṁ_⋆    =    ϵ_ ffm_ gas/t_ ff
        =     2.9×10^-5(ϵ_ ff/0.02) (Δ x/0.125  pc)^3 n_ H,6^3/2 M_⊙  yr^-1,

where we have normalized to the minimum cell size, Δ x,
relevant for the simulations in this paper.

The timesteps in the simulation are typically quite short, i.e., much
less than the signal crossing time of a cell, i.e., ≪
1.2×10^5yr for a signal speed of 1 km s^-1. Thus the
average mass of stars that are expected to be created in a given cell
in a given simulation timestep is often very small, i.e.,
<1 M_⊙.

To enable both a practical computation that does not involve too many
star particles, but also with the eventual aim of producing star
particles with masses that are characteristic of observed stellar
masses, the star formation sub-grid model also involves a parameter of
a minimum star particle mass, m_⋆,min. For the
density-regulated models we consider here, we set m_⋆,min=10 M_⊙. Thus in this case the star particles
represent small (sub-)clusters of stars, since the mean stellar mass
is ≲ 1 M_⊙ for realistic stellar initial mass functions
<cit.>
(see also discussion of star particle
dynamics in <ref>). With this value of m_⋆,min
we are almost always in a regime in which the mass of stars to be
created in a given timestep is smaller than m_⋆,min and so
the decision to form a star particle or not needs to be implemented
probabilistically, i.e., the “stochastic star formation” regime. In
this case, the star particle is formed with probability ṁ_⋆Δ t / m_⋆,min, where Δ t is the simulation
timestep. If on the other hand ṁ_⋆Δ t > m_⋆,min (which can occur in certain circumstances), then the
star particle is simply created with this mass.

Another factor affecting the choice of m_⋆,min is the
desire not to change the gas mass in a cell by too large a fraction
when the star particle is created, i.e., to avoid too large changes in
density, pressure, etc. In general we set an upper limit of this
fraction of 0.5. In the fiducial case a cell of size 0.125 pc at the
star formation threshold density contains a minimum gas mass of
m_ gas,min=63 M_⊙, so this fraction is ≲ 0.17 for
these models (≲ 0.34 for the lower threshold density case).

Overall there are three density (“d”)-regulated runs (i.e., three
choices of threshold density) for each of the noncolliding (“nocol”)
and colliding (“col”) simulation set-ups. The parameters of these
star formation models and simulations are listed in
Table <ref>.






  §.§.§ Magnetically-regulated star formation


We introduce a new “magnetically-regulated” star formation model
that takes into account magnetic criticality, i.e., star formation is
only allowed to proceed if a cell has a mass to flux ratio that is
greater than a certain value. If the cell is magnetically
“supercritical” by this criterion, then it forms stars at a fixed
efficiency per local free-fall time, ϵ_ ff, where we
will adopt the same value of 0.02 that was used in the
density-regulated models. Thus this magnetic criticality condition
acts to replace the density-threshold criterion of
<ref>. However, as we discuss below, the choice of m_⋆,min also introduces an effective minimum density for star
formation in this model also.

To assess the mass-to-flux ratio criterion, as an approximation, we
treat each grid cell individually and calculate the dimensionless
mass-to-flux ratio

    μ_ cell = ρΔ x √(G)/B c_1,

where ρ is the density within a cell of length Δ x, G is
the gravitational constant, B is the strength of the magnetic field
within the cell, and c_1 comes from defining the critical
mass-to-flux ratio as

    ( M/Φ)_ crit = c_1/√(G)

and is ultimately dependent of the geometry of the system. For an
infinite disk, the value is c_1=1/(2π)
<cit.>; for an isolated cloud it is roughly
1/√(63)∼ 0.126 <cit.>. We will
consider variations of c_1 of a factor of two to higher and lower
values. These and other parameters of the magnetically
(“B”)-regulated star formation models are also listed in
Table <ref>. We note that although the true mass-to-flux
ratio depends on the geometry of the entire flux tube and cannot be
completely confined to a localized quantity, this only acts as a
first-order correction.



If a cell is magnetically subcritical (i.e., μ_ cell<1), the
magnetic pressure is deemed strong enough to withstand gravitational
contraction, preventing any stars from forming within that cell. For
those cells that are magnetically supercritical (i.e., μ_
  cell>1), then star formation may be allowed to occur. However the
two other criteria introduced in <ref>, i.e., the cell is
resolved at the finest refinement level and T<3000 K, also must be
satisfied. In addition, we also limit the fraction of gas mass in a
cell that is turned into stars at a single timestep to be
<0.5. Given the minimum star particle mass, m_⋆,min,
this imposes an “effective density threshold” for the
magnetically-regulated star formation model. For this reason, we will
investigate magnetically-regulated models with m_⋆,min =
10 M_⊙ (as in the density-regulated models), but also with
m_⋆,min = 1 M_⊙. These choices correspond to
effective minimum threshold densities of n_ H,sf= 3.55×
10^5  cm^-3 for m_⋆,  min=10 M_⊙ and
n_ H,sf= 3.55× 10^4  cm^-3 for m_⋆, 
    min=1 M_⊙.

If all the above conditions are satisfied, then the star formation
process is allowed to occur at fixed efficiency per local free-fall
time, as described in <ref>. We will see that the
magnetically-regulated models with m_⋆,  min=1 M_⊙
can form significant numbers of star particles out of the stochastic
regime, but these masses should not be interpreted as being a
realistic assessment of the stellar initial mass function, since their
values depend on the size of the simulation timestep.

Overall there are four magnetically (“B”)-regulated runs (i.e.,
three choices of mass-to-flux threshold for m_⋆, 
    min=10 M_⊙ and one run with m_⋆, 
    min=1 M_⊙ at the fiducial mass-to-flux threshold) for
each of the noncolliding (“nocol”) and colliding (“col”)
simulation set-ups. The parameters of these star formation models and
simulations are also listed in Table <ref>.





  §.§.§ Star particle and star cluster dynamics


Once the star formation criteria are met, mass is removed from the
cell and placed into a point-like star particle. These evolve as a
collisionless N-body system. However, these are not treated as
accreting sink particles, so they do not gain additional mass from the
gas, which we expect to be realistic due to the action of stellar
winds from the young stars.  The particles still interact with the gas
gravitationally via a cloud-in-cell (CIC) 
algorithm which maps the particle positions onto the grid. This limits
the closest distances between star particles to the grid resolution,
ultimately resulting in softer mutual gravitational interactions.  As
a result, small scale, i.e., internal, star cluster dynamics is not
expected to be well-modeled.  However, the early stages and larger
scales of the spatial and kinematic distribution of the stars, should
be more accurately followed.

Note also that our ability to follow the true internal dynamics of the
formed star clusters is limited by the fact that we do not fully allow
for the presence of a range of stellar masses, including both low-mass
and high-mass stars, or the presence of binary or higher order
multiple star systems. However, since our ability to accurately follow
the dynamical evolution of the star cluster is mostly limited by the
fact that gravitational forces are not well resolved below the grid
scale of the simulation, our focus is mostly on the global
distribution of stars in the simulations and the large scale spatial
and kinematic distributions of the stars in the clusters, e.g., low
order spatial mode asymmetries.

The current modeling also does not include feedback from the formed
star particles. A goal of a future paper is to include protostellar
outflow momentum feedback in these models, but at the moment the star
formation that results should be considered a baseline estimate in the
limit of zero feedback.



§ RESULTS



We perform analysis of each the simulations, comparing and contrasting
star formation models as well as non-colliding vs. colliding cases. In
particular, we discuss: morphology of the clouds and clusters
(<ref>); properties of star-forming gas
(<ref>); global star formation rates
(<ref>); spatial clustering of stars
(<ref>); and star vs. gas kinematics
(<ref>).



 §.§ Cloud and Cluster Morphologies












The morphology of the gas and the stars are shown in
Figure <ref> for the non-colliding clouds and
Figure <ref> for the colliding clouds. In the
non-colliding cases, the gas evolution is essentially identical, where
turbulent velocities and self-gravity create a network of relatively
slowly growing filaments with increasing differentiation in mass
surface density. Evolution is relatively passive and quiescent.  In
general, the onset of star formation takes place near 4.5  Myr
for the m_⋆,  min=10 M_⊙ cases and 3.5  Myr
for the m_⋆,  min=1 M_⊙ case.

Both density-regulated and magnetically-regulated models result in 
pockets of localized star formation concentrated at density peaks 
within filamentary structures. The slightly more populated network of 
filaments in the northeast region forms a higher number of 
(small) star clusters, but overall, star formation is scattered 
sparsely throughout both GMCs and remains relatively spatially 
isolated.  The magnetically-regulated models form clusters with a 
higher degree of elongation, i.e., following the axes of the natal 
filaments. There is also slightly more widespread star formation 
activity compared with the density regulated case. By t=5.0  Myr, 
approximately 5-8 clusters have formed in the density-regulated cases, 
with the higher critical density models forming fewer clusters, whereas 
10-12 separate clusters have formed in the magnetically-regulated cases.

Differences are more pronounced in the m_⋆, 
    min=1 M_⊙ case, as stars form in elongated clusters along
the filaments instead of the more localized spherical clusters as in
the density regulated case or the slightly eccentric clusters as in
the m_⋆,  min=10 M_⊙ magnetically regulated
cases. 

However, one needs to bear in mind that the magnetically-regulated
models also involve an effective minimum density threshold as an
additional requirement for star formation. This threshold density
depends on the minimum star particle mass that is allowed in the model
via the requirement that no more than 50% of the cell's gas mass can
be converted to a star particle (see <ref>). These effective
threshold densities are n_ H,⋆ min= 3.55× 10^5 
  cm^-3 for m_⋆,  min=10 M_⊙ and n_
  H,⋆ min= 3.55× 10^4  cm^-3 for m_⋆, 
    min=1 M_⊙. Thus the variation in m_⋆,  min
is a way of investigating how varying this effective density threshold
influences the resulting stellar population. Recall that for
star-forming gas, the star formation activity in lower density regions
is suppressed because the rate scales inversely with the local
free-fall time, i.e., SFR =ϵ_ ff m_ cell/t_ ff∝ n_ H^3/2. Thus the overall SFR in these
magnetically-regulated models will depend on both the probability
distribution function (PDF) of densities of the gas above the
effective threshold density that achieves the magnetic criticality
condition. In the simple density-regulated models, the SFR will simply
depend on the PDF of densities above the threshold density.



Considering now the GMC collision cases (Figure <ref>),
we see that they produce a much more active and dynamic environment
that leads to creation of much denser gas structures (see also Paper
II). At the interface of the colliding flows, a primary high-density
filamentary structure is formed.  This relatively compact, sheet-like
structure lies predominantly in the plane perpendicular to the
collision axis, with smaller filaments extending outward in various
directions. Mass surface densities of ∼1  g cm^-3 are
reached much sooner compared to the non-colliding case.  This results
in earlier and more rapid star formation, generally beginning near
3.0 to 3.5  Myr for the m_⋆,  min=10 M_⊙
cases and earlier than 3.0  Myr for the m_⋆, 
    min=1 M_⊙ case of magnetically-regulated star formation.

In all such cases, the clusters form in the central colliding region,
at the peaks of filaments located in the primary filamentary
network. These sites often correspond with overdense clumps located at
filament junctions, potentially pointing toward star formation
triggered by filament-filament interactions on the smaller scale. By
t=4.5  Myr, the individual star clusters have grown and merged
into one dominant star cluster located near
(x,y)=(5 pc,10 pc), while stars continue to form from
dense clumps scattered throughout the post-shock colliding region.
This large cluster appears to contain multiple populations of smaller
star clusters that have merged together through a combination of
gravitational attraction and initial velocity inherited from the natal
gas of the collision.  The spatially separated clusters from earlier
times have grown in population and are moving toward the main cluster,
while a few smaller clusters are continuing to form along the
still-colliding dense filamentary gas. By t=5.0 Myr, the main
cluster (which has grown to a few thousand stars in the m_⋆,
   min=10 M_⊙ case and a factor of 10 higher in the
1 M_⊙ case) is co-located with the majority of the dense gas,
as more star clusters form in the vicinity. There exists a small
population of individual stars that
form in relative isolation and/or are dynamically ejected from the
denser regions.

The factor of two variations in n_ H,sf do not greatly alter
the overall cluster morphology. However, there are small differences
in total cluster number as well as cluster size corresponding to the
density threshold, with increasing thresholds leading to reduced star
and cluster formation. 
The magnetically-regulated models exhibit slightly earlier star
formation, initializing just prior to t=3.0 Myr in each case, and a
higher number of clusters formed, which culminates in a larger central
cluster at late times compared with the density-regulated models.
Within these models, increasing values of c_1 result in reduced
star formation overall, though the locations where star formation is
centered do not change.

The  model initiates star formation the earliest,
with a primary central cluster and 5-8 smaller clusters already formed
by t=3.0  Myr.  Stars form in elongated structures directly
corresponding to the dense gas filaments similar, on small scales, to
that of the  model.  By t=4.0  Myr, stars
are present throughout the primary filament, still generally following
the filamentary structure of the gas, with smaller clusters forming
elsewhere throughout the colliding region. By t=5.0  Myr, the
primary central cluster has grown directly as well as from
gravitational interactions with the nearby clusters. Outlying clusters
have continued to increase in size and number.


Within dense filaments, the B-field is generally aligned perpendicular
to the filament axis (see Paper II). Qualitatively, the mass-to-flux
is expected to be highest at the density peaks, locally, but is
expected to decrease when the entire flux tube is taken into account
due to the lower-density environment surrounding the filaments. 
We note also that although these are ideal MHD simulations, some
numerical diffusion of flux is expected to occur that may influence
the star formation activity. The effects of modeling non-ideal MHD
processes will be explored in a future paper in this series.


Figure <ref> shows the combined 3D structure of the gas 
density, magnetic field geometry, and star particles. The 
and  models are compared at the same time, t=5.0  Myr, 
revealing the denser and more compact structure created in the GMC collision. 
These figures show the contrasting global morphology of the gas and stellar structure. Detailed analysis of various aspects of these properties is performed
in subsequent sections.





The different star formation models 
can be directly compared via |B| vs. n_ H phaseplots
(Figure <ref>). 
The critical thresholds for each model are plotted over phaseplots in 
respective non-colliding and colliding runs without star formation
(i.e., fiducial runs from Paper II extended to t=5.0  Myr). 
In this manner, the total mass along with various properties of the 
gas affected by each star formation model can be estimated.
The colliding case forms regions of overall higher density and
magnetization, both enhanced by approximately an order of magnitude.
At a given density, gas in the colliding case generally contains
stronger field strengths due to the nature of the compressive
flows. The star formation thresholds for both the density-regulated
and magnetically-regulated star formation routines are overplotted as
blue and red lines, respectively.

The density-regulated star formation regime affects a greater total
gas mass in the colliding case. As the critical density threshold
decreases, the number of affected cells in both scenarios increases,
leading to increased star formation regardless of magnetic field
strength.

As the threshold for mass-to-flux ratio is lowered, a similar pattern
of increasing star formation occurs. 
Key differences from the density-regulated models become
apparent as star formation is now allowed to occur in regimes of low 
density, super-critical gas and is inhibited in high-density, 
sub-critical gas. Overall, the various models primarily create stars 
from the same gas, though narrow regimes exist in which stars form 
exclusively within certain routines.
As discussed above, in these
magnetically-regulated models, the effective minimum density
thresholds for star formation provide an additional bound.  In the
m_⋆,  min=10 M_⊙ models, much of the gas in the
colliding cases–and even more so in the non-colliding cases–is
limited by this effective density threshold. 
The 1 M_⊙ cases allow star formation from a larger amount of
locally super-critical gas. However, we will see below that the
overall mass of stars formed by 5 Myr in the colliding case depends
only weakly on this choice.





 §.§ Properties of Star-forming Gas






We examine the masses of young stars and properties of their
progenitor gas cells, just before a star particle is
created. Figure <ref> displays the cumulative histograms
over 5 Myr of the stellar masses and key properties of the
star-forming gas.

For the non-colliding cases, the stars form strictly at their
threshold masses, indicating purely stochastic star formation is
occurring. The density-regulated models form approximately 100-250
stars each, with higher critical density thresholds resulting in fewer
stars.  The distributions of cell densities peak at the thresholds of
0.5, 1, and 2 × 10^6  cm^-3 for the respective models
and extend above the cutoffs by factors of a few. Gas temperatures
range from 6 to 40 K, averaging approximately 10 to 20 K with higher
density thresholds resulting in slightly lower temperatures. The local
normalized mass-to-flux ratio of the star-forming cells in these
models is supercritical by factors of few for  to
few tens for . The velocities of these cells are
generally a few km s^-1, consistent with values expected
from decaying turbulence in the self-gravitating GMCs.

The magnetically-regulated non-colliding models exhibit slightly
higher numbers of stars (few hundred for the m_⋆, 
    min=10 M_⊙ models).  Across these three models,
distributions for density have peaks at the cutoff of ∼
3×10^5  cm^-3, temperatures primarily near 20 K
(slightly above equilibrium), μ_ cell near 3-5 (slightly
supercritical), and |v| near 2 km s^-1 (supersonic but
consistent with decay of the initial turbulence). There exist slight
trends of increasing density and decreasing temperature as the
thresholds for μ_ cell increase between models. For the
m_⋆,  min=1 M_⊙ model, approximately 20 times more
star particles are created. 
The distribution of cell densities also exhibit a cutoff at the
effective minimum density of ∼ 3×10^4  cm^-3 in
this model, with the spread of densities reaching up to a factor of 20
higher. The temperature and velocity distributions exhibit similar
peaks and spreads as their higher minimum stellar mass
counterparts. However, the cells generally have a lower μ_
  cell near 1 (magnetically critical) when stars are formed. This
suggests that the condition for criticality is reached before the cell
density grows to a point at which it can produce more massive stars in
a given timestep and so stays in the stochastic limit.

The colliding models with density-regulated star formation produce
∼ 6×10^3 stars, over 20 times the number formed from the
non-colliding models over the same time period. The density
distributions are similar to the non-colliding cases, peaking at the
cutoffs, but exhibit an increased spread with cells reaching densities
higher by factors of a few. The gas temperatures are also higher,
averaging near 30-40 K but with a few cells reaching ∼ 90 
  K. Temperatures are generally lower for the higher-density cutoff
models, but all peak at temperatures higher than equilibrium, likely
due to shocks produced throughout the primary colliding region. The
collision also produces high-density gas at a wide range of μ_
  cell, ranging from a few times subcritical up to ∼ 20 times
supercritical. The distributions peak near the value for magnetic
criticality, with lower higher-density cutoff models corresponding
with higher values of μ_ cell. 
Cell velocity distributions
are nearly identical, peaking near 10-20 km s^-1.


For the magnetically-regulated colliding models with
m_⋆, min=10 M_⊙, a similar star particle count
is seen, exceeding their respective non-colliding models by factors of
10 to 20. The  model also forms some stars outside of
the stochastic regime, as masses of ∼ 12-13M_⊙ are
created. It is important to recall that the expected mass of the star
particle to be created depends on the local SFR in the cell, i.e., on
ϵ_ ff and the cell density, but also on the timestep of
the simulation. Thus the presence and mass distribution of these
higher mass star particles should not be over-interpreted. The
presence of stars outside the stochastic regime simply indicates that
some very high density, high SFR cells are present, and this is
confirmed in the plots showing the density distributions, with some
densities up to n_ H=10^7  cm^-3. We note that
star-forming cells can also have higher temperatures near 30 to 40 K,
perhaps indicating creation of the dense gas in shocks, but recall
that the star formation sub-grid model does not assess the degree of
gravitational instability in the gas. The star-forming cells show a
concentration of gas at the minimum magnetic criticality cutoff.
They also exhibit generally higher velocities indicating strong 
turbulence and/or bulk motion associated with the GMCs. 


The  model has the greatest total number of stars
formed, i.e., ∼ 5×10^4 and forms a range of stellar masses
up to ∼ 7 M_⊙ (but again this mass function should not be
expected to be compared to a real IMF, rather being simply the way the
model ensures the total mass of stars formed is correct given the
model parameters). 
Cell number densities range from ∼ 3×10^4 to ∼
9×10^6  cm^-3. Temperatures are near 40 K and μ_
  cell reaches a few tens but increases in cell number towards the
critical value cutoff of 1. Velocities exhibit a similar trend as the
10 M_⊙ models, showing high levels of turbulence, bulk
motion,  and/or infall to the primary cluster.




 §.§ Star Formation Rates and Efficiencies





The star formation rate (SFR) and overall star formation efficiency
(SFE) of molecular clouds are important quantities that help determine
the global galactic star formation process.  The time evolution of
these quantities for both non-colliding and colliding cases, for each
star formation model, is shown in Figure <ref>.

The SFR is calculated as the time derivative of the total mass of the
star particles. The efficiencies are determined by normalizing the
total stellar mass by the combined gas mass of the two original
GMCs. The evolution of these quantities is measured in simulation time,
as well as relative to the freefall time of the initial
GMC density (t_ ff=4.35 Myr).

In the non-colliding density-regulated cases, star formation initiates
shortly after 1 t_ ff, i.e., at approximately t=4.4 Myr. The
higher-cutoff density models form stars at slightly later times
corresponding to when the critical density is achieved. At a given
time, SFR and SFE vary by factors of a few between models. Over the 
course of the next 0.5 Myr (until simulation completion), the SFRs 
increase to ∼ (3-6)×10^-3 M_⊙  yr^-1 and 
then generally level off. Total stellar masses of ∼
9×10^2-3×10^3 M_⊙ are created, corresponding to
ϵ≈ 1% by 1.15 t_ ff.

For the magnetically-regulated cases with m_⋆,min=10 M_⊙, star formation also starts after 1 t_
  ff, and evolves in a similar manner as the density-regulated
models except with slightly higher SFRs and efficiencies. The
differences between these models is also much smaller, as the three
10 M_⊙ magnetically-regulated models reach about
6×10^-3 M_⊙  yr^-1 by 1.1 t_ ff and
then level off. The SFEs also reach and slightly exceed 1% by the
simulation end time. The  case exhibits the most
dissimilar behavior of the non-colliding models, initiating star
formation approximately 1 Myr earlier and reaching
4×10^-3 M_⊙  yr^-1 and ϵ=1% by 1
t_ ff. By 5 Myr, the SFR and SFR are approximately 2-3 times
greater than the other magnetically-regulated models.

The above trends are mostly likely caused by the fact that all these
star formation models have effective density thresholds that need to
be met to allow star formation to proceed, even the
magnetically-regulated models (see <ref>). These thresholds
decrease monotonically as we consider the density-regulated models,
then the B-(0.5,1,2)-nocol models, and then the B-1-1M-nocol
model. The simulations are in a regime in which the total SFR and
eventual total SFE are set mostly by the fraction of gas in the GMCs
that can meet these density threshold criteria. For the particular
B-field strengths in these simulations (i.e., 10μ G),
the choice of the magnetic threshold parameter does not play a
significant role in setting the SFR.

The colliding cases produce much higher SFRs and SFEs during their
evolution. The density-regulated models begin forming stars 
at a rapid pace shortly after t=3.2 Myr, with higher density 
thresholds slightly delaying the onset of star formation. There is some
oscillation in the growth of the SFRs, but overall it increases from
onset until t_ ff near 0.08 M_⊙  yr^-1 and 
then levels off through the culmination of the simulations.
Star formation efficiencies reach 1% by 3.7 to 4 Myr and more than 20% by
t_ ff. While the early behavior differs slightly in time between
the density-regulated models, they appear to converge at later times.

The m_⋆,min=10 M_⊙ magnetically-regulated colliding
cases show very similar results to each other throughout the whole
evolution, which indicates that the SFR is not limited by the
mass-to-flux thresholds in this simulation set-up. Indeed, these
models also converge with the density-regulated models by about
4 Myr. 
 starts forming stars at the earliest
times, but it also shows convergence in SFR by about 4 Myr. 
These results indicate that the SFR is in fact not limited by the density
threshold criteria either. In the GMC-GMC collision the SFRs appear to
be set by the creation of structures that can place gas at densities
greater than any of the threshold densities, after which, even with
ϵ_ ff=0.02, it is turned quite efficiently into stars.


It is important to note that stellar feedback has not yet been
included in our star formation models. Our current treatment may be a
good approximation for initial 
SFRs, but mechanisms such as protostellar outflows that become 
important during formation, and ionization, winds and radiation pressure 
from massive stars soon after, will likely result in reduced SFRs.




 §.§ Spatial Clustering



We investigate various quantitative metrics for spatial structure
of the star clusters formed in our simulations. Global star and gas 
properties of the primary clusters are measured and the
angular dispersion parameter (ADP) and minimum
spanning tree (MST) methods are used to analyze cluster substructure.
The ADP is sensitive to angular substructure at chosen radii, while 
the MST determines the degree of overall centrally concentrated 
clustering.

In order to define the primary cluster within a given model, we use
the Density-based spatial clustering of applications with
noise (DBSCAN) algorithm <cit.>. This density-based
clustering algorithm is applied to our projected star particle data
and the median particle position of the highest population cluster is
used as our initial cluster center. A circular aperture with initial
radius of 0.4 pc is centered at this point, and a new center is 
determined by finding the center of mass using stars included only 
within this aperture. This process is repeated with aperture radii 
iteratively decreasing by factors of two, down to length scales of
0.1 pc. The ADP is found for the primary cluster using each of these
defined centers, while the MST is found for the entire domain.






  §.§.§ Global structure of primary cluster







We measure global structural properties of the primary star clusters
created in the simulations. Figure <ref> shows results for
these clusters that have formed by t=5  Myr in four models:
, , , and
.  Due to relatively sparse particle density, the
clusters in the non-colliding m_⋆,min=10 M_⊙ models
were not included in this analysis. Also note that while elliptical
annuli are displayed (see ADP discussion in Section <ref>),
global properties for cluster structure are calculated using circular
annuli.

The top row of Figure <ref> shows the mass surface density of
stars locally within each annulus, Σ_⋆, the enclosed
average mass surface density of stars, Σ_⋆,
and the enclosed average mass surface density of the gas,
Σ_ gas.  For the non-colliding case,
Σ_⋆ ranges from about 10^3 to 5×
10^4  M_⊙ pc^-2, whereas the colliding cases form
clusters with approximately 10^5 to 10^6 
  M_⊙ pc^-2. It can be seen that Σ_⋆ falls off
quickly, reaching ∼10^2 and ∼10^4 
  M_⊙ pc^-2 by R=0.4  pc for non-colliding and
colliding cases, respectively. We find best fit power law profiles of
the form:

    Σ_*(R) = A ( R/0.2  pc)^-k_Σ_*,

where A is a normalization factor, R is the distance from cluster
center, and k_Σ_⋆ is the power law exponent.
k_Σ_⋆ is found to be ∼2.3 to 2.6 for these primary
clusters.

We denote a half-mass radius, R_1/2, as the radius within which
half of the total mass of the cluster out to R=0.4  pc is
contained. Within our chosen clusters, R_1/2≃0.08-0.1 pc.
Within R_1/2, the total enclosed stellar masses are 4.0×
10^2 M_⊙, 2.3× 10^4 M_⊙, 2.3×
10^4 M_⊙, and 2.2× 10^4 M_⊙ for models
, , , and
, respectively. This shows that the properties of
the primary cluster in the colliding simulations are not much affected
by the choice of star formation subgrid model. Within R_1/2, the
averaged stellar mass surface density is
Σ_*,1/2=2.1× 10^4  M_⊙ pc^-2
(4.4  g cm^-2) for the non-colliding case and 1.1×
10^6  M_⊙ pc^-2 (2.28×10^2 
  g cm^-2), 8.0× 10^5  M_⊙ pc^-2
(1.7×10^2  g cm^-2), and 1.06× 10^6 
  M_⊙ pc^-2 (2.22×10^2  g cm^-2) for the
respective colliding GMC models.  The respective gas masses are
53 M_⊙, 5.7×10^2 M_⊙,
7.7×10^2 M_⊙, and 4.5×10^2 M_⊙.  Note that
the ∼0.1  pc scales are barely resolved in our simulation,
so only the average enclosed gas masses are measured.  At this stage,
the cluster in the non-colliding case has Σ_⋆ <
Σ_ gas. The colliding cases all have
Σ_⋆ > Σ_ gas.

We compare these cluster properties with those of observed young
clusters <cit.>.  The cluster
formed in the non-colliding simulation  that has
Σ_⋆≃ 2×10^4 M_⊙ 
  pc^-2 (4.4  g cm^-2) is much denser than any known
young cluster of comparable mass (i.e., with
M_*,1/2≲1000 M_⊙).  The colliding simulations
produce more massive clusters and these are also seen to have much
higher mass surface densities (by more than a factor of ten) at their
half-mass scale than the densest known Galactic clusters, such as the
Arches or Westerlund 1. We note that stellar feedback is not currently
included within our simulations and expect the implementation of
protostellar outflow feedback, planned in a future paper, to result in
a reduction of Σ_⋆.

The 1-D velocity dispersion, σ_z,1/2, is calculated for the
stars seen to be within projected radii of R_1/2.  The cluster
formed in the non-colliding case has σ_z,1/2=1.57 
  km s^-1, while those formed from GMC collisions have much
higher values of 9.94, 9.83, and 8.81  km s^-1,
respectively for these clusters shown in Figure <ref> (left to
right).  More detailed kinematic analysis is performed in
Section <ref>.

The dynamical state of the clusters is investigated via calculation of
the virial ratio,

    Q_i=-T_⋆/Ω    =    3σ^2R/2GM_⋆,

where T_⋆ and Ω are the total kinetic and gravitational 
potential energies of the stars, respectively. For a given radius R, 
M_⋆ is the total enclosed stellar mass and σ is the 1-D 
velocity dispersion of the enclosed particles.
Values of Q_i<1 indicate a bound cluster, while Q_i=0.5 represents
a state of virial equilibrium. 

We find virial ratios at R_1/2 of 0.34, 0.24, 0.29, and 0.19 for
the primary clusters from the simulations ,
, , and ,
respectively. These are all sub-virial, with collisions forming more
tightly bound clusters at this stage in their evolution. However, as
noted in <ref>, due to R_1/2 approaching the grid
scale, gravitational forces are not well resolved and thus accurate
dynamical evolution of the cluster is limited.  When the virial ratio
is calculated at better-resolved scales of R=0.4  pc, we find
increased values of 0.44, 0.38, 0.36, and 0.33 for the same clusters,
respectively.




  §.§.§ Angular Dispersion Parameter



The angular dispersion parameter (ADP), δ_
  ADP,N(R) <cit.>, is a technique for quantifing the
degree of substructure of a stellar distribution, especially designed
for application to centrally-concentrated star clusters. It is similar
to the azimuthal asymmetry parameter (AAP) developed by
<cit.>. In its simplest form, this technique
divides the distribution spatially into equal-area circular sectors
and compares the dispersion of the number counts contained within each
region. Further division using concentric annuli allows study of this
substructure as a function of radius. In order to account for a global
elongation or eccentricity of the cluster, best-fitted elliptical
annuli can be used. We will adopt this as our fiducial method. To
obtain the best-fit ellipse shape and orientation, a linear fit to the
stars projected in the central 0.2 pc of the cluster is used to set
the position angle, θ_e of the semi-major axis, a. Then the
dispersion in position in directions parallel and perpendicular to
θ_e are calculated to derive the eccentricity. For results in
a given annulus of semi-major axis, a, we display them at a radius,
R, for which the circular area would be equal to that of the
ellipse.

For a given annulus divided into a total of N equal sectors, each
ith sector contains n_i stars. The ADP is defined as:

    δ_ ADP,N=√(1/(N-1)n∑_i=1^N(n_i-n)^2) = √(σ^2/σ^2_ Poisson),

where σ is the standard deviation of the n_i values,
n is the average of the number of stars per sector in the
given annulus, and σ_ Poisson is the standard deviation
expected from a Poisson distribution. Thus, values of δ_
  ADP,N≃ 1 indicate nearly random distributions of sources,
azimuthally.


ADP analysis was performed using 20 equally-spaced concentric annuli
out to a maximum radius of 0.4 pc with N=6 equally-divided sectors.
δ_ ADP,6 is computed for twenty orientations of the sector
pattern at every 3^∘ angular rotation, and the final value is
averaged.  Both circular and elliptical annuli are used to calculate
δ_ ADP,6, using the same previously determined cluster
center.

The star cluster formed in the non-colliding model, which has a
relatively low number of stars and thus larger Poisson errors, has
δ_ ADP,6≃ 1.5 - 2.5 for circular annuli and ≃
1 - 2 for elliptical annuli.  In both cases, δ_ ADP,6
peaks near R=0.15  pc.

The primary clusters from the colliding models have similar
morphologies, especially the density-regulated and
magnetically-regulated m_⋆,min=10 M_⊙ cases.

In the  model, there exists a denser population of
lower mass stars and the location of the subcluster R=0.4  pc
is at a slightly different position. This slight deviation may be
attributed to the earlier onset of star formation from lower density
gas.

The circular and elliptical δ_ ADP,6 values are similar,
although again the latter is slightly smaller in size.  For the
clusters in the  and  simulations,
δ_ ADP,6≃2-3, while the  model has
overall higher values of δ_ ADP,6≃6-7.  The radial
behavior of the clusters from the three colliding models is similar as
well, in that the cluster has moderate values of δ_ ADP,6
out to R≈ 0.2  pc. These decrease at the outskirts of the
defined cluster, then increase relatively sharply out to R≈
0.4  pc upon the presence of the subcluster.

<cit.> carried out a similar ADP analysis of the Orion
Nebula Cluster (ONC). They found δ_ ADP,N with N=4,6,9
rises from below 1 in the very center of the ONC to reach fairly
constant values of about 1.5 to 2.5 from 0.1 pc to about
1.3 pc. Accounting for ellipticity
in the annuli brings these values of δ_ ADP,N down to
about 1 to 1.5, with some variations near 2. Da Rio et al. concluded
the projected spatial distribution of young, embedded stars in the ONC
is relatively smooth, which may be evidence for dynamical processing
if the cluster is older than a few orbit crossing times.

The primary cluster formed in the  model returns
similar values of δ_ ADP,6 as the ONC, while the m_⋆,min=10 M_⊙ colliding models return slightly higher
values.  The  case has much higher δ_
  ADP,6, which may be attributed to its much larger number of stars
(∼10 times higher) resulting in smaller Poisson errors.  As
discussed in <ref>, we caution that the simulation code is
not able to resolve small scale gravitational interactions between
stars, potentially affecting the outcome of this metric. Still, we
expect future application of this technique to better resolved
clusters, including models where final stellar concentrations are
reduced by including local feedback, will provide useful comparisons
with observed young, embedded clusters to help test different
formation scenarios.





  §.§.§ Minimum spanning tree


Another method of studying the hierarchical structure of stellar
distributions is through the use of the minimal spanning tree
(MST). The MST <cit.> is a technique borrowed from graph theory in
which all of the vertices of a connected, undirected graph are joined
such that the total weighting for the graph edges is minimized. In the
case of star clusters, the projected euclidean distances between the
individual stars acts as the edge weight.

To study the hierarchical structure of a collection of stars,
<cit.> introduced a dimensionless
parameter, Q, which can distinguish and quantify between smooth radial
clustering (i.e., more centrally concentrated) vs. multi-scale type
clustering (i.e., more substructure). Specifically,

    Q=s/m.


The numerator is the normalized correlation length

    s=d/R_ cluster,

where d is the mean pairwise separation distance between
the stars and R_ cluster is the overall cluster radius,
calculated as the distance from the mean position of all stars to the
farthest star.

The denominator is the normalized mean edge length

    m=∑_i=1^N_⋆-1e_i/√((N_⋆A)(N_⋆-1)) ,

where N_⋆-1 is the total number of edges, e_i is the
length of each edge, and A=π R_ cluster^2 is the cluster
area.

lcccccr
Parameters of Observed Clusters
0pt

Cluster           Q            s           m 


Taurus                   0.47          0.55         0.26  

IC2391                   0.66          0.74         0.49  

Chameleon                0.67          0.63         0.42  

ρ Ophiuchus         0.85          0.53         0.45  

IC348                    0.98          0.49         0.48  

2<cit.>



The threshold of Q_0=0.785 determines a quantitative threshold of
either smooth radial clustering (Q>Q_0) or multi-scale clustering
(Q<Q_0). Table <ref> lists Q, s, and
m for various observed clusters <cit.>.




We track the evolution of Q throughout our simulations (see
Figure <ref>). In the m_⋆,min=10 M_⊙
non-colliding cases, clustering begins with Q≈ 0.5, at the low
end of observed clusters (e.g., Taurus), and continues to decrease
monotonically. For , an initially much higher Q
is seen, but this quickly decreases to values even lower than the
other cases. This general behavior can be understood as the formation
of an initial cluster (more tightly concentrated in the m_⋆,min=1 M_⊙ case). However, the overall stellar
distribution soon appears very dispersed as other independent clusters
form throughout the GMCs.

The colliding cases exhibit very different behavior. The
density-regulated cases begin with similar Q parameter values and
initially decrease from 0.5 to 0.25. However, beginning near 3.4 Myr,
they experience a sharp increase in Q, reaching between 0.8 and 1.2,
corresponding with the high end of observed clusters (e.g., ρ
Ophiuchus and IC348). The higher-density cutoff models reach higher
maximum Q values and peak at later times. After this peak, the Q
values drop to fairly multi-scale-type clustering, but then rise again
toward 0.5 to 0.7 in the final 0.5 Myr. The m_⋆,min=10 M_⊙ magnetically-regulated cases have similar
qualitative behavior, with the Q peak occurring earlier in time,
near 3.6 Myr and reaching very large values, surpassing centrally
clustered observations. However, Q also drops down to ∼ 0.25
but again equalizes to values near 0.6. These can be understood as the
initial formation of a moderately distributed star cluster which
quickly becomes very centrally dominated as a result of new star
formation in the compressed gas due to the collision that forms a
primary cluster. However, the colliding region soon produces other
clusters separated from the primary star cluster, thus decreasing
Q. Then, beginning from t≈ 4.0 Myr, the gas and stars
continue to coalesce, growing in size, number of stars, and central
concentration. The primary cluster grows and accumulates more of the
surrounding clusters, leading to the growth into a slightly multi-scale
distribution overall. The  model begins near
Q=0.6 and experiences a lower peak near 0.8 before dropping off to
0.25. The final rise of Q is concurrent with the other colliding
models, but instead of settling near 0.6, Q continues to rise until
the end of the simulation t=5.0 Myr, reaching a very high central
clustering value of 1.5.

When comparing the results from our non-colliding vs. colliding
simulations, clusters formed by GMC collisions spend a much greater
fraction of the initial 5 Myr evolution with Q parameters within the
range of observed clusters. While this result should not be
overinterpreted, as clusters produced in the non-colliding cases may
evolve into more centrally peaked distributions at beyond 5 Myr, a
much stronger clustering of stars naturally arises from colliding gas,
and this behavior is quantitatively realized in our simulations.





 §.§ Gas and Star Kinematics






The relationship between the kinematics of young stars and their
surrounding gas has been studied in order to gain insight into the
formation and early evolution of young, embedded stellar
populations. For example, using data from the INfrared Spectra of
Young Nebulous Clusters (IN-SYNC) survey
<cit.>, which achieves radial velocity accuracies of
about 0.3 km s^-1, <cit.> have studied the
kinematic properties of young stars in NGC 1333, while
<cit.> have carried out a similar analysis of
IC 348. <cit.> <cit.> analyzed similar data for
the ONC and its extended southern filament, including comparison to
gas tracers such as ^13CO.  Our simulations allow a similar
investigation of the kinematic properties of both ^13CO-defined
gas and the young stars under various star formation scenarios.

Our gas structures are defined using synthetic ^13CO(J=1-0)
emission, based on the same observational assumptions as Paper II
(i.e., GMCs are at a distance d=3 kpc, the optically thin limit
applies, and we bin with a spectral resolution of 0.212 
  km s^-1). Figure <ref> shows position-velocity
diagrams for non-colliding and colliding cases for density and
magnetically-regulated star formation models 
<cit.>.
The mean gas velocity, mean stellar velocity, and difference in these
means as functions of position are also shown as profiles plotted
below their respective colormaps.  Mean values are taken using
positional bins of 0.5 pc (i.e., 9.5× 10^-3  deg).

The non-colliding cases show widely dispersed gas over the positional
space, with clumpy morphology in ^13CO(J=1-0). The velocity
v_z is fairly low, staying within ± 5  km s^-1. The
gas velocity gradient is relatively shallow, following the general
structure of the clouds. The gas and stellar kinematics in the
density-regulated and magnetically-regulated star formation models are
similar, with more star clusters present in the 
model. The star clusters can be seen localized in positional space
with a small scatter in velocity space. Generally, the stars are
positioned in the vicinity of other high-intensity ^13CO clumps.

The colliding cases show very different behavior in position-velocity
space. The gas is much more localized spatially in x, the direction
of the colliding flows. The structure is more concentrated in the
y-direction as well due to the higher central gravitational
potential formed. The average ^13CO-weighted velocity gradient is
much steeper in x for the colliding cases, but relatively similar in
magnitude to the non-colliding cases in the y line of sight.  Larger
clumps with higher intensities of ^13CO gas are seen in the
colliding cases, with a dense network of filamentary structures
present. Additionally, the gas velocity dispersion is much greater,
with portions reaching velocity dispersions of ± 20 
  km s^-1. Gas and stellar kinematic morphologies are also
similar among the different star formation models.  The central star
cluster is seen to have a very large velocity dispersion, with the
primary clusters in  and  models found
in  <ref> to have σ_z,1/2=9.94 and 9.83 
  km s^-1, respectively.  Separate, smaller clusters can also be
seen with their own stellar populations near high-intensity clumps of
gas.





From the position-velocity information, we calculate velocity
gradients of the gas (dv_ los/ds), the velocity
dispersion of the gas (σ_ gas), the velocity dispersion of
the stars (σ_ *), the ^13CO-weighted average
velocities of the gas (v_ gas), the mass-weighted
average velocities of the stars (v_*), and the velocity
offset between the two (Δv). Table <ref> summarizes these
properties for the four models as viewed from x, y, and z lines
of sight, as well as their RMS values.

We note that the average σ_* values (9.86 and 9.53 km s^-1) in the whole domain are only slightly lower than those
of the primary cluster (9.94 and 9.83 km s^-1). In these
cases, the central cluster contains the majority of the stars and thus
dominates the overall distribution.

Within the non-colliding cases, both the gas and stellar kinematics
agree fairly closely between different star formation models. The
velocity gradients are larger in the x and y directions due to
asymmetries from the impact parameter. A RMS value of
0.072 km s^-1 pc^-1 is recorded. The dispersion of the
gas and stars are similar in both models, with σ_ gas,
  RMS≈ 1.76 km s^-1 and σ_ *, RMS≈
1.68 km s^-1. The mean velocities of the gas and stars are
also similar, with v_ gas, RMS≈
0.91 km s^-1 and v_ *, RMS≈
0.94 km s^-1. The stellar velocities in the density-regulated
model have slightly lower dispersions, but higher mean
values. Overall, the velocity offset between the gas and the stars for
the non-colliding models is approximately 0.35 km s^-1.

For the colliding cases, the density-regulated and
magnetically-regulated star formation models exhibit similar kinematic
properties of the gas and stars. On average, dv_
    los/ds=0.26 km s^-1 pc^-1, σ_ gas,
  RMS=4.14 km s^-1, σ_ *,
  RMS=9.70 km s^-1, v_ gas,
  RMS=0.50 km s^-1, v_ *,
  RMS=1.16 km s^-1, and Δv=1.35 km s^-1. Relative to non-colliding clouds,
the collision induces a much larger velocity gradient (∼3-4 times
greater), a larger velocity dispersion in the gas (∼2 times
greater), and a much larger stellar velocity dispersion (∼5 times
greater).

As functions of position, the mean velocities of gas, stars, and their
offsets, are compared. Offsets exist in both non-colliding and
colliding cases, becoming most apparent in close proximity to star
clusters.  For the non-colliding case, these differences are
relatively small, at a few km s^-1. However, the colliding
case contains regions in which the offsets exceed 5 
  km s^-1.  Averaged over position space, colliding GMCs result
in velocity offsets a factor of ∼4 times higher than
non-colliding GMCs. This parameter may be an indicator for determining
the dynamical formation history of young star clusters, as results
from collisions show more disturbance kinematically.

We make a first, simple comparison with the results of the IN-SYNC
survey of the ONC and surroundings <cit.>. This
survey found offsets between gas and star velocities of approximately
Δ v_r∼ -0.5 km/s, but up to -1.0 to -1.5 km/s in some
regions. The magnitude of such offsets are in general consistent with
those seen in both the non-colliding and colliding models, especially
considering variations associated with the particular line of sight.
Given that the observational data for Orion is just a single example
of a star-forming region, viewed on a particular sight line, it is
difficult to draw definitive conclusions about whether the gas and
star kinematics favor one scenario over another. Larger numbers of
star-forming regions need to be studied with similar methods. In
addition, other metrics, such as the comparison of low and high
density gas tracers <cit.>,
need to be examined, which on the simulation side requires extension 
of the astrochemical modeling to include species such as N_2H^+.





§ DISCUSSION AND CONCLUSIONS



We have implemented two classes of star formation sub-grid routines
into the MHD code Enzo that we are using to study GMC collisions: a
density-regulated model based on a threshold density and a new
magnetically-regulated model based on a threshold mass-to-flux
ratio. Varying key parameters for each star formation routine, we
explored the large-scale morphology, properties of star-forming gas,
global star formation rates and efficiencies over time, spatial
clustering of the stars, and gas and stellar kinematics. For each
model, we investigated scenarios of non-colliding and colliding GMCs.

The non-colliding cases evolved in a relatively quiescent manner,
driven by the initial turbulence and interplay of self-gravity and
magnetic fields. Star clusters formed only in the very late stages of
the simulations, from overdense clumps located within filaments and
dispersed throughout the GMC complex. Generally, these clusters
contained hundreds of solar masses each and grew at a relatively slow
rate. Star clusters in the density-regulated star formation routines
were smaller and more isolated. The clusters formed in the
magnetically-regulated models exhibited slightly more elongated
morphologies. For this simulation set up, the level of star formation
activity appears to be regulated by the effective density threshold
that is used in each of these models, with the mass-to-flux criterion
not having a large influence.

During collisions between GMCs, stars formed earlier and in larger
clusters, from high-density gas produced in the primary filamentary
colliding region. 
While star formation rates level off by the completion of the 
simulations, extrapolation of future behavior is unclear. 
Nevertheless, by t=5  Myr, individual clusters have grown
and merged to form one large, dominant cluster with a total stellar
mass of 5× 10^4 M_⊙. For this particular set-up, the
final overall level of star formation is relatively independent of all
the explored star formation sub-grid models. Star formation appears to
be limited by the ability of the collision to direct mass into high
density regions, which then eventually form stars with high overall
efficiency.

Just prior to star formation, both density and magnetically regulated
star formation result in fairly similar gas properties of parent
cells. However, colliding cases experience relatively wider ranges of
densities, temperatures, μ_ cell, and velocity
magnitude. Higher mean values for density and temperature are found,
while gas is more magnetically subcritical and turbulent.

The primary star clusters formed in the various models were analyzed
and found to have much higher surface densities at their half-mass
scale than any observed cluster. We expect the future inclusion of
stellar feedback will reduce these surface densities. 
The angular dispersion parameter (ADP) analysis was carried out on the
primary clusters in the simulations. ADP values are generally greater
than those see in the ONC, which may indicate the ONC is dynamically
older than the simulated clusters.
The MST Q parameter was also used to investigate the global spatial
distribution properties of the star, with non-colliding cases
resulting in overall highly multi-scaled clustering due to the
scattered formation of independent clusters. Colliding GMCs produce
clusters with Q parameters that vary between those expected of
multi-scale and centrally clustered distributions.

Kinematically, our colliding GMC cases produce velocity gradients 3-4
times greater than those of the non-colliding cases. The velocity
dispersions also differ, with the gas in the colliding clouds having
approximately twice the velocity dispersion. Stellar velocity
dispersions in the simulations are dominated by the potentials of the
primary clusters that form, with this leading to much greater
dispersions in the colliding case. We find that the colliding cases
produce typically 4 times larger offsets between the mean gas and mean
star velocities compared to the non-colliding case.

Finally, we remind the reader of several important caveats. The young
stars do not inject feedback, especially protostellar outflow
feedback, into the surrounding gas. Internal star cluster dynamics are
not well followed because of gravitational softening at the grid scale
(∼ 0.1 pc) and because the star particles lack realistic mass and
multiplicity distributions. Still the conditions that are simulated
here may provide boundary conditions for more detailed models that are
able to follow full N-body evolution of the clusters
<cit.>. Finally, in the context of GMC
collisions, a wide range of cloud (e.g., degree of initial
magnetization) and collision (e.g., velocities; impact parameters)
parameters remain to be explored with these models. These items will
be addressed in subsequent papers in this series.



The authors would like to thank Nicola Da Rio and
Juan Pablo Farias for useful discussions.  Computations described in
this work were performed using the publicly-available 
code (http://enzo-project.org), This research also made use of the
yt-project (http://yt-project.org), a toolkit for analyzing and
visualizing quantitative data <cit.>. Volume renderings
were performed using 
(https://wci.llnl.gov/simulation/computer-codes/visit). The authors
acknowledge University of Florida Research Computing (www.rc.ufl.edu)
for providing computational resources and support that have
contributed to the research results reported in this publication.

Enzo (Bryan et al. 2014), Grackle (Bryan et al. 2014; Kim et al. 2014), PyPDR, yt (Turk et al. 2011), VisIt





[Anathpindika(2009)]Anathpindika_2009
Anathpindika, S. 2009, , 504, 437

[Balfour et al.(2015)Balfour, Whitworth, Hubber, &
  Jaffa]Balfour_ea_2015
Balfour, S. K., Whitworth, A. P., Hubber, D. A., & Jaffa, S. E. 2015,
  , 453, 2471

[Barrow et al.(1985)Barrow, Bhavsar, &
  Sonoda]Barrow_ea_1985
Barrow, J. D., Bhavsar, S. P., & Sonoda, D. H. 1985, , 216, 17

[Bryan et al.(2014)Bryan, Norman, O'Shea, Abel, Wise,
  Turk, Reynolds, Collins, Wang, Skillman, Smith, Harkness,
  Bordner, Kim, Kuhlen, Xu, Goldbaum, Hummels, Kritsuk, Tasker,
  Skory, Simpson, Hahn, Oishi, So, Zhao, Cen, Li, & The Enzo
  Collaboration]Bryan_ea_2014
Bryan, G. L., Norman, M. L., O'Shea, B. W., et al. 2014, , 211, 19

[Butler & Tan(2009)]Butler_Tan_2009
Butler, M. J., & Tan, J. C. 2009, , 696, 484

[Butler & Tan(2012)]Butler_Tan_2012
—. 2012, , 754, 5

[Butler et al.(2015)Butler, Tan, & Van
  Loo]Butler_ea_2015
Butler, M. J., Tan, J. C., & Van Loo, S. 2015, , 805, 1

[Cartwright & Whitworth(2004)]Cartwright_Whitworth_2004
Cartwright, A., & Whitworth, A. P. 2004, , 348, 589

[Chira et al.(2013)Chira, Beuther, Linz, Schuller,
  Walmsley, Menten, & Bronfman]Chira_ea_2013
Chira, R.-A., Beuther, H., Linz, H., et al. 2013, , 552, A40

[Cottaar et al.(2014)Cottaar, Covey, Meyer, Nidever,
  Stassun, Foster, Tan, Chojnowski, da Rio, Flaherty, Frinchaboy,
  Skrutskie, Majewski, Wilson, & Zasowski]Cottaar_ea_2014
Cottaar, M., Covey, K. R., Meyer, M. R., et al. 2014, , 794, 125

[Cottaar et al.(2015)Cottaar, Covey, Foster, Meyer, Tan,
  Nidever, Chojnowski, da Rio, Flaherty, Frinchaboy, Majewski,
  Skrutskie, Wilson, & Zasowski]Cottaar_ea_2015
Cottaar, M., Covey, K. R., Foster, J. B., et al. 2015, , 807, 27

[Crutcher(2012)]Crutcher_2012
Crutcher, R. M. 2012, , 50, 29

[Da Rio et al.(2014)Da Rio, Tan, & Jaehnig]DaRio_ea_2014
Da Rio, N., Tan, J. C., & Jaehnig, K. 2014, , 795, 55

[Da Rio et al.(2017)Da Rio, Tan, Covey, Cottaar, Foster,
  Cullen, Tobin, Kim, Meyer, Nidever, Stassun, Chojnowski,
  Flaherty, Majewski, Skrutskie, Zasowski, &
  Pan]DaRio_ea_2017arXiv
Da Rio, N., Tan, J. C., Covey, K. R., et al. 2017, ArXiv e-prints,
  arXiv:1702.04113

[Dedner et al.(2002)Dedner, Kemm, Kröner, Munz,
  Schnitzer, & Wesenberg]Dedner_ea_2002
Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of
  Computational Physics, 175, 645

[Dobbs et al.(2015)Dobbs, Pringle, &
  Duarte-Cabral]Dobbs_ea_2015
Dobbs, C. L., Pringle, J. E., & Duarte-Cabral, A. 2015, , 446,
  3608

[Duarte-Cabral et al.(2011)Duarte-Cabral, Dobbs, Peretto, &
  Fuller]Duarte-Cabral_ea_2011
Duarte-Cabral, A., Dobbs, C. L., Peretto, N., & Fuller, G. A. 2011,
  , 528, A50

[Ester et al.(1996)Ester, Kriegel, Sander, &
  Xu]Ester_ea_1996
Ester, M., Kriegel, H., Sander, J., & Xu, X. 1996, in  (AAAI Press),
  226–231

[Farias et al.(2017)Farias, Tan, &
  Chatterjee]Farias_2017arXiv
Farias, J. P., Tan, J. C., & Chatterjee, S. 2017, ArXiv e-prints,
  arXiv:1701.00701

[Few et al.(2016)Few, Dobbs, Pettitt, &
  Konstandin]Few_ea_2016
Few, C. G., Dobbs, C., Pettitt, A., & Konstandin, L. 2016, ,
  460, 4382

[Foster et al.(2015)Foster, Cottaar, Covey, Arce, Meyer,
  Nidever, Stassun, Tan, Chojnowski, da Rio, Flaherty, Rebull,
  Frinchaboy, Majewski, Skrutskie, Wilson, &
  Zasowski]Foster_ea_2015
Foster, J. B., Cottaar, M., Covey, K. R., et al. 2015, , 799, 136

[Fujimoto et al.(2014)Fujimoto, Tasker, &
  Habe]Fujimoto_ea_2014
Fujimoto, Y., Tasker, E. J., & Habe, A. 2014, , 445, L65

[Ginsburg et al.(2012)Ginsburg, Bressert, Bally, &
  Battersby]Ginsburg_ea_2012
Ginsburg, A., Bressert, E., Bally, J., & Battersby, C. 2012, ,
  758, L29

[Gutermuth et al.(2005)Gutermuth, Megeath, Pipher,
  Williams, Allen, Myers, & Raines]Gutermuth_ea_2005
Gutermuth, R. A., Megeath, S. T., Pipher, J. L., et al. 2005, ,
  632, 397

[Habe & Ohta(1992)]Habe_Ohta_1992
Habe, A., & Ohta, K. 1992, , 44, 203

[Hacar et al.(2016)Hacar, Alves, Forbrich, Meingast,
  Kubiak, & Großschedl]Hacar_ea_2016
Hacar, A., Alves, J., Forbrich, J., et al. 2016, , 589, A80

[Haworth et al.(2015a)Haworth, Shima, Tasker,
  Fukui, Torii, Dale, Takahira, & Habe]Haworth_ea_2015a
Haworth, T. J., Shima, K., Tasker, E. J., et al. 2015a,
  , 454, 1634

[Haworth et al.(2015b)Haworth, Tasker, Fukui,
  Torii, Dale, Shima, Takahira, Habe, &
  Hasegawa]Haworth_ea_2015b
Haworth, T. J., Tasker, E. J., Fukui, Y., et al. 2015b,
  , 450, 10

[Heitsch et al.(2009)Heitsch, Stone, &
  Hartmann]Heitsch_Stone_Hartmann_2009
Heitsch, F., Stone, J. M., & Hartmann, L. W. 2009, , 695, 248

[Henshaw et al.(2014)Henshaw, Caselli, Fontani,
  Jiménez-Serra, & Tan]Henshaw_ea_2014
Henshaw, J. D., Caselli, P., Fontani, F., Jiménez-Serra, I., &
  Tan, J. C. 2014, , 440, 2860

[Henshaw et al.(2013)Henshaw, Caselli, Fontani,
  Jiménez-Serra, Tan, & Hernandez]Henshaw_ea_2013
Henshaw, J. D., Caselli, P., Fontani, F., et al. 2013, , 428,
  3425

[Inutsuka et al.(2015)Inutsuka, Inoue, Iwasaki, &
  Hosokawa]Inutsuka_ea_2015
Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, ,
  580, A49

[Jaehnig et al.(2015)Jaehnig, Da Rio, &
  Tan]Jaehnig_ea_2015
Jaehnig, K. O., Da Rio, N., & Tan, J. C. 2015, , 798, 126

[Klein & Woods(1998)]Klein_Woods_1998
Klein, R. I., & Woods, D. T. 1998, , 497, 777

[Kong et al.(2017)Kong, Tan, Caselli, Fontani, Wang, &
  Butler]Kong_ea_2017arXiv
Kong, S., Tan, J. C., Caselli, P., et al. 2017, ArXiv e-prints,
  arXiv:1701.05953

[Krumholz & McKee(2005)]Krumholz_McKee_2005
Krumholz, M. R., & McKee, C. F. 2005, , 630, 250

[Krumholz & Tan(2007)]Krumholz_Tan_2007
Krumholz, M. R., & Tan, J. C. 2007, , 654, 304

[Lim & Tan(2014)]Lim_Tan_2014
Lim, W., & Tan, J. C. 2014, , 780, L29

[Ma et al.(2013)Ma, Tan, & Barnes]Ma_ea_2013
Ma, B., Tan, J. C., & Barnes, P. J. 2013, , 779, 79

[Mouschovias & Spitzer(1976)]Mouschovias_Spitzer_1976
Mouschovias, T. C., & Spitzer, Jr., L. 1976, , 210, 326

[Nakano & Nakamura(1978)]Nakano_Nakamura_1978
Nakano, T., & Nakamura, T. 1978, , 30, 671

[Parravano et al.(2011)Parravano, McKee, &
  Hollenbach]Parravano_ea_2011
Parravano, A., McKee, C. F., & Hollenbach, D. J. 2011, , 726, 27

[Pillai et al.(2006)Pillai, Wyrowski, Carey, &
  Menten]Pillai_ea_2006
Pillai, T., Wyrowski, F., Carey, S. J., & Menten, K. M. 2006, ,
  450, 569

[Rathborne et al.(2006)Rathborne, Jackson, &
  Simon]Rathborne_ea_2006
Rathborne, J. M., Jackson, J. M., & Simon, R. 2006, , 641, 389

[Sakai et al.(2008)Sakai, Sakai, Kamegai, Hirota,
  Yamaguchi, Shiba, & Yamamoto]Sakai_ea_2008
Sakai, T., Sakai, N., Kamegai, K., et al. 2008, , 678, 1049

[Scoville et al.(1986)Scoville, Sanders, &
  Clemens]Scoville_Sanders_Clemens_1986
Scoville, N. Z., Sanders, D. B., & Clemens, D. P. 1986, , 310, L77

[Smith et al.(2016)Smith, Bryan, Glover, Goldbaum, Turk,
  Regan, Wise, Schive, Abel, Emerick, O'Shea, Anninos, Hummels,
  & Khochfar]Smith_ea_2016arXiv
Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2016, ArXiv
  e-prints, arXiv:1610.09591

[Stutz & Gould(2016)]Stutz_Gould_2016
Stutz, A. M., & Gould, A. 2016, , 590, A2

[Takahira et al.(2014)Takahira, Tasker, &
  Habe]Takahira_ea_2014
Takahira, K., Tasker, E. J., & Habe, A. 2014, , 792, 63

[Tan(2000)]Tan_2000
Tan, J. C. 2000, , 536, 173

[Tan et al.(2014)Tan, Beltrán, Caselli, Fontani,
  Fuente, Krumholz, McKee, & Stolte]Tan_ea_2014
Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, Protostars
  and Planets VI, 149

[Tasker & Tan(2009)]Tasker_Tan_2009
Tasker, E. J., & Tan, J. C. 2009, , 700, 358

[Truelove et al.(1997)Truelove, Klein, McKee, Holliman,
  Howell, & Greenough]Truelove_ea_1997
Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, , 489,
  L179

[Turk et al.(2011)Turk, Smith, Oishi, Skory, Skillman,
  Abel, & Norman]Turk_ea_2011
Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, , 192, 9

[Van Loo et al.(2013)Van Loo, Butler, &
  Tan]Van_Loo_Butler_Tan_2013
Van Loo, S., Butler, M. J., & Tan, J. C. 2013, , 764, 36

[Van Loo et al.(2015)Van Loo, Tan, &
  Falle]VanLoo_ea_2015
Van Loo, S., Tan, J. C., & Falle, S. A. E. G. 2015, , 800, L11

[Wang & Abel(2008)]Wang_Abel_2008
Wang, P., & Abel, T. 2008, , 672, 752

[Wang et al.(2008)Wang, Zhang, Pillai, Wyrowski, &
  Wu]Wang_ea_2008
Wang, Y., Zhang, Q., Pillai, T., Wyrowski, F., & Wu, Y. 2008, ,
  672, L33

[Wu et al.(2017)Wu, Tan, Nakamura, Van Loo, Christie, &
  Collins]Wu_ea_2017
Wu, B., Tan, J. C., Nakamura, F., et al. 2017, , 835, 137

[Wu et al.(2015)Wu, Van Loo, Tan, &
  Bruderer]Wu_ea_2015
Wu, B., Van Loo, S., Tan, J. C., & Bruderer, S. 2015, , 811, 56

[Zuckerman & Evans(1974)]Zuckerman_Evans_1974
Zuckerman, B., & Evans, II, N. J. 1974, , 192, L149






