{'fraction_non_alphanumeric': 0.03778344877870944, 'fraction_numerical': 0.016594969012030624, 'mean_word_length': 4.4369301514310635, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular cloud formation in high-shear, magnetized colliding flows

The colliding flows (CF) model is a well-supported mechanism for generating molecular clouds. However, to-date most CF simulations have focused on the formation of clouds in the normal-shock layer between head-on colliding flows. We performed simulations of magnetized colliding flows that instead meet at an oblique-shock layer. Oblique shocks generate shear in the post-shock environment, and this shear creates inhospitable environments for star formation. As the degree of shear increases (i.e. the obliquity of the shock increases), we find that it takes longer for sink particles to form, they form in lower numbers, and they tend to be less massive. With regard to magnetic fields, we find that even a weak field stalls gravitational collapse within forming clouds. Additionally, an initially oblique collision interface tends to reorient over time in the presence of a magnetic field, so that it becomes normal to the oncoming flows. This was demonstrated by our most oblique shock interface, which became fully normal by the end of the simulation.


INTRODUCTION
The interstellar medium (ISM) is a dynamic environment where material cycles from a warm, tenuous phase to a cold, dense phase and back again. The evolution of the gas through these phases establishes limits on star formation in the galaxy. Understanding the details of such dynamics is crucial as theoretical scenarios for molecular cloud evolution (and hence star formation) have shifted away from models based on long molecular cloud lifetimes (i.e. steady states where τ cloud >> τ f f ; cf. Shu, Adams & Lizano (1987); Mouschovias (1991)) . Over the last decade or more, evidence has grown that instead supports a scenario in which clouds are transient structures born out of the cold, dense phase of the ISM. For example, starless clouds appear to be scarce, meaning clouds do not slowly evolve towards conditions in which star formation begins (Beichman et al. 1986;Dame, Hartmann & Thaddeus 2001). Also, the majority of stars in clouds are, in general, young showing ages < 5 M yr (Hartmann, Ballesteros-Paredes & Bergin 2001). This implies that star formation begins soon after the cloud itself forms and that clouds are not generally long-lived due to the absence of older stars (Fukui & Yonekura 1998 Carpenter 2000). Finally, high degrees of hierarchical structure in molecular clouds should have short lifetimes due to star-star scattering or tidal interactions (Lada & Lada 1995;Eisenhauer et al. 1998;Beck, Kelly & Lacy 1998).
While observational support for short cloud lifetimes has grown steadily, we note that scenarios invoking rapid cloud and star formation are not new (Hunter 1979;Larson 1981;Hunter et al. 1986;Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999a). Given the complexity that comes with a dynamical theory of molecular cloud evolution, however, exploring the full dynamics of these scenarios depends heavily on high performance computational methods. One scenario that generates transient molecular clouds and whose exploration has been possible with modern numerical simulations is the 'colliding flows' model of molecular cloud formation. In this model, molecular clouds are formed in the shocked collision layer between two largescale colliding streams of gas. Simulations of colliding flows have shown that nonlinear density structures readily form in the shocked collision region between the flows via a variety of instabilities (cf. Heitsch et al. (2005) for discussion of the unstable modes). Further, these structures develop column densities high enough for effective UV-shielding, and thus, H2 formation (N ≈ 1 − 2 cm −2 ; van Dishoeck & Black (1988); van Dishoeck & Blake (1998)). Gravitational instabilities then cause the turbulent, shocked gas in the dense structures to collapse and form stars. Given the highly dynamical environment, the transition from the beginning of molecular cloud formation to star formation to cloud destruction occurs in roughly a dynamical time (Audit & Hennebelle 2005;Heitsch et al. 2006;Vázquez-Semadeni et al. 2006;Heitsch et al. 2008), matching observations (Elmegreen (2000); Ballesteros-Paredes & Hartmann (2007), and references therein). The clouds produced in colliding flows simulations are similar in many ways to those seen in observation.
While colliding flows models tend to be idealized, they are not without motivation. Coherent large-scale streams of gas are plausible in many situations in the ISM, such as, expanding bubbles driven from energetic OB associations and/or supernovae, turbulent motions arising from gravitational instabilities, density waves in the spiral arms of galaxies, and cloud-cloud collisions Inutsuka et al. 2015). Observational evidence supporting these various scenarios, as well as their association with molecular cloud formation, takes a number of forms. Atomic inflows surrounding molecular gas have been observed in Taurus (Ballesteros-Paredes, Hartmann & Vázquez-Semadeni 1999b) and other molecular clouds (Brunt 2003). Looney et al. (2006) show that the active star forming molecular cloud core in BD +40 4124 likely arose from cloud-cloud collisions (a localized version of a colliding flow). Molecular clouds have been found at the edges of supershells (Dawson et al. 2011(Dawson et al. , 2013, which are a form of colliding flows driven by multiple supernovae. If molecular clouds form via the accumulation of gas on dynamical timescales (Hartmann, Ballesteros-Paredes & Bergin 2001), then magnetic fields are expected to be dynamically important on similar timescales. Some colliding flows simulations have studied the role of magnetic fields in molecular cloud formation (Heitsch et al. 2007;Hennebelle et al. 2008;Banerjee et al. 2009;Heitsch, Stone & Hartmann 2009;Vázquez-Semadeni et al. 2011;Chen & Ostriker 2014;Körtgen & Banerjee 2015). Some of these simulations (Heitsch et al. 2007;Hennebelle et al. 2008;Banerjee et al. 2009) have found power law relations between the magnetic field strength and density, similar to what is seen observationally, i.e. B ∝ n k , where 1/2 < k < 2/3 for n > 100 cm −3 , and k = 0 for n < 100 cm −3 (Troland & Heiles 1986;Crutcher 1999;Crutcher et al. 2010;Tritsis et al. 2015). The inclusion of fields in the models also reduces the star formation rates to values more in agreement with observations (Vázquez-Semadeni et al. 2011;Chen & Ostriker 2014). Such a reduction has been claimed to result from lower degrees of turbulent substructure found in the simulations (Heitsch et al. 2007;Hennebelle et al. 2008;Heitsch, Stone & Hartmann 2009;Chen & Ostriker 2014).
While the colliding flows model has proved useful for understanding transient molecular cloud formation, the great majority of work has focused on head-on collisions, with no obliquity in the initial shocks formed in the flow. It is therefore worthwhile to explore the consequences of allowing the flows to interact at an interface initially inclined relative to incoming velocity vectors. Strong shear in the interaction region can lead to stretching of embedded field lines and the possible generation of turbulence from KH modes. The role of shear in molecular cloud and star formation has begun to be investigated numerically. Hydrodynamic simulations by Rey-Raposo, Dobbs & Duarte-Cabral (2015) show that clouds can inherit shear velocity fields during their formation in spiral arm galaxies, and that these shear flows impair subsequent star formation. Körtgen & Banerjee (2015) find that varying the intersection angle of magnetized colliding flows reduces the star formation efficiency of the gas. This was attributed to a post-shock shear flow disrupting the formation of high density structures.
Continuing along these lines, we address the role of shear and magnetic fields in molecular cloud formation, with a focus on the bulk dynamics of the flow. Shear is generated in our models by keeping the flows parallel and varying the angle of the collision interface. We present four adaptive mesh simulations, at a peak resolution of 0.05 pc, varying the collision interface from a normal incidence to highly inclined. All of our simulations include a uniform magnetic field aligned with the flows to simulate idealized ISM conditions. We also include self gravity and time-dependent cooling. We present our numerical model in Section 2, our results in Sections 4-9, and our discussion in Section 10.

NUMERICAL MODEL
We conducted our simululations using AstroBEAR 1 (Carroll-Nellenback et al. 2013), a publicly available, massively parallelized, adaptive mesh refinement (AMR) code that contains a variety of multiphysics solvers (i.e. selfgravity, magnetic resistivity, radiative transport, ionization dynamics, heat conduction, and more). Our setup was of two, 40 pc diameter cylinders colliding in a 3D domain under the influence of gravity, magnetic fields, and cooling ( Fig.  1). Gravity and cooling source terms were solved using a Strang-split corner upwind transport (CTU) scheme that was 2nd-order accurate in time. This was combined with a directionally unsplit CTU scheme for the 3D ideal magnetohydrodynamics (MHD) equations, using the HLLD Riemann solver. Gravitational interactions included both the self-gravity of the gas, as well as the gravitational acceleration due to sink particles, which were implemented following Federrath et al. (2010a). To solve Poisson's equation for the gravitational potential of the grid, AstroBEAR uses HYPRE 2 (Falgout & Yang (2002); see the appendix in Kaminski et al. (2014) for a description of AstroBEAR's self-gravity algorithm).
A uniform magnetic field was initialized everywhere in the box, parallel to the flows. The field was initially dynamically weak, with β = 10 and βram ≈ 38 at the start of the simulations (β is the ratio of thermal to magnetic pressure, and βram is the ratio of ram to magnetic pressure). The field had an initial strength of B = 1.3 µG, which is at the lower end of current ISM magnetic field estimates (Beck 2001;Heiles & Troland 2005). Cooling and heating were included using a parametrized cooling curve adapted from Inoue & Inutsuka (2008) to include the effects of UV shielding. The Figure 1. Model diagram. Two oppositely driven cylindrical flows intersected at an angle θ, rotated about the y− axis. The interface was perturbed with a random sequence of sine waves. The flows were homogeneous and embedded in a stationary, uniform ambient medium of the same density and pressure. The setup was initialized with a uniform magnetic field in x, parallel to the direction of the colliding flows. modification allowed the gas to cool to T = 10 K for densities greater than n > 1000 cm −3 (Ryan & Heitsch in prep).
Our simulations tracked the formation of gravitationally collapsed objects (sink particles) which were created only after a series of checks had been satisfied, e.g., (i) a cell and surrounding region was Jeans unstable (ii) the center cell (i.e. location of potential sink) in this Jeans unstable region was collocated with a gravitational potential minimum (iii) the surrounding region was exhibiting a converging flow into this center zone etc., as given in Federrath et al. (2010a). These checks ensured that sinks formed only when appropriate, i.e. when a region would go on to form a gravitationally bound object if better resolution were available. The simulations had 5 levels of AMR, giving a finest cell size of ∆xmin = 0.05 pc. Given this resolution, we identified single sink particles as protoclusters rather than protostars. We will therefore use the words 'sink' and 'protocluster' interchangeably.
Sink particles interacted with the gas (and other sinks around them) through gravitational interactions, and the ability to accrete surrounding material. They remained individual objects throughout the course of the simulation (i.e. did not merge), and did not provide any form of energy or momentum feedback into the surrounding medium (i.e no winds or radiation). Gas around a sink particle was accreted only when the density in the surrounding zones exceeded a given threshold, dictated by the Truelove condition (Truelove et al. 1997), where λJ is the cell-centered Jeans length. Then, only the excess gas was removed (Federrath et al. 2010a). The suite of simulations consisted of four runs, each with a different orientation of the collision interface, specified by the inclination angle θ ( Table 1). The inclination angle was varied between θ = 0 • and θ = 60 • , i.e. between a head-on and highly inclined collision. This translated into the difference between a normal shock at the collision layer and an oblique shock. All parameters for the present suite of runs were the same as the 'smooth' model of Carroll-Nellenback, Frank & Heitsch (2014), with the exception of the magnetic field and the variation of the collision interface (the smooth run was a hydro, head-on collision). We will therefore compare the θ = 0 • run to the smooth run to study the effect of the magnetic field alone. Hence, the smooth run from here on out will be called the 'hydro version of the θ = 0 • case' (or, the 'hydro run', for short).
The flows were injected into a stationary ambient medium at a velocity of v = 11 km s −1 and an isothermal mach of M = 1.5. The mean molecular weight was set to µ = 1.27, and the adiabatic exponent to γ = 5/3. The gas was initially in thermal equilibrium at a uniform number density of n = 1 cm −3 , corresponding to the (linearly) stable, warm neutral medium (WNM) phase of the cooling curve. This set the thermal pressure and temperature everywhere inside the box to be P therm k −1 B = 4931 K cm −3 , and T = 4931 K, where kB is the Boltzmann constant. The ram pressure of each flow was Pramk −1 B = 18, 500 K cm −3 , giving a total mass flux into the collision region of M f lux ≈ 886 M M yr −1 .
As in previous simulations, the collision interface was seeded with a set of random sinusoidal perturbations to excite the nonlinear thin shell (NTS) and Kelvin Helmholtz (KH) instabilities (Heitsch et al. 2006;Carroll-Nellenback, Frank & Heitsch 2014). These perturbations had a maximum amplitude of A = 2 pc, spectral index α = −2.0, and maximum wave number kmax = 16 pc −1 . Boundary conditions on the box were set to inflow-only on the faces where the flows were injected, and extrapolating on all other faces. The boundary conditions for the gravity solver were set to multipole expansion.
The simulations were initialized to 3 levels of AMR, with the finest cells centered on the collision interface within a cylindrical volume 40 pc in diameter and 20 pc long, and the coarser meshes nested outwards from there. This made for an initial effective resolution of ∆x ef f = 0.2 pc. Two additional levels were triggered throughout the simulation based on gradients in the fluid variables, as well as resolution of the local Jeans length (λJ ), such that if a cell's λJ was smaller than 64 zones on that cell's level, another level of AMR would be added. This brought the finest cell size to ∆xmin = 0.05 pc, as previously stated. The final simulation time for all of the runs was tsim = 27.3 M yr, with the exception of the θ = 60 • . This was the only run that did not form any sink particles by this time, and so was extended out until tsim = 32.8 M yr.
In what follows, our analysis will focus on the θ = 0, 15, and 60 • cases, as the θ = 30 • case did not significantly differ from the θ = 15 • case.

GENERATING SHEAR VIA AN OBLIQUE COLLISION INTERFACE
The jump conditions across an oblique, 1D, shock-bounded slab convert initially intersecting velocity vectors into a shear flow field. As a measure of this shear for our setup here, we produced mass weighted histograms of the magnitude of the vorticity (||∇ × v||) in cylindrical analysis regions centered on the collision region. Each of these "hockey pucks" encompassed the corresponding collision region by tracing the collision interface and extending out to 5 pc on either side of the interface. They were normalized to contain the same amount of mass (1, 000M ). Histograms were generated at t = 1 M yr, which is approximately the time it would take a post-shock sound wave to travel from the center of the CF cylinder to the outer boundary. As can be seen in Fig. 2, lower shear runs (i.e., θ = 0−15 • ) have significantly more mass at lower vorticity (||∇× v|| < 100), than at higher. As the collision angle increases to θ = 30 • , we see two changes occur. First, the amount of mass at ||∇ × v|| ≤ 1 decreases considerably. Second, more of the mass moves to higher vorticity. As the collision angle steepens to 60 • , this trend strengthens. For the θ = 60 • case, there is neglibible mass at the lowest ||∇ × v||. That is to say, nearly all of the mass has acquired vorticity in this run. Moreover, most of the mass has acquired high vorticity (||∇ × v|| > 10).
Thus, as the inclination angle increases, more vorticity is generated in the collision layer. This vorticity can be associated with the solenoidal mode of the post-shock turbulence (Federrath et al. 2010b). Turbulent solenoidal fluid motions are efficient in providing support against collapse. Compared to compressive modes of turbulence, Federrath & Klessen (2012) show that solenoidal modes greatly reduce the star formation rate in turbulent, magnetized clouds. The passage of gas through the oblique shocks of our simulations thus transforms the compressive nature of the flows into solenoidal. In this way, our higher shear cases should exhibit greater degrees of turbulent support.

MASS-TO-FLUX RATIO
Before we present the results of the simulations, we will first discuss some general issues associated with gravitational stability and magnetic fields that are relevant to our studies. We begin with the critical mass-to-flux ratio (M2FR) of cylindrical uniform flows. The M2FR compares the relative strength of the magnetic field to gravity. It does not take into account other forces which oppose gravity, such as thermal or ram pressure forces. The M2FR is given by, where Σ is the mass column density and B is the magnetic field threading the cylinder (Nakano & Nakamura 1978). Equation 2 can be rearranged for the critical length of the cylinder. In terms of typical ISM values this is given by, where n is the number density (Vázquez-Semadeni et al. 2011). For the initial WNM values of our flows (n = 1 cm −3 , B = 1.3 µG), Eqn. 3 gives a critical length of, for the WNM component of the gas. Using this, we can ask whether we would expect the warm gas to be magnetically supercritical. We first have to define a length scale over which we will check stability. Naturally this would be the collision region, as this is where the molecular clouds will be forming. As we will see, the width of this region (L coll ) is on the order of 10 pc. This means that over the length scale of the collision region, the flow would be sub-critical. That is, given Lcrit >> L coll , the magnetic field, at least initially, should be strong enough to withstand gravity. Given the mass flux into the collision region along the colliding flows, one can ask 'how long will it take to accumulate enough mass into the collision region for the collision region to go unstable?' For this, it is helpful to recast Eqn.
(3) in terms of the critical column number density, -Semadeni et al. 2011). For simplicity, we imagine the collision region to be initially massless. Its column density as a function of time is then, where v is the speed of the flows. Equating this to Eqn.
(5) and solving for t gives the timescale for the collision region to become magnetically supercritical: tcrit ≈ 6 M yr. Thus, the collision region is expected to very quickly (i.e. tcrit << tsim) become magnetically supercritical. While this is a crude estimate, when taken in an average sense (over the collision region) it implies that the mean field is unable to support the collision region against collapse. However, despite this prediction of a weak field, we did not see a largescale, global collapse occur in our simulations. This, we will argue, was due to the kinetic energy of the flows themselves, which also prevented global collapse in the hydro version of the θ = 0 • case (Carroll-Nellenback, Frank & Heitsch 2014), rather than the magnetic field. Thus, the M2FR should be used with caution for estimating global stability within the colliding flows model. We next consider smaller scales associated with the cold, dense phase out of which molecular clouds form. The cold neutral medium (CNM) in our simulations had number densities of approximately n ≈ 500 cm −3 . For the same initial uniform magnetic field, the critical length for the CNM is, As we will see, column density structures associated with the CNM were typically much larger in width than this. Thus, widespread local collapse (i.e. over length scales associated with the cold component of the gas) might be expected to occur in our simulations. However, we did not see widespread local collapse. Clearly, the magnetic field assumption in the M2FR calculation was in error. Turbulent instability behind the shocks would have deformed the magnetic field, thereby producing strong field fluctuations. We suspect that field amplification within the cold, dense gas inhibited collapse locally, and that only after the excess magnetic energy was lost (e.g. through numerical reconnection), could collapse proceed. Thus, we expect the cold clumps were actually largely sub-critical and that this could have been shown through a more rigorous calculation of the local M2FR. While such a calculation was beyond the scope of this paper, others have worked on calculating the local M2FR in simulations (e.g. Banerjee et al. Thus, in what follows we will be concerned with issues of local vs. global collapse and the mechanisms by which other processes, such as, turbulence and field amplification, can inhibit collapse. In particular, we must consider that for large values of the inclination angle θ, post-shock flows can retain a significant fraction of their pre-shock velocity. As shear leads to turbulence, we expect the turbulent velocity to be some fraction of the incoming velocity (v turb ∝ f vo). If turbulence provides support to the cloud against collapse, then we would expect enhancements of the local Jeans length due to turbulence (λ turb ) to be of order, where cs(x) is the local sound speed. Thus, depending on the fraction of inflow velocity retained by the turbulence, we expect the collision region in flows with shear to be more stable to collapse than those without shear. In addition, any degree of turbulence will produce local field amplifications (i.e. additional support against collapse), which as stated, was not accounted for by the M2FR analysis presented above.

PROTOCLUSTER FORMATION AND EVOLUTION
In this section we summarize our principle result by demonstrating how shear and magnetic fields directly affected the star formation in terms of the creation of sink particles. While the sink particles in our simulations had masses that should be associated with clumps (i.e. protoclusters), were we to include higher levels of resolution we would expect the clumps to form cores which would then collapse to form individual stars. Fig. 3 presents the final masses of the sink particles as a function of their formation time for the various runs. This plot relays three pieces of information that show the effect of shear. First, the total number of protoclusters formed in each of the runs decreased with shear. The number decreased from four protoclusers in the θ = 0 • case to one in the θ = 30 • case to zero in the θ = 60 • case (in the run time given to the other simulations). Only after extending the θ = 60 • case out another 5 M yr did a protocluster eventually form at t sink (θ = 60 • ) ≈ 32 M yr, as shown in Fig. 2. While this trend did not hold for the θ = 15 • case, which formed a couple of low mass protoclusters late in the simulation, it is clear that the initial inclination angle of the colliding flows interface was directly related to the amount of post-shock support against collapse.
Second, higher levels of shear delayed the formation of protoclusters. The θ = 0 • case produced a protocluster by t ≈ 11 M yr, whereas the θ = 15 • case did not produce a protocluster until t ≈ 18 M yr. As mentioned above, the (θ = 60 • ) was inhospitable enough that protocluster formation was delayed until t ≈ 32 M yr. This shows the various simulations were evolving under different timescales. This can be understood by considering tcrit -the timescale to acquire enough material into the central region to go magnetically supercritical. As the inclination angle was varied, the velocity in Eqn. 6 would have become v = v cos θ, due to deflection at the collision interface (i.e. the generation of shear). Thus, the critical timescale for shear environments is related to tcrit by, increased time to form protoclusters and the overall reduction in protocluster number. Finally we note that while the θ = 30 • case does disagree with this result (forming its single protocluster before the θ = 15 • case), this sink particle was extremely low-mass despite having ample time to grow. This indicates that while local collapse occurred in the region where the protocluster formed, the shearing motions generated by the initial inclination angle of the interface inhibited further growth. Third, shear slowed growth of protoclusters, i.e. higher shear cases had lower average accretion rates (<Ṁ >= M f inal /(tsim − t sink )). For example, consider the two protoclusters that formed at t ≈ 24 M yr. For this pair, the θ = 0 • protocluster grew to be ≈ 3× more massive than its θ = 15 • counterpart. This indicates that the environments surrounding the protoclusters were less bound in the higher shear runs (recall, only gravitationally bound and unstable gas in the surrounding zones can be accreted onto sink particles). The higher accretion rates of the lower shear simulations translated into more massive protoclusters.
The role of magnetic fields on dynamics can also be extracted from Fig. 3, as there was a significant reduction in the number of protoclusters formed in the MHD, no-shear case (θ = 0 • ), compared to the hydro version of this run (Carroll-Nellenback, Frank & Heitsch 2014). For the same final simulation time and resolution, the hydro run formed a total of 27 protoclusters compared to the 4 formed in the θ = 0 • case. Moreover, Carroll-Nellenback, Frank & Heitsch (2014) found differences in the mass distributions of protoclusters depending on whether there was global or local collapse. The present protocluster masses are consistent with the mass distribution of the hydro run, which only exhibited local collapse. Taken together, this demonstrates the local support magnetic fields are providing against collapse, since . Masses of the various sink particles at 1 M yr after their formation with corresponding final masses. Each set of points (early and final mass), is connected by a line. The θ = 0, 15, 30, 60 • cases are given by the blue, black, red, and grey lines, respectively. Note that the last sink particle for the θ = 15 • case is not plotted, as that sink's lifetime was less than 1 M yr. The 60 • is a single point as that sink lived for exactly 1 M yr.
neither the hydro or MHD runs showed evidence of global collapse.
We now turn to comparing the sink particles at a similar evolutionary time. Figure 4 shows the mass distribution of the sink particles at 1 M yr post-formation. This time was chosen as it maximizes the accretion time of the θ = 60 • sink. However, this time loses the final sink of the θ = 15 • run, which had a lifetime < 1 M yr. In the figure, each sink's mass at 1 M yr post-formation is connected by a line to the final mass of that sink. This provides a quick reference of the average accretion rates of the sink particles. Figure 4 shows that both the initial peak mass, as well as the average initial mass, decline with increasing shear (focusing on the leftmost points of each pair). It is interesting that the θ = 60 • case seems to contradict this behavior. However, this has to do with the reorientation of the collision interface (discussed in more detail below). By the time this sink particle forms, the initially steep collision angle had become largely normal to the oncoming flows. This would have reduced the degree of post-shock shear, and consequently the amount of support afforded to the surrounding envelope. Thus, the θ = 60 • sink particle was able to accrete at a rate comparable to lower shear runs. Mass-weighted histograms of the vorticity in the collision region for the θ = 60 • case shows that this has indeed occurred (Fig. 5).
To conclude this section, we note that the accretion rates of the various sink particles were highly environmentally dependant. Some sinks were positioned in large regions of gravitationally unstable, collapsing gas compared to others, and this would have caused their accretion rates to be larger. This is why some of the sinks that have formed later in time have grown to be more massive than older sinks.

MORPHOLOGY
We now turn to column density maps (CDMs) for a more detailed comparison of the flow evolution. We begin by discussing a morphological artifact common to MHD colliding flows, and then move on to the main morphological features of the flow.
Column density maps of the θ = 0 • case are shown in Fig. 6. As can be seen in the left hand panel, there is a large, ring-like structure surrounding the flows. This ring is an artifact of the simplified initial conditions and its formation can be understood as follows. Initializing colliding flows as cylinders naturally produces a region where shocked gas is expelled laterally with respect to the cylindrical axis. Any colliding flows geometry will produce a characteristic two shock structure (one to decelerate each flow), separated by a contact discontinuity. For finite-sized flows, there must also be a region where high pressure post-shock material is driven out of the collision region. Analysis of a similar process in time variable protostellar jets (which, with a change in reference frame are similar to the configuration studied here), shows that lateral motions on the order of the post-shock sound speed (cps) carry material away from the interaction region along radial streamlines into the ambient gas (Raga 1992). In this way, converging flows along the length of the cylindrical regions are converted into radial flows expanding away from the axis of the cylinders.
When magnetic fields are present, tension forces can restrict these lateral motions. The length and time scales for this restriction depend both on the field strength and geometry. Studies of magnetized time variable jets with strong cooling show that even initially weak toroidal fields can lead θ = 0 • Figure 7. Streamline plot of Bx and Bz components of the magnetic field averaged along y as shown in the x − z plane for the θ = 0 • case. Note the bending of the field lines as material is expelled from the collision region (see text for description).
to the collapse of post-shock flows onto the axis (De Colle, Raga & Esquivel 2008;Hansen, Frank & Hartigan 2015). Thus, the ring-like structure present at approximately 30 pc away from the center of the colliding flows can be attributed to the effect of magnetic tension, since it was not seen in our pure hydro case (Carroll-Nellenback, Frank & Heitsch 2014). Such rings were also seen in other MHD colliding flows runs (Vázquez-Semadeni et al. 2011;Banerjee et al. 2009). We note that in the simulations presented here, the initial field was parallel to the colliding flows. Thus, as post-shock gas was driven outward from the interaction region, the flow drove arcs in the magnetic field (Fig. 7) whose tension eventually halted further expansion, thereby producing a ring of high density material in the collision plane. The position of this ring can be estimated by assuming the ring has reached a steady-state, as we will do next.

Magnetized Ring Model
To make calculations simplest, we envision the following scenario. At t = 0, the magnetic field is B = B0x for r > R, where R is the colliding flows radius. For r < R, we take B = 0. This is a fine approximation, given the field is dynamically weak within the flows (recall, βram ≈ 38). As material enters the collision region, it is shocked and then expands away from the collision region, as described above. We approximate this expansion as being spherically symmetric. Now, the ram pressure of the ejecta pushes outward on the surrounding low density, magnetized ambient medium. In 2D, this leads to a 'ring' of flux that moves outward (in 3D, a spherical 'shell' ; Fig. 8). This ring has two boundaries, an outer radius, ro, and an inner radius, ri (Fig. 8, bottom Figure 6. Column density map for the θ = 0 • case. The integration length is the 62.5 pc for the left-hand column and 75 pc for the right-hand column in these and all subsequent column density maps. Sink particles are given as fuchsia points. Each tick mark represents 10 pc. panel). Once this ring comes to steady state, the magnetic pressure at ro will be balanced by the unperturbed, ambient magnetic pressure, Pmag, which is ∝ B 2 0 . At ri, Pmag will be balanced by the ram pressure of the ejecta. We approximate this ram pressure to be the pre-shock ram pressure, given the low Mach of the flows. Note, we will ignore edge effects near the colliding flows themselves, and instead focus on the dynamics of the ring perpendicular to the flows (Fig. 8, top  panel).
Assuming the magnetic field is the dominant force in the Figure 8. Diagram of the magnetic ring model. Top panel shows the incoming flows, spherical post-shock expansion region and corresponding magnetic field arcs, and the direction along which we are integrating the momentum equation (r). Also given is the unperturbed field value (B 0 ), the cylindrical radius of the colliding flows (R), and the normal unit vector (n) discussed in the text. Bottom panel shows the ring of high density material (in the mid-plane of the flows) with inner radius r i and outer radius ro. ambient (therefore, ignoring gravity and pressure forces), the steady-state ideal MHD momentum equation for the ring is given by, where µ is the magnetic permeability, r is the radius of curvature, which we take to just be the distance from the center of the colliding flows to a position in the ring, andn is a unit vector normal to the field line that is anti-parallel to r (i.e. n =−r). Defining the radius of curvature in this way is equivalent to assuming the field is being bent along spherical arcs. Equation 10 says that in steady state, the magnetic tension of the curved field lines (RHS) is balanced by the gradient in magnetic pressure (LHS). Projecting Eqn. 10 onto the r-axis and integrating gives the following expression for B 2 (r) in the ring, Using this, we balance the magnetic and ram pressures at the ring's inner radius ri, which reduces to, Now, the inner radius is where ejected material is being decelerated, so we wish to find an estimate for this edge of the ring to compare with our simulations. For this, we will use flux-freezing. Equating the flux in the ring at steady state, to the flux in the ring initially, gives, Plugging in for ro (Eqn. 13) yields, which is approximately, or r ≈ 3R. As seen from the left hand panel of Fig. 6, this simple model reproduces the ring's position to within a factor of two.

Main Features in Column Density
Moving past considerations of the ring, we now focus on the details of the flow in the collision region. Note that the incoming flows corresponded to WNM, which was at a column density of N ≈ 2 × 10 20 cm −2 . After passing through the shocks, the gas entered the collision region where the thermal instability drove the gas into the CNM state, which the pdfs of Section 8 show was ≈ 500× denser. We identify the different phases in column density loosely, based on morphology and relative mass fraction from previous studies (Heitsch et al. 2007;Hennebelle et al. 2008;Banerjee et al. 2009). In addition, we use the minimum HI column density for UV-shielding (NHI ≈ 1 − 2 × 10 21 cm −2 ) to label regions that have effectively become 'molecular' (van Dishoeck & Black 1988; van Dishoeck & Blake 1998). Note that to have accurately tracked molecular gas, we would have needed to include UV shielding in the code. In the images, dark blue regions are the initial WNM, cyan regions are thermally unstable gas, from which the CNM phase and subsequent molecular clouds form (green-yellow), and denser gas (i.e. clumps) within the molecular phase is shown in orange and red.
The structures seen within the interaction region of the θ = 0 • case (Fig. 6) were morphologically similar to those of the hydrodynamic version of these runs (Carroll-Nellenback, Frank & Heitsch (2014), fig. 2). We note that the degree of heterogeneous or "clumpy" structure appears to be lower in the present magnetized run. This impression is in line with previous findings for other MHD studies (Heitsch et al. 2007;Hennebelle et al. 2008;Heitsch, Stone & Hartmann 2009;Chen & Ostriker 2014), which find that fields tend to produce larger, more coherent filamentary structures, compared to the hydro versions of these simulations.
The first frame of Fig. 6, taken at t = 10 M yr, shows the gas shortly before the first protoclusters formed. By the next frame (t = 20 M yr), three protoclusters had formed close enough to each other that two of them overlapped and nearly overlapped the third (as can be seen by zooming in on the figure). The location of these sink particles was offcenter from the flow axis, indicating they did not form out of a global collapse mode. Instead, they were within local potential minima, associated with regions of high density gas that condensed out of the background turbulent environment and had become gravitationally unstable. By the last time panel of Fig 6, another protocluster had formed offaxis, roughly 25 pc away. This sink particle was in a large region of high N . Given this sink had the fastest average accretion rate of the θ = 0 • run protoclusters (Fig. 3), this region represented a large and deep potential well. Finally, we note that there was significant widening of the collision region with time, due to the NTS, KH and cooling instabilities.
We now shift our attention to column density maps of the θ = 15 • case ( Fig. 9). At t = 10 M yr, the collision region was structurally similar to the 0 • case, except for regions of lower maximum N (see the left hand column of the plot). Thus, even at low inclination angles, the shear generated by the oblique shocks at the interface disrupted substructure formation to some degree. However, the flows were capable of assembling some localized structures dense enough to become molecular by this time.
By t = 20 M yr, the first protocluster can be seen, having formed at t sink ≈ 17 M yr. It is located near the axis as seen in the left hand panel of Fig. 9, but off-center, in one of the NTSI nodes as seen in the right-hand panel. Over time, the fingers of the NTSI grew, and by t = 27.3 M yr the instability had become 'z-shaped', with a large component aligned with the flow axis. Also by this time, four additional protoclusters had formed. Two of these were in the same NTSI node mentioned above, and the other two were roughly 10.2 and 21.9 pc away, near the edge of the cylinder. The positions of all of these sink particles was again consistent with local instability being triggered in the flow.
Finally, we discuss the extreme θ = 60 • case, which exhibited very different behavior than the other runs. In particular, we found that the initial steep angle of the collision interface reoriented over time and became normal to the incoming flows. Animations 3 of the runs suggested the reorientation was due to perturbations at the collision interface seeding the NTSI. However, field line tension may have also played a role, as a similar, albeit weaker realignment is seen in the hydro version of this run (Haig & Heitsch in prep). The reorientation appears to begin in a region where the flow was nearly normal to the surface of an NTSI node (see the right panel of Fig. 10 at t = 10.1 M yr). Because of the high angle of the collision interface, most incoming material was deflected by the oblique shocks to flow parallel along the shock face. The NTSI node created a local distortion of the shock, allowing oppositely driven material to meet at lower obliquity. The gas behind this stronger shock had higher thermal pressure and began expanding in the (y,z) plane. This expansion increased the area of the low obliquity regions of the flow, which led to more high pressure post-shock gas as seen by t = 20.1 M yr (Fig. 10, right). In this way, it appears the initially highly oblique shock region was transformed into a shock that was more normally directed, relative to the incoming flows. By t = 32.8, the shock had become nearly fully normal.
We conclude with comments on the general characteristics of the θ = 60 • collision interface. First, the collision region was less dense than in the other runs at t = 10.1 M yr (Fig. 10, left). This suppression of growing high density regions arose from the strong shear. Incoming material was sharply deflected away from the axis of the cylinder, and thus less gas accumulated at early times. Second, the collision interface was thinner compared to the other runs at early times, because of this deflection of material away from the interaction region. Material streamed out into the ambient medium, rather than building up along the flow axis. Only by t = 20 M yr, given the reorientation of the shocks bounding the interaction region (and corresponding weaker degree of shear), did material begin to accumulate in the collision region (Fig. 10, left). By t = 32.8 M yr, local collapse had set in and a single protocluster is visible on the CDM, again having formed away from the global potential minimum.

MAGNETIC FIELDS AND DYNAMICS: β −1 MAPS
We next discuss maps of average β −1 , which were generated by discretely summing β −1 through the grid along the same lines of sight used in the column density maps of the previous section, i.e.,β where L, mx, β −1 i , and dxi are the length of the box along the given dimension, number of cells along that dimension, value of β −1 at the ith cell, and the ith cell's width, respectively.
Beginning with the θ = 0 • case, field amplification associated with the radial ejection of gas from the collision Figure 9. Column density map for the θ = 15 • case. Sink particles are given as fuchsia points. Each tick mark represents 10 pc. region (discussed in Section 6) produced a ring of high average β −1 at r ≈ 30 pc (Fig. 11, left). Additionally, there was another, inner ring present from t = 10.1 M yr on at a distance of r ≈ 15 pc from the center of the colliding flows that did not appear in the CDMs. This ring was not contained in the collision region itself, as can be seen from the right-hand frames of Fig. 11. Rather, it stretched down the length of the cylinder, peaking in brightness near the outer edge. This indicates it was generated by the strong shear present at the boundary between the colliding flows and the stationary ambient medium. This shear layer would have caused strong cooling and therefore a decrease in thermal pressure (hence the increased β −1 along this boundary).
Looking down the barrel of the flows at t = 10.1 M yr (Fig. 11, left), there are regions inside of the collision region where average β −1 had increased above its initial value of Figure 10. Column density map for the θ = 60 • case. The single sink particle is given as a fuchsia point. Each tick mark represents 10 pc. Figure 11. Average β −1 map for the θ = 0 • case (see Eqn. 19 for definition). Sink particles are given as black points. Each tick mark represents 10 pc. β −1 = 0.1. A glance back at Fig. 6 shows these regions are co-located with regions of high column density N . In some of these regions, average β −1 had increased by a factor of 10 or more. Since this material is associated with the CNM phase (as discussed in Section 6), its thermal pressure was at least equal to its initial value (cf. Section 8). Thus, in order for β −1 to have increased by this amount, the magnetic pressure must have increased by a factor of 10 or more. This supports that the field was being strongly amplified in the high density gas.
In addition, we see regions where average β −1 was reduced below its initial value. These 'voids' (shown in the lowest color on the color bar) generally correspond with regions of low N . Since flux freezing implies the accumulation of flux is co-extensive with regions of high density, the voids were formed as material was swept away by both turbulence and the collapse of gas into neighboring potential minima. In other words, voids are associated with regions that have a net positive ∇ · v. Thus, regions of low N should not be strongly magnetized, consistent with the voids in Fig. 11. As the simulation evolved, there was an increase in the size of both the magnetic voids, as well as, regions of high average β −1 .
The ramifications of magnetic fields on the dynamics of the flow were already discussed in Section 5. The presence of high average β −1 regions co-located with regions of high N support that magnetic fields suppressed star formation locally. This is further supported by the location of forming protoclusters. As seen in the t = 20.1 and 27.3 M yr panels, protoclusters formed away from regions of highest average β −1 . That is, protocluster formation was inhibited where the field was the strongest. Instead, protoclusters formed where the gas was dense, but average β −1 0.6.
Average β −1 maps of the θ = 15 • case (Fig. 12) are similar to the case with no shear. Regions of high average β −1 are associated with high N structures. There was also a similar formation of magnetic voids. The shear angle present in these runs does, however, leave a signature in the average β −1 maps early in the simulation. Looking down the axis of the flows at t = 10.1 M yr, it is evident that there are weaker regions of enhanced average β −1 (Fig. 12, left). This likely occurs because of the way incoming flows were redirected away from the axis, due to the oblique shocks. Whereas material was driven away from the collision region with radial symmetry in the no-shear case, material here picked up positive and negative vz components across the contact discontinuity. The field lines threading the interaction region must have therefore undergone a local stretching and may have been shorted out by numerical diffusion. This is in contrast to the large scale arcs of field which formed in the purely radial flows of the θ = 0 • case. Fig. 13 shows a very different set of average β −1 maps for the θ = 60 • case. First we note that projected along the axial view we see almost no regions of high average β −1 . This is because the interaction region remained quite thin early on as gas passed through the oblique shock and was quickly shunted away from the central region. When seen from the side, however, (Fig. 13, right) we do find significantly higher values of average β −1 in the collision region. This reflects both the projection through the thin interaction region, as well as the strong local field amplification that occurred due to the shear. Field lines must have assumed a "z" shaped configuration as gas moving from the left on one side of the contact discontinuity was driven upward by the oblique shock it encountered, and gas moving from the right on the other side of the contact discontinuity was driven downward by the oppositely oriented oblique shock it encountered.
Note that as the collision region reoriented from t = 20.1 M yr onward (cf. Section 6), more material and field collected inside of the interaction region in the same manner as was discussed in the lower shear cases. This accounts for the increased regions of high average β −1 seen in the left hand side of Fig. 13 at later times, although overall these regions were much smaller than their lower shear counterparts.

THERMODYNAMICS -PROBABILITY DISTRIBUTION FUNCTIONS
We now discuss the thermodynamic evolution of the gas as it was first shocked in the flow collision and then cooled into a cold, dense state, which could then undergo further compression (or expansion) due to gravity and/or magnetic fields. The probability distribution functions (pdfs) in this section give the amount of mass at a given pressure and number density. The y-axis gives both scaled magnetic and thermal pressure (marked by the color and gray scale, respectively), and the x-axis gives number density n (cm −3 ). Isotherms are straight lines of slope m = 1 on these pdf log-log plots, and increase in temperature from right to left. We begin with the θ = 0 • run. The thermal pressure distribution (i.e. grey-scale pdf, Fig. 14) in this and the remaining runs was identical to previous colliding flows studies. Namely, gas was shocked and heated as it entered the collision region and via the thermal instability, cooled down until it reached the equilibrium curve. For a more detailed discussion of the dynamics relayed by the thermal pressure pdf, we refer the reader to previous work (e.g. Carroll-Nellenback, Frank & Heitsch (2014)). Presently, we consider the simultaneous thermal and magnetic evolution of the gas as shown by overlaying the Pmag, n pdf in color. Note, however, it is not always possible to exactly correlate values in one pdf with values in the other (i.e. to simultaneously know the P therm of a given Pmag, n combination).
As can be seen in Fig. 14, most of the mass at low densities (i.e. near log n = 0) had β > 1, consistent with the notion of magnetic voids discussed in Section 7. This follows from comparing the Pmag distribution with the P therm distribution at these densities. While the distribution of Pmag had a much higher spread, this spread was mostly below the P therm distribution. In contrast, gas at higher densities (log n > 2.5) was mostly on top of the equilibrium curve, implying most of the densest gas had β < 1. This is in agreement with the densest gas having had enhanced magnetic support (and thus greater support against collapse), as was also discussed in Section 7. Despite the spread in Pmag at these high n, the gas was entirely at T = 10K, as the P therm distribution was constrained to lie on the equilibrium curve at these densities. Lastly, the Pmag, n pdf tended to shift toward the upper right of the plot while becoming more tightly distributed. If a line were drawn through the middle of the entire distribution, it would be roughly linear until about log n ≈ 1, and would then shift upwards, which Figure 12. Column β −1 map for the θ = 15 • case. Sink particles are given as black points. Each tick mark represents 10 pc. Figure 13. Column β −1 map for the θ = 60 • case. The single sink particle is given as a black point. Each tick mark represents 10 pc.
on this log-log plot would correspond to a power law relationship similar to what is observed in other work (Heitsch et al. 2007;Hennebelle et al. 2008;Banerjee et al. 2009). Now we discuss the dynamics attributable to these various features of the Pmag, n pdf. As the simulation proceeded, turbulence at the collision interface (produced by the KH and NTS instabilities) generated magnetic field fluctuations. Reconnection diffusion, effective at removing excess magnetic energy from magnetized, turbulent flows (Lazar-ian, Santos-Lima & de Gouveia Dal Pino (2010); see also Klessen, Heitsch & Mac Low (2000); Federrath et al. (2011)), created the large (i.e. many orders of magnitude) spread in Pmag at low n. This spread in Pmag persisted as the gas was cooled and compressed. Flux-freezing, which accompanied this compression, then led to a sharp increase in Pmag, as seen by the rather large island of material (in red shading) up and to the right of the initial values. Note that colors correspond to a log scale, and so the the amount of mass Figure 14. Pressure vs. density histograms for the θ = 0 • case. The y-axis gives both the thermal (grey scale distribution) and magnetic (color scale distribution) pressure as a function of number density (the logs of these quantities, that is). The color bar gives the amount of mass at a given pressure and density. The thermodynamic equilibrium curve is given on the plot, as well as the initial magnetic pressure of the flows (bottom horizontal line), ram pressure of the flows (top horizontal line), and T = 10 K isotherm (diagonal line).
within this population was significantly greater than any other Pmag, n combination of the flow.
As time continued, a high Pmag, n "tail" of the pdf formed that, 1) followed the equilibrium curve, 2) was above the ram pressure line of the flows, and 3) only appeared after protoclusters had formed in the flow. Thus, this tail traced gravitational instability in the flow. Most of the mass in the tail was at a β < 1, consistent with the discussion from Section 7 that regions of high density were correlated with regions of high β −1 . Additionally, there was a spread in β for the population of gas parcels in this tail, as was also discussed in Section 7. In particular, protoclusters formed in dense enough gas that was at lower average β −1 than the surroundings. Numerical reconnection within this high Pmag, n tail reduced the local field strength, thus enabling protocluster formation.
The pdfs for the θ = 15 • case (Fig. 15) were similar to those for the θ = 0 • case. The largest difference was the delayed growth of the high Pmag, n tail in the higher shear run. Thus, the plot offers further evidence that shear suppressed compression in the flows. Eventually collapse was triggered in localized pockets of the flow, allowing gas parcels to begin to climb the equilibrium curve, as seen by t = 20.1 M yr. As in the previous case, this coincided with the presence of newly formed protoclusters in the flow.
Keeping in mind the trends discussed previously, Fig.  16 shows an even longer delay in collapse/compression in the highest shear angle case, θ = 60 • . This is evident by the diminished upper-right island of high Pmag, n material early on. It is also shown by the greater delay in the high Pmag, n tail compared to the lower shear runs, which was still not visible by t = 27.3 M yr. Only by t = 32.8 M yr did the high Pmag, n tail emerge in the pdf, again corresponding with the presence of a newly formed protocluster in the flow. Table 2. Analysis region for the different runs. Each were centered on the collision interface and had lengths in x, y, and z, given by Lx, Ly, Lz. The maximum wavelength λmax (∝ k −1 min ) for each run is also given, which was equal to the longest dimension of the analysis region.

