The properties of clusters, and the orientation of magnetic fields relative to filaments, in magnetohydrodynamic simulations of colliding clouds

We have performed Smoothed Particle Magneto-Hydrodynamics (SPMHD) calculations of colliding clouds to investigate the formation of massive stellar clusters, adopting a timestep criterion to prevent large divergence errors. We find that magnetic fields do not impede the formation of young massive clusters (YMCs), and the development of high star formation rates, although we do see a strong dependence of our results on the direction of the magnetic field. If the field is initially perpendicular to the collision, and sufficiently strong, we find that star formation is delayed, and the morphology of the resulting clusters is significantly altered. We relate this to the large amplification of the field with this initial orientation. We also see that filaments formed with this configuration are less dense. When the field is parallel to the collision, there is much less amplification of the field, dense filaments form, and the formation of clusters is similar to the purely hydrodynamical case. Our simulations reproduce the observed tendency for magnetic fields to be aligned perpendicularly to dense filaments, and parallel to low density filaments. Overall our results are in broad agreement with past work in this area using grid codes.


INTRODUCTION
Recent works have shown that young massive clouds (YMCs) can form through the collision of molecular clouds . Dobbs et al. (2020) showed that YMCs are able to form on timescales of 1-2 Myr, in line with observed age spreads (Longmore et al. 2014). Observationally, there is evidence of cloud cloud collisions in our Galaxy from red and blue shifted CO velocities in molecular clouds along the line of sight, in some cases at the sites of massive young clusters (Furukawa et al. 2009;Fukui et al. 2014;Fukui et al. 2017;Kuwahara et al. 2020). Dobbs et al. (2015) also found in galaxy scale simulations that such collisions of massive clouds, although infrequent, do occur.  carried out a parameter study showing high density, low turbulence and high velocities promote YMC formation. They also determined the properties of clusters which formed, showing that for the cloud masses used (10 4 − 10 5 M ) the properties are E-mail: C.L.Dobbs@exeter.c.uk comparable to lower mass YMCs in our Galaxy (Portegies Zwart et al. 2010), though high velocities led to more elongated clouds and larger cloud radii, at least in the earliest stages of evolution. These previous simulations however are all purely hydrodynamical. Whether such clusters still form when magnetic fields are present, and still have the same properties, is an open question.
A number of studies have examined the effects of magnetic fields in simulations of colliding flows of interstellar gas (Heitsch et al. 2009;Inoue & Inutsuka 2009;Körtgen & Banerjee 2015;Wu et al. 2017Wu et al. , 2020Klassen et al. 2017;Fogerty et al. 2016Fogerty et al. , 2017Zamora-Avilés et al. 2018;Seifried et al. 2020). Heitsch et al. (2009) carry out simulations investigating molecular cloud formation, and show that there is a strong dependence on the initial field direction, with the collision only inducive to producing molecular clouds when the field is parallel to the direction of flow. More recent works have included self gravity and shown the impact of magnetic fields on core and star formation. Körtgen & Banerjee (2015) find that magentic fields delay core and star formation, although Zamora-Avilés et al. (2018) find that star formation occurs earlier with strong magnetic fields, due to suprresion of the non-linear thin shell instability. Sakre et al. (2020) model core formation and suggest that stronger fields provide support to allow the formation of more massive cores (see also Inoue et al. 2018). Sakre et al. (2020) also investigate field direction and find that starting with an initial field parallel to the collision produces a more disordered field compared to when the initial field is perpendicular. Wu et al. (2020) also find that less fragmentation occurs in models with stronger fields, and there are fewer stars. Although some studies now explicitly include sink particles (Inoue et al. 2018;Fogerty et al. 2016;Zamora-Avilés et al. 2018;Fukui et al. 2020), few investigate the role of magnetic fields on cluster formation.
Filaments are widespread both in the neutral ISM and star forming regions, and as such many of the above studies also investigate the relation of magnetic fields to filaments. Recent observations now reveal the alignment of magnetic fields with structures in the ISM. In particular, observations appear to show that the magnetic field is typically aligned parallel to filaments in HI (McClure-Griffiths et al. 2006;Clark et al. 2014;Planck collaboration XXXII 2016), and low density molecuar gas (Heyer & Brunt 2012;Heyer et al. 2020). Whereas in higher density molecular gas, the field is more likely to be perpendicular to the filaments (Alves et al. 2008;Heyer & Brunt 2012;Planck collaboration XXXV 2016). Fissel, et al. (2019) show examples of both parallel and perpendicular alignment within the Vela cloud, traced by 12 CO and 13 CO respectively. Simulations of turbulence (Soler et al. 2013;Klassen et al. 2017;Xu et al. 2019) and shock compressed layers (Chen et al. 2016) also show a tendency for the magnetic field to be aligned perpendicular to the field at high density, and parallel otherwise. The dependence of the field orientation may be simply a density criterion (Soler et al. 2013;Chen et al. 2016), or additionally related to the mass to flux ratio (Seifried et al. 2020). Inoue & Inutsuka (2016) find that the alignment of the magnetic field with the filament depends on the initial angle between the shock wave and the magnetic field.
In this paper we perform simualtions of colliding clouds with magnetic fields, though our focus is on the formation of massive clusters rather than individual stars. In Section 2 we describe our method and initial conditions, and in particular the timestep constraint we apply to ensure the magnetic divergence remains low. We descibe the morphologies of the collisions, star formation rates, the relation of magnetic field to filaments and the properties of clusters formed in Section 3. In Section 4 we compare to previous work, and in Section 5 we present our conclusions.

Details of Simulations
We have performed these calculations using Phantom (Price et al. 2018), which is a publicly available Smoothed Particle Magneto-Hydrodynamics (SPMHD) code. Sink particles are included according to the method described in Bate et al. (1995). Magnetic fields are evolved as the magnetic variable B/ρ. Stability of the magnetic fields is ensured using the source term correction (Børve et al. 2001) and the divergence is constrained using hyperbolic divergence cleaning (Tricco & Price 2012;Tricco et al. 2016). Unlike previous work, we apply a modified timestep criteria based upon the divergence cleaning method, which we describe in Section 2.2. This timestep constraint, which is in addition to the usual Courant and acceleration criteria (Price et al. 2018), ensures that the timesteps are small enough to prevent large increases in the divergence. For simplicity, and so that when we vary the velocity field or magnetic field everything else is unchanged, we employ an isothermal equation of state, adopting a temperature of 20 K.
We set up the initial conditions for our simulations in a similar way to Dobbs et al. (2020) and , though with a few differences. We simulate two ellipsoidal colliding clouds, which are colliding head on along their minor axes. The clouds have dimensions of 80×30×30 pc. Both clouds have masses of 1.5 × 10 5 M . All the simulations we present use 6.1 million particles. Although not shown, we initially ran simulations with one tenth the number of particles, which show very similar results to those we display here. The initial setup of the clouds differs from our previous simulations in two main ways. The first is that the two clouds lie within a low density medium, which is one hundredth of the density of the clouds. This is the same approach as Wurster et al. (2019), who modelled isolated clouds within a low density medium. The magnetic field permeates both the clouds and this surrounding medium, preventing the magnetic field becoming unstable at the edge of the cloud, and removing the need for more complex magnetic boundary conditions. For simplicity, the surrounding medium has periodic boundary conditions (satisfying the need for magnetic boundaries), where magnetohydrodynamic forces are periodic across the boundary but gravitational forces are not. Including the low density medium ensures that any increase in the field does not occur at the edge of the simulation region, since the field does not evolve significantly in the low density region. The extent of the low density medium is ±120 pc in the x dimension, and ±46 pc in the y and z dimensions. Secondly, following the setup of the initial conditions in Wurster et al. (2019), the particles are initially allocated on a grid in both the clouds and the low density surrounding medium, rather than randomly.
As for the previous simulations ), we apply a turbulent velocity field to each cloud. The velocity field is set up to follow a Gaussian distribution, which produces a power spectrum consistent with P (k) ∝ k −4 , Burger's turbulence. In our previous work, we showed that quite high velocities are required to form massive clusters over short timescales. Here we only consider one set of collision velocities, and set up each cloud with a velocity of 21.75 km s −1 , such that the total relative velocity between the clouds is 43.5 km s −1 . This velocity is chosen so that the collision has a significant effect on star formation, i.e. the collision enhances the star formation rate above that which would occur for isolated clouds ), but is still consistent with the highest velocities observed for colliding streams in the Milky Way (Motte et al. 2014;Fukui et al. 2015Fukui et al. , 2018. In  we show results with different cloud dimensions and collision velocities, but here we focus on varying the magnetic field strength and orientation. We apply two different field strengths, of 2.5 × 10 −6 and 2.5 × 10 −7 G, and align these fields either parallel or perpendicular to the collision. We note that particularly our weaker field strength is unrealistically low compared to observations, but is intended for comparisons. Our higher field strengths are at the low end of observed field strengths (e.g. Heiles & Troland 2005;Crutcher et al. 2010). We discuss in Section 5 how we expect the trends we observe to extend to higher strengths. The Alfvén velocity is ∼0.06 and 0.6 km s −1 for the weak and strong fields respectively, so significantly lower than both the turbulent velocity field, and collision velocity. We also carry out purely hydrodynamical simulations as well for comparison. We vary the velocity dispersion of the turbulence, which produces clouds which are unbound and bound initially. We show the simulations presented in Table 1. As indicated in Table 1 the main variables in the simulations are the level of turbulence, which changes the virial parameter, the magnetic field strength and the orientation of the magnetic field. We insert sink particles once the density reaches a critical density of 10 −18 g cm −3 and the criteria in Bate et al. (1995) (e.g. converging flows) are fulfilled, using an accretion radius of 0.001 pc. With this resolution, each sink particle typically represents a small group of stars. Artificial viscosity is included with a switch for the α parameter (Cullen & Dehnen 2010). As recommended for strong shocks (Price & Federrath 2010), we take β = 4. The artificial resistivity is described in Price et al. (2018).