ENERGY SPECTRA
We next turn to power spectra of the gas in the collision region for kinetic, gravitational, and magnetic energies (E kin , Egrav, and Emag, respectively), following Carroll-Nellenback, Frank & Heitsch (2014). To accommodate steeper shear angles, the x dimension of the analysis region increased with θ and a 10 pc buffer was added to either side. The dimensions of the analysis region for the different runs are given in Table 2. Spectra for the different runs are shown in Fig. 17, for three different times. Note, the x-axes in the spectra are equivalent to λmaxλ −1 , where λmax is given in Table  2. Thus, each unit on the x-axis represents that fraction of λmax. We include a k −2 line (solid line, red) to compare to the spectra of v 2 (dotted-dashed line, black) for t = 10.1 M yr. We do not include the v 2 spectrum for other times, as it was largely unchanged throughout the course of the simulations.
As can be seen by the v 2 spectrum, gas displayed the Burger's turbulent spectrum (v 2 ∝ k −2 ) in all of the runs over some range of k. For the θ = 0 • case, this was 3 < kk −1 min < 40, corresponding to length scales of roughly 13 pc > λ > 1 pc. At the low k end of the spectra, the driving scale of the turbulence is apparent and was on the order of the colliding flows radius (λ ≈ 40 pc). At higher k, both dis- Figure 15. Pressure vs. density histograms for the θ = 15 • case. As before, the y-axis gives both the thermal (grey scale distribution) and magnetic (color scale distribution) pressure as a function of number density (the logs of these quantities, that is). The color bar gives the amount of mass at a given pressure and density. The thermodynamic equilibrium curve is given on the plot, as well as the initial magnetic pressure of the flows (bottom horizontal line), ram pressure of the flows (top horizontal line), and T = 10 K isotherm (diagonal line). Figure 16. Pressure vs. density histograms for the θ = 60 • case. As before, the y-axis gives both the thermal (grey scale distribution) and magnetic (color scale distribution) pressure as a function of number density (the logs of these quantities, that is). The color bar gives the amount of mass at a given pressure and density. The thermodynamic equilibrium curve is given on the plot, as well as the initial magnetic pressure of the flows (bottom horizontal line), ram pressure of the flows (top horizontal line), and T = 10 K isotherm (diagonal line). Note, the pdfs were similar enough between t = 20.1 and 27.3 M yr that we do not include the t = 27.3 M yr panel in the figure. sipation and gravitational collapse limited the inertial range captured in the grid.
The dominant energy on all size scales for all of the runs was E kin (smooth, black lines), and remained so throughout the course of the simulations. This was due to the power being generated in the colliding flows themselves. We note that in the hydro version of the θ = 0 • case (Carroll-Nellenback, Frank & Heitsch 2014), E kin > Egrav at late times in the flow on large scales (λ > 10 pc). This indicated a lack of global collapse. The much stronger suppression of Egrav (fine-dashed lines) in the MHD runs (i.e. E kin >> Egrav) shows that large-scale global collapse was also not occur-ring in the MHD cases. Also, E kin < Egrav on small scales (λ < 10pc) in the hydro case. This indicated strong local collapse was occurring. That E kin > Egrav on small scales in the MHD cases further supports that magnetic fields have impaired local collapse. A diminished degree of local collapse occurred in the MHD runs, and only in regions where Egrav > Emag.
Indeed throughout the runs, Egrav was comparable to Emag, on most scales. However, in the θ = 0 • case, Egrav was generally higher than Emag for 10 < kk −1 min < 100. Note, this corresponds to a length scale of .4 < λ (pc −1 ) < 4, which is the size scale of protoclusters. This is consistent with weak, localized collapse, due to a stronger gravitational field on these scales. By t = 27.3 M yr, there was a slight increase in Emag relative to Egrav, which may reflect an enhancement of the field due to compression from gravitational collapse.
The θ = 15 • spectra show similar behavior to the θ = 0 • case. At intermediate scales, Emag and Egrav were comparable. At smaller scales, as gravitational collapse set-in, Egrav > Emag. The effect of shear in this case did not make itself readily apparent in the energy spectra.
The spectra for the θ = 60 • case, however, show clearly the effects of imposed shear on the distribution of energies at different scales. At 10 M yr, Emag > Egrav, on all size scales. In the range of 9 < kk −1 min < 90, corresponding now to size scales of 1 < λ(pc −1 ) < 10 (cf. Table 2), we see Emag >> Egrav, with Emag approaching E kin near kk −1 min ∼ 60 (λ ∼ 1.5 pc). Note that this is the only run in which such an equipartition occurred. This enhanced Emag relative to E kin occurred with a simultaneous decrease in power of Egrav on all scales compared to the other runs. Thus, the signature of high shear was a relative amplification of the field on protocluster scales, as well as a decrease in Egrav on all scales, which we attribute to the disruption of dense structures forming in the flow as already discussed. Only by the end of the simulation (t = 32.8 M yr) did Egrav become comparable with the energy locked up in the magnetic field (which was itself comparable to E kin ). Thus, the spectra are consistent with the formation of a sink particle at late times.