Divergence cleaning timestep contraint
The magnetic field is evolved as where ρ is the gas density, B is the magnetic field, v is the velocity and ψ is a scalar field used for divergence cleaning. We assume units of the magnetic field such that the Alfvén speed is vA = |B| / √ ρ (Price & Monaghan 2004). As per Tricco et al. (2016), the evolved cleaning parameter is ψ/c h , where c h is the characteristic speed, referred to as the 'wave cleaning speed.' The evolution of the parameter is given by 1 1 We have modified this slightly from Tricco et al. (2016) to explicily include the overcleaning parameter σ, which has a slightly d dt where c h = v 2 A + c 2 s with cs being the sound speed and h is a scale length (equal to the smoothing length in SPH). To ensure that the cleaning is resolved, a new timestep criteria is introduced. In SPH, the new timestep for particle i, given by where Ccour = 0.3 is the tradional coefficient for the Courant condition. Tricco et al. (2016) introduced the 'overcleaning' parameter σij ≡ σ to control the cleaning. Optimally, σ = 1, however, larger values could be chosen to reduce divergence errors, albeit at the accompying cost of shorter timesteps (recall (3)). The divergence error is monitored by the dimensionless value and they suggest increasing σ if the mean value is > 10 −2 .
In most simulations where the magnetic field is reasonably well-behaved (e.g. Orszag & Tang 1979;Ryu & Jones 1995;Wurster et al. 2019), divB | mean < 10 −2 is satisfied using the default value of σ = 1. However, in very dynamic regions, such as the interface between colliding flows as presented here, this criteria is violated; away from the interface, however, this criteria is satisfied. Therefore, σ > 1 is required for the stability of the magnetic field at the interface. This requires a careful choice of σ prior to starting the simulation such that it is large enough to properly clean the magnetic divergence, but low enough such that computational resources are not wasted.
To circumvent this and to prevent extra computational expense away from the interface, we dynamially calculate different definition here. We explicitly note that this σ is not the velocity dispersion. σij based upon a particle's local environment. Specifically, where σmax is a parameter defining the maximum permitted σij, f is a scalar, and this minimising operation occurs over all of i's neighbours, particles j. Note that we keep a constant value of σ in the wave cleaning equation, (2). The advantage to this method is that we do not need to guess a value of σ prior to starting the simulation, and σij will only increase as needed and where it is needed, which will save computational resources given our use of individual timestepping (Bate et al. 1995). Emperical tests of colliding flows similar to those presented here suggest f = 10 and σmax = 512. Although σij is dynamically calculated, careful consideration must be made of σ and σmax. In the Ryu & Jones MHD shock tube tests (Ryu & Jones 1995), divB | max < 10 −2 meaning that the new algorithm has no impact and it is safe to use the default values. However, in the Orszag-Tang vortex (Orszag & Tang 1979) with 128 particles in the x-direction, the mean value of divB is 0.005, while In this test, setting σmax = 512 has a trivial affect on the results (including the value of divB | max ), except the simulation runs ≈ 10 times slower when performed with global timestepping due to the decreased dt clean ; in this case, it is optimal to set σmax = 1, essentially turning off the dynamic calculation of σij.
Therefore, this new timestep criterion is required when modelling magnetic fields at strong, chaotic shocks (such as the colliding flows presented here) to permit a reliable evolution of the magnetic field. However, it should be disabled when modelling well-behaved magnetic fields.

Morphology of shocked region, clusters and magnetic field
In this section we look at the overall evolution for the collisions of clouds with different magnetic field strengths and orientations. We first discuss the results from the simulations with the clouds with higher virial ratios. The evolution of the BWXHighturb model is shown in Figure 1, top row. The clouds start colliding at around 0.5 Myr. The collision leads to a few main central filaments which are perpendicular to the direction of the collision, as seen in the left panel of Figure 1. The shapes, number and structure of these initial filaments are due to the initial turbulent velocity fields of the colliding clouds. These filaments are gravitationally unstable and as such sink particles form along the filaments. As the collision progresses, a more substantial central structure emerges, the number of sink particles increases, and the distribution of sink particles becomes more clustered rather than filamentary. At the final time frame, 2.4 Myr, the sink particles become particularly concentrated to the uppermost region of the collision interface, leading to a more evident cluster here. The evolution of the other simulations with higher virial parameter is very similar to that shown in Figure 1, with the exception of the run with a strong field perpendicular to the collision (BSYHighturb). The morphology of the collision for these MHD runs is shown in Figure 2, at a time of 2.3 Myr. As shown in Figure 2, the simulations with fields parallel to the collision (aligned along the filament) show very similar morphologies (top row), as does the simulation with a weak field perpendicular to the collision (lower left), although there are some small differences in the morphology of the gas, and the spatial distributions of the sink particles. However the model with a strong field perpendicular to the collision (BSYHighturb) shows a very different morphology.
Here the presence of a filament along the collision interface is less clear, and instead sink particles are strongly grouped into distinct clusters. Most of the sink particles are congregated in a cluster in the upper region where the clouds have collided. The evolution of the star formation in this model is also quite different compared to the others, with fewer stars forming earlier compared to the other simulations.
We show the magnetic field for these models at the same time frame as Figure 2 in Figure 3. As expected, the field is stronger in the models where the initial field strength is higher. However we also see that the field has evolved to higher values in the case where the field is originally perpendicular to the collision. The field is clearly strongest, and has been considerably more amplified in the simulation with the strong field perpendicular to the collision. It seems likely that the high magnetic field in the shocked region where the clouds collide is the reason for the difference in morphology, and the resulting difference in star formation. By contrast the model with a weak magnetic field parallel to the collision (BWXHighturb) shows little amplification of the field in the denser, shocked regions. The field also becomes more disordered in the region of the shock. The field becomes most random in the cases where it is initially parallel to the shock, and the shock appears to increase the component of the field along the shock whereas in the models where the field is already aligned with the shock, the effect is more simply to amplify the field in that direction.
We show the equivalent simulation without magnetic fields, HydroHighturb, in Figure 4 (left panel). The morphology of the gas and distribution of sink particles is fairly similar to the models with weak magnetic fields, and the model with a strong field parallel to the shock, although the sink particles appear more concentrated to one main cluster in the hydrodynamical model compared to the MHD cases. The concentration of sink particles into one cluster is more similar to the model with the strong field perpendicular to the shock, BSYHighturb, although otherwise the morphology of the gas and stars is quite different, and there appears to be much more star formation and many more sink particles in the hydrodynamical model (and indeed the other MHD models) compared to the BSYHighturb model.
We now present results for the models where the clouds have lower virial parameters. In Figure 1 (lower row) we show the equivalent evolution for the simulation with a low magnetic field parallel to the direction of the collision (BWXLowturb). Compared to the case with the higher virial parameter clouds, there is a clearer shocked region, and clearer filaments where sink particles are forming. The morphology of the gas, is not too dissimilar to previous grid code simualtions (e.g. Körtgen & Banerjee 2015). The distribution of sink particles shows a clear elongated structure. In the last panel, the cluster of sink particles appears to have Figure 1. The evolution is shown for the collision of the clouds with the higher virial parameter (BXWHighturb, top panels), and lower virial parameter (BXWLowturb, lower panels). In both cases, the magnetic field is 2.5 × 10 −7 G and parallel to the direction of the collision.
contracted and is less elongated in the direction perpendicular to the collision, due to gravity acting on the sinks and gas. The evolution of this model appears fairly similar to those presented in  with α < 1. Similar to , with a higher virial parameter, the sinks are more dispersed, although in the models here there is still a fairly clear cluster at the location of the shock interface.
In Figure 5 we show the gas surface density plots for the simulations with lower virial parameters and different magnetic field strengths and directions. Again the morphologies of the gas and the sink particles distributions appear similar for three of the simulations, but different for the case with a strong field perpendicular to the direction of the collision (BSYLowturb). Similar to the higher virial parameter simulations, with the exception of BSYLowturb, there is an elongated distribution of sink particles along the shock. By contrast for the BSYLowturb model, far fewer sink particles appear to form, and they tend to be concentrated into a smaller cluster region. The hydrodynamical model, shown in Figure 4 is very similar to the magnetic field models with the exception of BSYLowturb, suggesting that the morphology of these other MHD models is closer to the hydrodynamical case than BSYLowturb.
We do not show the magnetic field vectors for the lower virial parameter models, however the behaviour of the magnetic field is very similar to that shown in Figure 3. The field is strongest in the simulation with a strong field perpendicular to the direction of the collision. There is some amplification of the field for the models with a strong field parallel to the collision, and for the weak field perpendicular to the collision, but the field strength appears relatively unchanged with the weak field parallel to the collision.
In reality, the magnetic field may not be aligned either directly perpendicular or parallel to the collision. We ran a further model where instead the field is aligned at 45 degrees to the direction of the collision. This model shows behaviour that is in between those presented here, i.e. the morphology and distribution of sinks appears in between the cases with a strong field parallel, and perpendicular to the collision, and the magnetic field is amplified to a level in between these two cases.