DISCUSSION
We have presented 3D, adaptive mesh refinement simulations of magnetized, shear colliding flows including selfgravity, sink particles, and cooling. Feedback was not included in the present work and thus our simulations only followed the formation and early evolution of molecular clouds. The simulations were of two colliding flows, intersecting at an inclined interface under a dynamically weak magnetic field. The inclined interface was used to generate shear at the collision layer, and was varied from not-inclined (nor-mal incidence of the colliding flows) to highly inclined. As molecular clouds are unlikely to form from perfectly head-on collisions between two large-scale streams of gas, breaking the symmetry of the collision interface in this way was a previously unexplored next step in the colliding flows model exploration. The purpose of this set of experiments was to study the effect of shear on molecular cloud formation in magnetized colliding flows.
Our work has shown that shear, imposed by an inclined collision interface, impacts cloud dynamics and reduces protocluster formation. In particular, as shear increased (i.e. the inclination angle steepened), it took longer to form protoclusters, they formed in lower numbers, and they were lower in mass due to diminished accretion rates. These effects are consistent with recent work by Körtgen & Banerjee (2015), who found that higher degrees of inclination between colliding flows leads to a reduced number of sink particles, as well as a delay in their formation. Additionally, Chen & Ostriker (2014) show that under ideal MHD conditions, increasing the angle between upstream colliding flows and the magnetic field (which, with a change of reference frame is similar to our setup here) decreases the number of gravitationally bound cores that form in the post-shock region.
We also showed that without shear, protocluster formation was greatly impeded in the presence of even a weak field (β = 10, βram ≈ 38) alone. This was evidenced by the stark difference in protocluster number between the noshear, MHD case (θ = 0 • ) and the hydro version of this simulation in Carroll-Nellenback, Frank & Heitsch (2014). More than 5 times as many protoclusters formed in the hydro run (n=27) compared to the MHD run (n=4). This result fits in with other colliding flows studies that have shown magnetic fields impair gravitational collapse. Hennebelle et al. (2008) showed that it takes longer to form self-gravitating objects in MHD colliding flows (aligned field, B0 = 10 µG), compared to flows without a field. Heitsch, Stone & Hartmann (2009) also compared colliding flows with and without a magnetic field and found that the degree of post-shock turbulence decreases with a magnetic field. This led to a suppression of dense cores forming in the magnetized cases. Addition-ally, Vázquez-Semadeni et al. (2011) found that increasing the magnetic field strength in colliding flows simulations decreases the SFE of forming clouds.
To understand the flow dynamics responsible for the weaker protocluster formation, we began by distinguishing between 'global' and 'local' collapse, and showed that collapse in our simulations was due to localized regions becoming unstable, rather than the entire collision region. Indicators of local collapse included the positions of sink particles away from the global potential minimum, weak mass 'fall-back' into the collision region, and power spectra that showed E kin > Egrav on large scales, but E kin < Egrav on small scales (Carroll-Nellenback, Frank & Heitsch 2014). Our results were consistent with these indicators: sinks formed away from the center of the collision region, they had variable accretion rates, and both Egrav and Emag were << E kin on large scales, but Egrav ∼ E kin > Emag on small scales. All of the runs (θ = 0−60 • ), as well as the hydro run, exhibited local collapse only (i.e. did not collapse globally).
Compared to the hydro run, collapse was weakened in our θ = 0 • case due to the interplay between the turbulent post-shock velocity field and the magnetic field. Turbulence in the interaction region introduced distortions to the magnetic field, which produced regions of local field amplification. These were shown to be co-located with regions of high density (Sections 7 and 8). Thus, local collapse had an additional form of support in the MHD case(s) -magnetic support, due to field amplification (see also Banerjee et al. (2009) ;Heitsch, Stone & Hartmann (2009)). Over time, numerical reconnection would have reduced the field strength in these pockets of dense gas (Lazarian, Santos-Lima & de Gouveia Dal Pino 2010;Vázquez-Semadeni et al. 2011;Federrath et al. 2011;Chen & Ostriker 2014), thereby allowing highly localized regions to become magnetically supercritical and produce protoclusters (should the gas also be Jeans unstable).
As shear increased in the flows, a higher degree of postshock turbulent velocity resulted. This disrupted the formation of dense structures in the flow, which further inhibited local collapse and protocluster formation (by locally decreasing the M2FR, as well as, increasing the Jeans length). This seemed to be the dominant effect of shear, as we did not find that stronger shear led to greater field amplification. The reasons for this were likely two-fold. While the shear generated by the oblique shocks at the collision interface would have led to greater distortions of the field (Hartmann 2001;Heitsch et al. 2007;Chen & Ostriker 2014), they were also likely shorted out faster in the higher shear cases. Additionally, given the disruption of high density structures forming in the flow, regions of increased field strength due to flux freezing would have also been diminished. Both the pdfs and spectra (Sections 8 and 9, respectively) support the conclusion that field amplification did not increase with shear. They do, however, show decreases in high density structures. This is consistent with Körtgen & Banerjee (2015), who attribute declining sink particle formation to shear disrupting the formation of high density structures.
Lastly, we discussed the large-scale realignment of the collision interface that took place in our highest shear angle case, θ = 60 • . This unexpected result may have arisen from the orientation of an NTSI node with respect to the oncoming flows, as discussed in Section 6. However, tension in the magnetic field may have also contributed to the realignment. This is supported by current work by Haig & Heitsch (in prep), which shows that the steep 60 • inclination angle reorients to a lesser degree without a magnetic field. A realignment of the collision interface was also mentioned in Körtgen & Banerjee (2015).

ACKNOWLEDGMENTS
We thank Ralph Klesson, Simon Glover, Mordecai-Mark Mac Low, Eric Mamajek, Alice Quillen, and Andy Pon for interesting discussions on the present work. We thank Christoph Federrath and Clare Dobbs for their valuable comments and suggestions, and Bastian Körtgen for proposing including an early mass distribution of protoclusters. We thank the referee, Robi Banerjee, for his insightful review and recommending a vorticity analysis. We thank the University of Rochester's Center for Integrated Research Computing (CIRC) for time on their supercomputers, data storage, and visualization at the VISTA collaboratory, and the Texas Advanced Computing Center and XSEDE for computing time. Lastly, we thank Rich Sarkis of the University of Rochester for maintaining our data systems. This work was supported by the U.S. Department of Energy through grant GR523126, the National Science Foundation through grant GR506177, and the Space Telescope Science Institute through grant GR528562. This paper has been typeset from a T E X/ L A T E X file prepared by the author.