Star formation rates
We show in Figure 6 the star formation rates for the different models. The top panel shows the star formation rates in the models with the higher virial parameter. The figure shows that the magnetic field does not appear to impede star formation in most of the models, with the star formation rates extremely similar to the purely hydrodynamical case. This is not so surprising since the morphology of these runs is quite similar. The hydrodynamical model has a slight increase compared to the other models but it is not clear this is particularly significant. The model with a strong field perpendicular to the collision (BSYHighturb) however shows quite a different behaviour in the star formation rate. The star formation increases at 0.5−1 Myr later compared to the other models. The star formation rate still reaches values as high as the models though, and actually appears to accelerate faster than the other models after the initial delay. We see from Figure 6 that for the panels in Figures 2-5, shown at a time of 2.3 Myr, star formation has been ongoing for around 1.3 Myr.
In the lower panel of Figure 6 we show the star formation rates for the models with a lower virial parameter. Again Figure 2. The column density, and distribution of sink particles is shown for the collisions with higher virial parameter clouds. The magnetic field is parallel (top row) and perpendicular (lower) to the direction of the collision, and initial field strength is weaker (2.5 × 10 −7 G) in the left panels, and stronger (2.5 × 10 −6 G) in the right hand panels. the star formation rates are very similar (almost identical) for all the models except the model with a strong field perpendicular to the collision, even the purely hydrodynamic simulation. This again reflects that the morpohology, and sink distributions of all these models are very similar. Again the model with a strong field perpendicular to the collision (BSYHighturb) shows a delay in the star formation rate, but then the star formation rate increases to values as high, or even higher than the other models. As previously, for a model with a strong field which lies neither parallel or perpendicular to the collision, the star formation rate lies between the BSY models and the other MHD and hydrodynamical models.

Magnetic field density relation
In Figure 7 we show the magnetic density relation for the higher virial parameter (top) and lower virial parameter models (lower) at the same times as shown in Figures 2  and 5. The figures show a region between densities of 10 −21 and 10 −17 g cm −3 where the magnetic field scales with the density with a relation slightly shallower than B ∝ ρ 1/2 . The initial cloud densities are ∼ 2 × 10 −22 g cm −3 , with the background density around 100 times lower. Below densities of 10 −21 g cm −3 , the magnetic field is roughly constant with density, although the behaviour is more noisy at low densities. Above densities of 10 −17 g cm −3 , there are relatively few particles to reliably infer a relation. The behaviour of the magnetic field with density is consistent with both Mocz et al. (2017) and Wurster et al. (2019), even though they model different environments of a turbulent molecular cloud, and disc formation around protostars in a smaller scale region of a molecular cloud. They find a B ∝ ρ 1/2 correlation at higher densities, whilst the magnetic field is relatively independent at lower densities. The transition however occurs at lower densities in our simulations, where the gas exhibits a lower range of densities, compared to Wurster et al. (2019). Similar to Wurster et al. (2019) and Mocz et al. (2017), we see that the magnetic field density correlation is largely independent of our initial conditions, where we are varying initial magnetic field strength, and the initial level of turbulence. We do see a slight tendency for the relation to extend to lower densities in the case with the weakest field parallel to the collision, which perhaps suggests in environments with relatively weaker fields and lower densities the relation extends to lower densities, in agreement with seeing an offset   Unlike Wurster et al. (2019), where the magnetic field strengths converge above densities of 10 −18 g cm −3 , we do see offsets in the magnetic field strength for the different models, although the models with a strong field parallel (BSX) and weak field perpendicular (BWY) to the shock are quite similar. The difference is perhaps not surprising, since the evolution is dominated by a strong shock, and the simulations are far from reaching any equilibrium in terms of the gas density distribution, dynamics, and magnetic field. The field strength is highest when the initial magnetic field is strong, and where the field is perpendicular to the collision and experiences amplification. The models which satisfy both, or neither of these properties, represent the outliers for the range of field strengths we use, and across the full range of magnetic field orientations.

Magnetic field and filaments
In this section we consider further the evolution of the magnetic field, the impact of the magnetic field on the evolution of the collision, and the relation of the magnetic field to the filaments formed at the shock interface. The filaments formed in our simulations typically have lengths in the range 1-10 pc, and large aspect ratios (> 10). In Figure 8, we show close ups of the shock region for the models with a strong field initially parallel (left, BSY), and perpendicular (right, BSX) to the direction of collision. Both panels are from the models with the higher virial parameter and higher level of turbulence. At this time (1.1 Myr), the clouds have collided, but there is relatively little fragmentation of the filaments formed from the shock at this point. In the case where the field is perpendicular to the collision, but parallel to the shock, the shock induces an increase in the strength of the field in this direction, as expected for a fast shock (Ryu & Jones 1995;Fukui et al. 2020). Theoretically, we would expect the magnetic field component parallel to the shock to increase by the same amount as the density. We see that both the density and magnetic field strength increase by a factor of several 10's. The large increase in magnetic pressure leads to a broader central filament or shocked region, rather than the narrow dense filaments seen in the other models. Unlike the other models, no sink particles have formed at this point, with the magnetic pressure preventing gravitational collapse (though sink particles do form later). In the left hand panel, the field strengths for the model with a strong magnetic field initially parallel to the collision, are much weaker. As expected from theory, there would be no increase in the field if it is perpendicular to the shock (parallel to the collision). As seen by comparing the panels in Figure 8, the field is more perpendicular to the filament in the left hand panel, but parallel to the filament in the right hand panel. For the left hand panel, there is likely some component of the field parallel to the shock, simply because the turbulent velocity field means the filaments formed from the shock are not completely perpendicular. As such this component will experience some amplification, and in some places the field acquires some component along the direction of the filament. For the left hand panel, the central filament has undergone some fragmentation and a few sink particles have formed along the central filament.
A similar phenomenon whereby the field is parallel to the weaker filaments, and perpendicular to the denser filaments, was noted in Wurster et al. 2019. Here we look at this a little more quantitatively. We do this simply by selecting the gas in the main central filaments in Figure 8, and comparing the magnetic field in the y direction with the magnitude of the magnetic field. For the BSX model (left), the density in the filaments increases to 10 −18 g cm −3 , the y component of the field is on average ∼ 10 −5 G, whilst the magnitude of the field is ∼ 3 × 10 −5 G. For the BSY model (right), the density in the filaments increases to 10 −20 g

Field parallel
Field perpendicular to filament to filament Low Strong field Not an outcome density High Weak field Weak field, or strong field density (field unimportant) with no parallel component Table 2. Possible outcomes for our simulations are shown according to the column and row labels, and the implication for the field strength.
After the time shown in Figure 8, the field becomes more aligned with the filaments in the BSX model, but it is not that long before there is more widespread fragmentation, the filaments seen in Figure 8 are less clear, and the field generally becomes more disordered. Similarly for the BSY model, the field becomes more disordered as many sink particles start to form. For the weaker field models, the field direction is similar to the original orientation of the field.
In Table 2 we have summarised possible outcomes, in terms of magnetic field orientation relative to filaments, and the relative density of the filament, and what initial setup or conditions correspond to this outcome. The Table assumes that the filaments are formed as the result of a high velocity collision, so if the filament formed by an alternative mechanism, it is possible that other inferences could be made. However we see that for our models, a low density filament with magnetic field parallel to the filament occurs if there is a strong magnetic field. The field cannot be perpendicular because that would lead to a high density filament. If the filament is high density, and the field parallel, then the field must be weak, because in this orientation, the field will be amplified, which will lead to a lower density filament if the field is initially strong. For a high density filament with magnetic field perpendicular to the filament, the most likely case is that the field is initially weak, so has little effect on the formation of the filament. Alternatively if the field is strong, the field must be strongly perpendicular to the filament. Interestingly we cannot obtain a relation between magnetic field orientation and filament density, because there is no one-to-one mapping between these two parameters, as in the high density case there are multiple initial conditions which produce a field perpendicular to the filament.

Properties of clusters
In this section we study the properties of the clusters formed in the different models. We use the DBSCAN program (Ester et al. 1996) to find the 3D distribution of sink particles, as described in . DBSCAN is a clustering technique which groups together points with similar densities. We use the same maximum separation distance as , 0.5 pc. We list the properties of the most massive cluster found in each simulation in Table 3. These properties are listed at a time of 2.3 Myr, which is the same time as shown in Figures 2 and 5. The properties of the clusters for the clouds with different initial magnetic field configurations are easiest to describe for the lower virial parameter cloud cases (lower 5 entries). We see that the properties of the clusters are very similar for all the models, except that with a strong field perpendicular to the collision (BSYHighturb), where a 5 times smaller cluster has formed. This is not surprising since by eye (Figure 2), the distribution of sink particles appears very similar in all models except the BSYHighturb model, where the distribution is completely different. By eye, there appear fewer sink particles in the main cluster in BSYHighturb, which although we need to take into account their masses, would suggest a lower mass, smaller radius cluster. In Figure 9 we plot the distribution of sink particles for the BSXLowturb, BWYLowturb, and BSYLowturb models, and again it is clear that the BSXLowturb and BWYLowturb models have very similar distributions of sink particles, and likewise the HydroLowturb and BWXLowturb look very similar to these two panels although not shown.
For the models with the higher virial parameter clouds (top 5 rows in Table 3), there is more variation in the cluster properties. Again the model with the strong field perpendicular to the collision, BSYHighturb, forms the smallest cluster, which is again unsurprising given the differences in morphology between this and other models, and the lower star formation rates for most of the duration of the simulations. Again, the distribution of sink particles in this model (Figure 9, right panel) is very different from the other models (top left and top centre panels). There is surprising variation in the cluster properties for the other runs. Figure 9 shows the sink particle distribution for the BSX-Highturb and BWYHighturb models. Here we see that the main accumulation of sink particles towards the top of the panel is identified as a single cluster in the BSXHighturb model, but only a subsection of this region is identified in the BWYHighturb model, hence a less massive cluster is picked out.
The likely difference between the two sets of models is that for the lower turbulence clouds, the sink particles tend to be located close together for all the models, and the DBSCAN algorithm readily finds a similar mass and size cluster in each case. For the higher turbulence clouds, the higher velocities leads to sink particles, and groups of sink particles which are slightly more disparate to each other. As we see for the BWYHighturb model ( Figure 9 top middle panel), the large distribution at the top of the panel is separated into about 3 different groups, one of which is picked out as the most massive cluster. In the BSXHighturb (Figure 9, top left panel), and likewise BWXHighturb models, these groups are selected as a single massive cluster. So for the higher turbulence cases, there is a similar distribution of sink particles, but the substructuring within them is slightly different. For the models with the magnetic field parallel to the collision, the substructuring is less pronounced, whereas for the BWYHighturb model with the field perpendicular to the collision, the substructuring is more evident. However it is difficult to say whether the difference is due to the magnetic field, as for the purely hydrodynamical model, the DBSCAN algorithm also only picks out a smaller subcluster (though larger than the BWYHighturb model), or in part due to the random variations in the distributions of the sink particles in different models.
Overall we see that the addition of a magnetic field tends to make a large difference to the cluster properties if the field is stronger and perpendicular to the collision. Otherwise, if the field is parallel, or weaker, the distribu- tion of sink particles is relatively unchanged on large scales, although we do see differences on smaller scales which may manifest in producing different substructure. For our models with bound clouds (lower virial parameter), we still readily produce massive clusters (> 10 4 M ) in all the models, which is probably not surprising given that they are strongly bound clouds. We also see that these clusters form on a timescale of ∼ 1.3 Myr. For the higher virial parameter clouds (except for the strong field parallel to the collision, BSYHighturb), we still see large congregations of sink particles with ∼ 10 4 M , though they are not always detected as a single cluster by our algorithm. For the strong field parallel to the collision cases, we see less massive clusters, but we note that the star formation is delayed in these instances. If we instead compare at the same time since star formation commences, for the higher turbulence model, the mass of the cluster is 1.4 ×10 4 M , so more comparable to the other models. However the morphology of the BSY model is still very different, and the resulting most massive cluster is denser and has a smaller half mass radius compared to the other models.

DISCUSSION -COMPARISON WITH PREVIOUS WORK
In this section we compare the results of our work to previous simulations. Previous work in this area including magnetic fields has tended to use grid based methods. Our finding of a large dependence on the initial direction of the magnetic field is in agreement with a number of previous works (Heitsch et al. 2009;Inoue & Inutsuka 2009;Fogerty et al. 2016). The earlier works (Heitsch et al. 2009;Inoue & Inutsuka 2009) did not include self gravity and focused on molecular cloud formation, but showed that gas densities only reached values comparable to molecular clouds when the field was parallel to the direction of collision. Similar to our work, the morphology of the gas resembles the hydro case when the field is parallel to the collision, but very different when the field is perpendicular. Körtgen & Banerjee (2015) and Fogerty et al. (2016) also see a delay in star formation with an inclined field compared to a field parallel to the collision, similar to the delay we see. Unlike our work, previous simulations see clearer differences with stronger and weaker fields when the field is parallel to the collision (Körtgen & Baner-jee 2015; Sakre et al. 2020;Heitsch et al. 2009), although we do not probe such high magnetic field strengths, where such differences may become more apparent. We also see that in our model with a strong field perpendicular to the collision, fragmentation is strongly compressed, in agreement with Wu et al. (2020).
In terms of the mass and time of sinks that form, previous simulations find different results. Inoue et al. (2018) find that the additional magnetic pressure leads to massive stars forming, whereas Fogerty et al. (2016) and Körtgen & Banerjee (2015) find the opposite. Zamora-Avilés et al.
(2018) find stars form earlier but Körtgen & Banerjee (2015) and Fogerty et al. (2016) see a delay. We clearly see a delay in agreement with the latter works. Our sink particles represent clusters rather than individual stars. At equivalent times, we see lower mass clusters in the runs with a strong perpendicular field. However if we take the time since the first sinks formed, the situation is less clear, and there is some indication that in the perpendicular case, denser if not necessarily more massive clusters can form.
Our simulations naturally produce filaments where the shock occurs from the colliding clouds. We see a tendency for magnetic fields to be more aligned with filaments when the magnetic field impedes the formation of dense gas and stars. The field is instead perpendicular to the filaments formed when the magnetic field has little effect. The densities of the filaments tend to be lower in the first case, and higher in the second case, in agreement with observations. This finding is comparable to Soler et al. (2013), who determined the alignment of the magnetic field in different density filaments formed from shocks in turbulent gas, and related the alignment of the field to the divergence of the velocity field (Soler & Hennebelle 2017). In our simulations we have a much simpler setup whereby we are modelling individual colliding streams of gas. However we may expect that the outcome for each filament formed in turbulence will have a similar dependence on the initial field strength and angle in our simulations. This idea was approached more rigorously by Inoue & Inutsuka (2009) and Inoue & Inutsuka (2016) who analytically relate the resultant density of the shock look to the initial angle and strength of the magnetic field prior to a shock. We at least see qualitatively similar behaviour in our models, although we note that there is not necessarily a one to one relationship between filament density and magnetic field properties.

CONCLUSIONS
We have performed SPMHD simulations of colliding clouds with magnetic fields to investigate the formation of massive clusters. We apply a timestep criterion which prevents large divergence errors. Although this criteria is not usually necessary in SPMHD calculations, we found it was required in the more extreme conditions of colliding clouds. Our simulations show that magnetic fields do not inhibit the formation of massive clusters from cloud cloud collisions, and YMCs can still form on timescales of 1-2 Myr, and the conclusions from our previous work  hold. Even in the case where the impact of the magnetic field is strongest, whilst we see a delay in star for- Figure 9. The sink particles are shown for the high virial parameter models (top) and low virial paramters (lower). The red points show the sink particles which are identified as the most massive cluster, as picked out by the DBSCAN algorithm. mation, we then see star formation rates which are similar to the other models.
Similar to previous work, we find that the initial orientation of the field has a strong effect on the outcome of the collision, and in our case the resultant clusters which form. As shown in Inoue & Inutsuka (2009), this can be related to the conditions of the magnetic field across a shock. If the field is initially parallel to the collision, it has little effect on the evolution even in our stronger field case. Thus dense filaments can form in the gas, and the magnetic field is aligned perpendicular to the filament, as seen in observations. If the field is initially perpendicular to the collision (parallel to the shock), then it is significantly amplified by the shock. As such the magnetic pressure prevents dense filaments forming and delays star formation, and leads to comparably lower density filaments. The magnetic field is aligned along the filament, again in agreement with observations. Thus the formation of the filaments determines the orientation of the field with respect to the filament. Despite the differences with orientation and field strength, we still see a B ∝ ρ 1/2 relation across our models, though the relations are offset from each other.
The influence of the magnetic field on the filaments leads to a corresponding impact on the clusters which form. In the cases where the field has little effect, namely when the field is parallel to the collision, or perpendicular but low strength, the star formation rates and clusters which form have similar properties to the purely hydrodynamic case. In our model where the field has a strong effect (perpendicular, higher strength), the formation of sink particles is delayed, and the resulting appearance and properties of the clusters which form are quite different. At the same absolute time, the clusters are considerably smaller compared to the other models. At the same duration since star formation commences, the clusters are still less massive though not by so much, but also quite dense.
Our simulations do not include stellar feedback, which we leave to future work. We would not expect feedback to have a large effect on the gas over the short timescales which our clusters form (e.g. Howard et al. 2018), though ionisation may start operating at relatively early times and may also have an impact on the magnetic field (Troland et al. 2016). We also have used fairly modest magnetic field strengths. We would expect the trends we find to continue to higher magnetic field strengths, although we don't necessarily relate our simulations to more extreme environments such as the Galactic Centre, where the dynamics and interstellar radiation field are also very different.
Finally, we are not aware of simulations similar to those presented here which have been carried out with SPH, rather than a grid code. It is encouraging that our simulations pro-duce results, and even filament morphologies, which are in agreement with previous grid code results.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.