There is no magnetic braking catastrophe: Low-mass star cluster and protostellar disc formation with non-ideal magnetohydrodynamics

TL;DR

This study uses advanced non-ideal MHD simulations to show that magnetic fields do not prevent protostellar disc formation or cause a magnetic braking catastrophe in low-mass star cluster formation.

Contribution

First radiation non-ideal MHD simulations resolving fragmentation down to the opacity limit, showing magnetic fields do not hinder disc formation or cause magnetic braking catastrophe.

Findings

Magnetic field strength in dense cores is independent of large-scale environment.

Star formation rate decreases with increasing initial magnetic field strength.

Protostellar discs form regardless of magnetic field strength, indicating turbulence influences disc formation.

Abstract

We present results from the first radiation non-ideal magnetohydrodynamics (MHD) simulations of low-mass star cluster formation that resolve the fragmentation process down to the opacity limit. We model 50~M$_\odot$ turbulent clouds initially threaded by a uniform magnetic field with strengths of 3, 5 10 and 20 times the critical mass-to-magnetic flux ratio, and at each strength, we model both an ideal and non-ideal (including Ohmic resistivity, ambipolar diffusion and the Hall effect) MHD cloud. Turbulence and magnetic fields shape the large-scale structure of the cloud, and similar structures form regardless of whether ideal or non-ideal MHD is employed. At high densities ($10^6 \lesssim n_{\rm H} \lesssim 10^{11}$~cm$^{-3}$), all models have a similar magnetic field strength versus density relation, suggesting that the field strength in dense cores is independent of the large-scale environment. Albeit with limited statistics, we find no evidence for the dependence of the initial mass function on the initial magnetic field strength, however, the star formation rate decreases for models with increasing initial field strengths; the exception is the strongest field case where collapse occurs primarily along field lines. Protostellar discs with radii $\gtrsim 20$~au form in all models, suggesting that disc formation is dependent on the gas turbulence rather than on magnetic field strength. We find no evidence for the magnetic braking catastrophe, and find that magnetic fields do not hinder the formation of protostellar discs.

Figures

Tables

Table 1. Table 1: Summary of the stellar population properties at the common time of t = 1.45 𝑡 1.45 t=1.45  t ff ff {}_{\text{ff}}  (denoted by superscript c), and at the final time of each simulation (the final times are listed in the t final subscript 𝑡 final t_{\text{final}} column, and denoted by superscript f). Entries are not duplicated for the models with t final = subscript 𝑡 final absent t_{\text{final}}= 1.45 1.45 1.45  t ff ff {}_{\text{ff}} . M stars subscript 𝑀 stars M_{\text{stars}} is the total mass in stars, N stars subscript 𝑁 stars N_{\text{stars}} is the total number of stars, N mergers subscript 𝑁 mergers N_{\text{mergers}} is the number of mergers prior to the common or final time, N systems subscript 𝑁 systems N_{\text{systems}} is the number of stellar systems with 1, 2, 3 or 4 stars, respectively, and N discs subscript 𝑁 discs N_{\text{discs}} is the number of discs around single, binary, triple, and quadruple stars; the number in brackets includes the total number of discs of each classification that are also part of higher order stellar systems.

Name	$M_{\text{stars}}^{\text{c}}$	$N_{\text{stars}}^{\text{c}}$	$N_{\text{mergers}}^{\text{c}}$	$N_{\text{systems}}^{\text{c}}$	$N_{\text{discs}}^{\text{c}}$	$t_{\text{final}}$	$M_{\text{stars}}^{\text{f}}$	$N_{\text{stars}}^{\text{f}}$	$N_{\text{mergers}}^{\text{f}}$	$N_{\text{systems}}^{\text{f}}$	$N_{\text{discs}}^{\text{f}}$
	[M⊙]			1, 2, 3, 4	1, 2, 3, 4	[t${}_{\text{ff}}$]	[M⊙]			1, 2, 3, 4	1, 2, 3, 4
N03	7.51	10	0	3, 0, 1, 1	2 (5), 0 (3), 1, 1	1.46	7.72	10	0	3, 0, 1, 1	2 (7), 0 (3), 1, 1
N05	3.83	17	0	4, 1, 1, 2	4 (13), 1 (4), 1, 2	1.53	7.57	27	0	8, 5, 3, 0	3 (17), 5 (8), 3, 0
N10	4.48	8	0	2, 1, 0, 1	2 (7), 1 (2), 0, 1	1.61	9.29	10	1	3, 0, 1, 1	2 (9), 0 (2), 1, 1
N20	6.51	10	1	2, 1, 2, 0	0 (8), 1 (3), 2, 0	1.59	9.99	11	1	2, 3, 1, 0	0 (7), 3 (4), 1, 0
Hyd	7.93	19	0	6, 0, 3, 1	4 (14), 0 (4), 3, 1
I03	7.10	11	1	2, 1, 1, 1	1 (3), 0 (3), 1, 1
I05	3.34	11	1	5, 3, 0, 0	4 (10), 3 (3), 0, 0	1.55	7.64	19	2	10, 1, 1, 1	4 (11), 1 (3), 1, 1
I10	4.14	7	0	3, 0, 0, 1	3 (6), 0 (1), 0, 1	1.52	6.37	9	0	3, 1, 0, 1	3 (5), 1 (2), 0, 1
I20	6.00	11	2	4, 0, 1, 1	1 (6), 0 (2), 1, 1

Table 2. Table 2: Summary of the stellar population properties for the models with μ 0 = 5 subscript 𝜇 0 5 \mu_{0}={5} . The models are compared at t ≈ 1.36 𝑡 1.36 t\approx 1.36  t ff ff {}_{\text{ff}}  (for all models, superscript c) and at t ≈ 1.45 𝑡 1.45 t\approx 1.45  t ff ff {}_{\text{ff}}  (for the ideal models, superscript f). The columns are the same as in Table  1 .

Name	$M_{\text{stars}}^{\text{c}}$	$N_{\text{stars}}^{\text{c}}$	$N_{\text{mergers}}^{\text{c}}$	$N_{\text{systems}}^{\text{c}}$	$N_{\text{discs}}^{\text{c}}$	$M_{\text{stars}}^{\text{f}}$	$N_{\text{stars}}^{\text{f}}$	$N_{\text{mergers}}^{\text{f}}$	$N_{\text{systems}}^{\text{f}}$	$N_{\text{discs}}^{\text{f}}$
	[M⊙]			1, 2, 3, 4	1, 2, 3, 4	[M⊙]			1, 2, 3, 4	1, 2, 3, 4
N05h	0.78	4	0	4, 0, 0, 0	4 (4), 0 (0), 0, 0
N05	1.43	7	0	0, 0, 1, 1	0 (7), 0 (2), 1, 0
I05h	0.61	3	0	3, 0, 0, 0	3 (3), 0 (0), 0, 0	2.13	14	0	5, 1, 1, 1	5 (14), 0 (2), 1, 1
I05	1.28	5	0	0, 1, 1, 0	0 (5), 1 (1), 0, 0	3.25	12	0	5, 2, 1, 0	4 (9), 1 (1), 1, 0

Equations

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathrm{t}}$$

$$\rho \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\mathbf{B}}{\rho }\right)$$

$$\rho \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{E}{\rho }\right)$$

$$+$$

$$\rho \frac{\mathrm{d}u}{\mathrm{d}t}$$

$$\rho {\Lambda }_{\text{ISR}}$$

$${\nabla }^2\Phi$$

$$\text{d}{t}_{\text{non-ideal}}=C_{\text{non-ideal}}\frac{h^2}{\left|\eta \right|},$$

$${\mu }_0=\left(\frac{M}{\pi R_0^2{B}_0}\right){\left(\frac{M}{{\Phi }_{\text{B}}}\right)}_{\text{crit}}^{-1},$$

$${\left(\frac{M}{{\Phi }_{\text{B}}}\right)}_{\text{crit}}=\frac{c_1}{3}\sqrt{\frac{5}{G}},$$

$$B_{\text{mean}}=10^{\left({\Sigma }_i^N\log B_i\right)/N},$$

$$B_{\text{mean}}^j=\frac{N10^{\left[{\sum }_{(B_i^j>0)}^N\log (B_i^j)\right]/N}-n10^{\left[{\sum }_{(B_i^j<0)}^n\log (-B_i^j)\right]/n}}{N+n}$$

Full text

††pagerange: There is no magnetic braking catastrophe: Low-mass star cluster and protostellar disc formation with non-ideal magnetohydrodynamics–19††pubyear: 2019

1 Introduction

Stars are born and evolve in a dynamic environment, and are neither born nor evolve in isolation. Therefore, how a star forms and evolves is directly influenced by its environment. Environmental influences in star forming molecular clouds include turbulence (e.g. Padoan & Nordlund, 2002; McKee & Ostriker, 2007), radiative feedback (e.g. Hatchell et al., 2013; Sicilia-Aguilar et al., 2013; Fall et al., 2010; Skinner & Ostriker, 2015; Kim et al., 2016), magnetic fields (e.g. Heiles & Crutcher, 2005) and low ionisation fractions (e.g. Mestel & Spitzer, 1956; Nakano & Umebayashi, 1986; Umebayashi & Nakano, 1990).

Molecular clouds are large, structured complexes, with typical masses of $10^{3}-10^{5}$ M⊙, sizes of $2-15$ pc (e.g. Bergin & Tafalla, 2007, and references therein) and magnetic field strengths of a few times the critical mass-to-flux ratio (e.g. Crutcher, 2012, and references therein). The clouds often contain a few, nearly parallel, large filaments, whose length may be comparable to the length of the cloud (e.g. Tachihara et al., 2002; Burkert & Hartmann, 2004). Magnetic fields are observed to be approximately perpendicular to dense filaments (Heyer et al., 1987; Alves et al., 2008; Franco et al., 2010), while in the less dense regions, magnetic fields are approximately parallel to the structure (i.e. to striations; e.g. Goldsmith et al., 2008; Planck Collaboration et al., 2016a, b). Dense filaments, and especially the places where they intersect, are home to regions of intense star formation. Individual young stars and their protostellar discs have been resolved in these regions (e.g. Dunham et al., 2011; Lindberg et al., 2014; Tobin et al., 2015; Gerin et al., 2017), and in many cases, the magnetic field geometry has been inferred in and near the discs (e.g. Stephens et al., 2014, 2017; Cox et al., 2018).

There have been numerous numerical studies of the formation and evolution of large ($M\geq 10^{3}$ M⊙) molecular clouds, that both exclude (e.g. Walch et al., 2012; Dale et al., 2014; Dale, 2017; Lee & Hennebelle, 2018a, b) and include (e.g. Gendelev & Krumholz, 2012; Myers et al., 2014; Geen et al., 2015; Girichidis et al., 2016; Geen et al., 2018; Lee & Hennebelle, 2019) magnetic fields. These studies allow for a direct comparison with observations on global (whole cloud) scales. However, limited numerical resolution means that stellar properties must be inferred from sub-grid models since models on these scales currently do not resolve the fragmentation of individual stars.

Simulating the formation and evolution of individual stars requires resolving the ‘opacity limit for fragmentation’ — the minimum Jeans mass as a function of density where radiation is trapped by dust opacity, leading to gas heating and a stable fragment (Rees, 1976; Low & Lynden-Bell, 1976). This minimum mass is several Jupiter masses and, when using smoothed particle hydrodynamics, must be resolved by 50-100 particles (Bate & Burkert, 1997). At our current computational capabilities, this restricts us to clouds of 50-500 M⊙while modelling the formation of individual stars without sub-resolution models.

To date, many studies on this scale either ignored magnetic fields or neglected radiation by employing a barotropic equations of state (where temperature is simply prescribed as a function of density). Simulations of 50 M⊙ clouds by Bate et al. (2003) using a barotropic equation of state produced as many brown dwarfs as stars — hence over-producing the number of low-mass objects. Larger simulations of 500 M⊙ clouds (Bate, 2009a), or with better statistics from multiple realisations of 50 M⊙ clouds (Liptai et al., 2017), confirmed this discrepancy with the observed initial mass function (IMF). Once radiative feedback was included (Bate, 2009b, 2012; Offner et al., 2009), the proportion of brown dwarfs decreased, bringing the mass function into good agreement with the observed IMF. Other low-mass cluster simulations have investigated the effect of different cloud density and velocity structure (e.g. Bate, 2009c; Liptai et al., 2017; Jones & Bate, 2018a), metallicity (e.g. Bate, 2014, 2019) and stellar radiation algorithms (e.g. Jones & Bate, 2018b).

Here, we focus on the effect of magnetic fields on low mass clouds (as previously studied in, e.g., Price & Bate, 2008, 2009; Hennebelle et al., 2011; Myers et al., 2013); for reviews on the role of magnetic fields in molecular clouds, see Crutcher (2012) and Hennebelle & Inutsuka (2019). These studies concluded that magnetic fields influenced the global evolution of the cloud, with magnetic pressure-supported voids and filamentary structures that followed the magnetic field lines forming in the strong field models. The star formation rate was found to decrease with increasing initial magnetic field strength. Our initial conditions most closely match that of Price & Bate (2008, 2009), although they used Euler potentials to model the magnetic field, whereas we directly evolve the induction equation. Euler potentials enforce $\bm{\nabla}\cdot\bm{B}\equiv 0$, however, they cannot capture the full physics of magnetohydrodynamics (MHD) since helicity is identically zero (e.g. Price, 2010, 2012); this means that complex magnetic fields structures, such as winding or combined poloidal/toroidal fields, can not be modelled with Euler potentials, but they can be with the direct evolution of the magnetic field (and divergence cleaning) that we use in this paper.

In this paper, we present results from low-mass radiation non-ideal magnetohydrodynamic star cluster simulations that resolve the fragmentation process down to the opacity limit. In Section 2, we present our methods, and in Section 3 we present our initial conditions. We present our results in Section 4 and discuss them with respect to the literature in Section 5. Our conclusions are presented in Section 6. Supplementary discussion of the dependence of the calculations on resolution and of the disc populations that are produced in the calculations are presented in Appendices A and B, respectively.

2 Methods

2.1 Radiation non-ideal magnetohydrodynamics

We solve the equations of self-gravitating, radiation non-ideal magnetohydrodynamics in the form

[TABLE]

where ${\rm d}/{{\rm d}t}\equiv\partial/\partial t+\bm{v}\cdot\nabla$ is the Lagrangian derivative, $\rho$ is the density, ${\bm{v}}$ is the velocity, $p$ is the gas pressure, ${\bm{B}}$ is the magnetic field, $\bm{J}=\bm{\nabla}\times\bm{B}$ is the current density, $\Phi$ is the gravitational potential, $E$ is the radiation energy density, $F$ is the radiative flux, P is the radiation pressure tensor, $c$ is the speed of light, $G$ is the gravitational constant, and $\mathbb{I}$ is the identity matrix. The non-ideal MHD coefficients for Ohmic resistivity, the Hall effect and ambipolar diffusion are given by $\eta_{\text{OR}}$, $\eta_{\text{HE}}$ and $\eta_{\text{AD}}$, respectively. We assume units for the magnetic field such that the Alfvén speed is $v_{\rm A}=|B|/\sqrt{\rho}$ (see Price & Monaghan, 2004a).

For the treatment of the equation of state and radiative transfer, we use the combined radiative transfer and diffuse interstellar medium model of Bate & Keto (2015). The underlying method of radiative transfer is grey flux-limited diffusion, using the implicit algorithms devised by Whitehouse et al. (2005) and Whitehouse & Bate (2006). Bate & Keto extended this method to treat gas and dust temperatures separately, and included various heating and cooling mechanisms that are particularly important at low densities. In Eqns. 5 to 7, $B_{\text{P,d}}$ is the frequency-integrated Planck function at the dust temperature, and $B_{\text{P,g}}$ is the frequency-integrated Planck function at the gas temperature, where the Planck function depends on temperature as $B_{\text{P}}=(\sigma_{\text{B}}/\pi)T^{4}$ where $\sigma_{\text{B}}$ is the Stefan-Boltzmann constant. The quantities $\kappa_{\text{d}}$ and $\kappa_{\text{g}}$ and the frequency-integrated opacities of the dust and gas, respectively. The quantity $\Gamma$ includes cosmic ray, photoelectric, and molecular hydrogen formation heating terms. The quantity $\Lambda_{\text{line}}$ provides gas cooling due to atomic and molecular line emission, and $\Lambda_{\text{gd}}$ treats thermal energy transfer between the gas and the dust. Eqn. 7 sets the dust to be in thermal equilibrium with the total radiation field, taking into account the energy transfer between the gas and the dust. The total radiation field is comprised of two components: the external interstellar radiation field ($\Lambda_{\text{ISM}}$ is the dust heating rate per unit volume due to the interstellar radiation field), and the radiation energy density of the flux-limited diffusion radiation field.

The equation of state is described by Whitehouse & Bate (2006). Briefly, we employ an ideal gas equation of state that assumes a 3:1 mix of ortho- and para-hydrogen (see Boley et al., 2007) and treats the dissociation of molecular hydrogen and the ionisations of hydrogen and helium. The mean molecular weight is taken to be $\mu_{\rm g}=2.38$ at low temperatures, and we use dust opacity tables from Pollack et al. (1985) and gas opacity tables from Ferguson et al. (2005).

2.2 Smoothed particle radiation non-ideal magnetohydrodynamics

We use the three-dimensional smoothed particle hydrodynamics (SPH) code sphNG that originated from Benz (1990), but has since been substantially extended to include individual particle timesteps and sink particles (Bate et al., 1995), variable smoothing lengths (Price & Monaghan, 2004b, 2007), radiative transfer (Whitehouse et al., 2005; Whitehouse & Bate, 2006; Bate & Keto, 2015), magnetohydrodynamics (Price, 2012) and non-ideal MHD (Wurster et al., 2014, 2016). The code is parallelised using both OpenMP and MPI.

The density and smoothing length of each SPH particle are calculated iteratively by summing over its nearest neighbours and using $h=1.2(m/\rho)^{1/3}$ where $h$, $m$ and $\rho$ are the SPH particle’s smoothing length, mass and density, respectively (Price & Monaghan, 2004b, 2007). The remainder of the terms, including our magnetic variable of $\bm{B}/\rho$ (Eqn. 3), are solved using the standard smoothed particle magnetohydrodynamics (SPMHD) scheme (for a review, see Price, 2012). We use the artificial viscosity described in Price & Monaghan (2005) to capture shocks. For magnetic stability, we use the Børve et al. (2001) source-term approach and constrained hyperbolic divergence cleaning method of Tricco et al. (2016). Artificial resistivity is included using the Tricco & Price (2013) resistivity switch, which produces higher artificial resistivity (Wurster et al., 2017a) and hence weaker magnetic fields strengths than the algorithm from Price et al. (2018) that we have used in our recent studies of isolated star formation (Wurster et al., 2018a, b, c). Our SPMHD algorithm allow us to capture the complex magnetic field structures, unlike in previous studies that used the Euler potential method (Price & Bate, 2008, 2009). To evolve the simulation, we employ a second-order Runge-Kutta-Fehlberg integrator (Fehlberg, 1969).

We use explicit timestepping for the Hall effect and ambipolar diffusion, with

[TABLE]

where $C_{\text{non-ideal}}=1/2\pi$ is a constant equivalent to the Courant number (Wurster et al., 2014). Ohmic resistivity is calculated using an implicit solver (see Appendix A of Wurster et al., 2018a).

3 Initial conditions

We initialise our clusters as spheres of cold, dense gas embedded in a warm medium. The sphere has an initial radius of $R_{0}=0.1875$ pc and mass of 50 M⊙, yielding an initial uniform density of $1.22\times 10^{-19}$ g cm*-3*; the characteristic timescale is the free-fall time of $t_{\text{ff}}=1.9\times 10^{5}$ yr. The initial temperature structure of the sphere is set such that it is in thermal equilibrium with the interstellar radiation field; the initial temperature increases from $T\sim 8.8$ K at the centre of the sphere to $\sim 13$ K at its edge. The sphere is embedded in a cube whose side-length is 0.75 pc and density is 30 times lower, with the warm and cold media in pressure equilibrium. We use boundary conditions such that the hydrodynamic forces are periodic but gravitational forces are not. This ‘sphere-in-box’ method prevents the need for magnetic boundary conditions at the edge of the initial sphere, which would immediately become deformed.

We impose an initial supersonic ‘turbulent’ velocity field as described by Ostriker et al. (2001) and Bate et al. (2003). We generate a divergence-free random Gaussian velocity field with a power spectrum $P(k)\propto k^{-4}$, where $k$ is the wavenumber. This is generated on a $128^{3}$ uniform grid, and the velocities of the particles are then interpolated from this grid. The velocity field is normalised so that the kinetic energy of the turbulence is equal to the magnitude of the gravitational potential energy of the initial sphere (i.e. excluding the warm medium). This equates to an initial rms Mach number of the turbulence of $\mathcal{M}=6.4$. Each model uses the same initial velocity field, which uses the same random seed as in many of our previous simulations (e.g. Bate et al., 2003; Price & Bate, 2008; Bate, 2009b, 2014; Liptai et al., 2017). We caution that a different random seed may yield different results (e.g. Liptai et al., 2017; Geen et al., 2018). For example, in their study, Geen et al. (2018) found the star formation efficiency to vary between 6 and 23 per cent simply by changing the random seed. Ideally, we would perform many realisations of each model with different seeds, but this is beyond our current computing resources.

We thread the entire domain with a uniform magnetic field initially parallel to the $z$-axis. The non-ideal MHD coefficients of Ohmic resistivity, ambipolar diffusion and the Hall effect are computed using version 1.2.1 of the Nicil library (Wurster, 2016). This includes a heavy ion represented by Magnesium (Asplund et al., 2009), a light ion represented by the hydrogen and helium compounds, and a single grain species, $n_{\text{g}}$. We assume a single grain species with radius $a_{\text{g}}=0.1~{}\mu$m and bulk density of $\rho_{\text{b}}=3$ g cm*-3* (Pollack et al., 1994), respectively, that includes three populations with charges $Z=-1,0,+1$, respectively, where $n_{\text{g}}=n_{\text{g}}^{-}+n_{\text{g}}^{0}+n_{\text{g}}^{+}$ to conserve grain density. The grain number density is calculated from the local gas density, assuming a dust-to-gas ratio of 0.01.

To model the long-term evolution of the cluster, we include sink particles (Bate et al., 1995) with accretion radius $r_{\text{acc}}=0.5$ au. When a gas particle reaches $\rho=10^{-5}$ g cm*-3*, near the end of the second collapse phase, it is replaced with a sink particle; all particles within $0.25$ au of the centre of the sink particle are automatically accreted, and all particles within $0.5$ au are tested to determine if they are to be accreted. Since we resolve the opacity limit, each sink particle represents one star or brown dwarf, and the properties of the sink reflect the properties of the star. We merge sink particles that approach within 27 R⊙ of each other.

In this study, we model a 50 M⊙ cloud in order to resolve the minimum local Jeans mass throughout the simulation, which is $\approx$0.0011 M⊙ during the isothermal collapse phase (Bate & Burkert, 1997; Bate et al., 2003, and discussion therein). Resolving this limit requires $\sim$75 particles per Jeans mass, or a minimum of $3.5\times 10^{6}$ particles for a 50 M⊙ cloud, as used in our previous studies. By default, we use $5\times 10^{6}$ equal mass SPH particles in the cold dense sphere, and an additional $2.3\times 10^{6}$ particles in the warm medium. To test the effect of resolution, we also perform two simulations using $35\times 10^{6}$ equal mass SPH particles in the sphere; see Appendix A.

We test four magnetic field strengths: $B_{0}=6.48,3.89,1.94\text{ and }0.972\times 10^{-5}$ G, which correspond to normalised mass-to-flux ratios of $\mu_{0}=3,5,10$ and 20. The normalised mass-to-flux ratio is given by

[TABLE]

where

[TABLE]

is the critical mass-to-flux ratio in CGS units; here, $\Phi_{\text{B}}$ is the magnetic flux threading the surface of the (spherical) cloud, and $c_{1}\simeq 0.53$ is a parameter determined numerically by Mouschovias & Spitzer (1976).

For each magnetic field strength at our fiducial resolution ($5\times 10^{6}$ particles), we model both an ideal and non-ideal MHD cluster. We also model a purely hydrodynamical model, for a total of nine simulations. The magnetised models are labelled Ixx and Nxx, where xx is the normalised mass-to-flux ratio. We refer to the hydrodynamical model as ‘hydro’ or Hyd.

4 Results

The different dynamics and physical processes included in each model led to variations in computational expense, so not all models had reached the same final time. Our primary analysis is performed at $t=1.45$ t${}_{\text{ff}}$, the minimum time reached by all simulations.

4.1 Large scale structure

4.1.1 Gas density evolution

The top five rows of Fig. 1 show the time evolution of the non-ideal MHD and hydro models (left to right, respectively). The bottom row shows the ideal MHD models at $t=1.45$ t${}_{\text{ff}}$, which may be compared to the row above showing the corresponding non-ideal models at the same time.

Magnetic fields have an immediate effect on the evolution of the large scale cloud structure, with visible differences in the large scale structure appearing as early as $t=0.15$ t${}_{\text{ff}}$. The weakly magnetised models become more filamentary at an earlier time, yielding slightly higher maximum densities. Stronger magnetic fields, however, can be seen to smooth out the large-scale shock structure, as also found by Price & Bate (2008); the effect is similar to increasing the gas pressure in hydrodynamical simulations (Bate & Bonnell, 2005). By $t=0.40$ t${}_{\text{ff}}$, the evolution of the hydro model visibly differs from the magnetised models with $\mu_{0}={20}$, showing that even a weak magnetic field can influence the large scale dynamics.

During the early evolution ($t\lesssim 0.40$ t${}_{\text{ff}}$), there are many transient features in the gas, with filaments forming, colliding and dispersing. Persistent, dense clumps begin to form at $t=0.40$ t${}_{\text{ff}}$ (top row) in the hydro model first. Due to small evolutionary changes already caused by the magnetic field, there is no obvious relation between the formation time of the first clumps and the initial magnetic field strength. As the evolution continues, the clumps become denser and filamentary, which is more pronounced in the more strongly magnetised models. Each model contains two primary groups of clumps/filaments, and these ultimately merge into a single, primary filament in the models with $\mu_{0}={3}$ and $5$.

At $t=0.90$ t${}_{\text{ff}}$ (third row), there is a weak qualitative trend regarding the dense gas near $x\sim 0.05$ pc and $z\sim 0$: As the magnetic field strength is increased, this clump becomes more compact and denser. However, this trend is broken by N03, in which this clump does not exist. Rather, the denser clump occurs at $x\sim-0.05$ pc and $z\sim 0$ and is more elongated. Although this qualitative difference may appear trivial, it demonstrates the effect of magnetic fields on the large-scale cloud structure. This qualitative difference also suggests that the properties of N03 may not conform to the trend of the remaining three non-ideal models plus the hydro model. This deviation from the trend will permeate most of our subsequent analysis. This difference is a result of the clump preferentially forming along the magnetic field lines in N03, thus becoming filamentary early. In N05, there are multiple regions where the gas has collapsed vertically into small clumps, and then the turbulent motion ultimately joins these clumps into a filament. In the remaining models, filaments typically form without any preference to the $z$-axis (i.e., to the initial magnetic field direction).

As the clusters evolve to $t=1.15$ t${}_{\text{ff}}$ (fourth row), the dense regions in N05 to Hyd collapse and host the birth of the first stars; thus the first stars form in these models in similar locations. The dense clump in N03 also collapses to form stars, but this dense structure is already more elongated and denser than in the remaining models.

As the evolution continues, the filament becomes denser in N03, while in the remaining models, the turbulent motion of the gas brings the various clumpy regions together to form the filament. Thus, at $t=1.45$ t${}_{\text{ff}}$, each model contains one large scale ($\sim$0.2 pc) filament, but their formation mechanisms differ between the strong and weak field cases.

This primary filament is more vertical in the $yz$-plane than the $xz$-plane shown in Fig. 1, but this is simply a result of the initial turbulent velocity field. In addition to the primary filament, there are many lower density filaments in all models, most of which emanate from the primary filament. The low-density filaments are denser but shorter in the more strongly magnetised models, since only dense objects can overcome the magnetic field strength in these models.

The evolutionary sequence described above occurs in both our ideal and non-ideal MHD models; non-ideal effects have little effect on the large scale evolution.

Fig. 2 shows the gas+stellar mass above density thresholds of $\rho>10^{-17}$ (top panel) and $10^{-15}$ g cm*-3* (lower panel) as a function of time for the hydro and non-ideal MHD models. The majority of the gas remains at low densities, with 76 to 87 per cent having $\rho\leq 10^{-17}$ g cm*-3* (recall that $\rho_{0}\approx 10^{-19}$ g cm*-3*). Prior to $t\approx{1.32}$ t${}_{\text{ff}}$, there is a trend of less gas with $\rho>10^{-17}$ g cm*-3* for models with stronger magnetic field strengths (top panel). This is reasonable since magnetic fields support gas from gravitationally collapsing, and was also found by Price & Bate (2008). The exception to this trend is N03, which has the same amount of gas above $10^{-17}$ g cm*-3* as the hydro model. This is a result of the rapid collapse of gas along the field lines to form the dense filament, as discussed above.

4.1.2 Magnetic field evolution

Fig. 3 shows the gas column density over-plotted with the magnetic field direction in both the $xz$- and $xy$-planes.

Given the initial magnetic field of $\bm{B}=-B_{0}\hat{\bm{z}}$, the mostly vertical magnetic field in the $xz$-plane is to be expected (top row). For stronger initial magnetic field strengths, the field remains closer to the initial conditions as the cluster evolves, especially on scales larger than shown in the figure. In N03, the magnetic field crosses the primary filament approximately perpendicularly, in agreement with Planck observations (e.g. Planck Collaboration et al., 2016a, b) and optical polarisation measurements in Taurus molecular cloud (Goldsmith et al., 2008; Chapman et al., 2011) and the Pipe Nebula (Alves et al., 2008), indicating that such magnetically controlled filament formation is likely to occur in the ISM; we note that the Planck observations (Planck Collaboration et al., 2016b) typically resolve down to 0.4 pc, whereas we can resolve the gas density in the filaments down to $\rho\sim 5\times 10^{-3}$ pc. In the models with weaker magnetic fields, the magnetic field also crosses the filaments perpendicularly, but there is also a considerable parallel component, especially in the regions surrounding the filament. This parallel component is more prominent in the less dense filament and in the weaker magnetic field models.

From observations, it is expected that the magnetic field is perpendicular to the dense filaments and parallel to the low-density filaments, and our results broadly agree with this, although our agreement is more tenuous in the initially weakly magnetised models. This tenuous agreement may also be a result of the initial preferential direction of the magnetic field. When we observe the cluster from the $xy$-plane (bottom row of Fig. 3), the magnetic field better follows the observationally expected paths, since in this plane, there is no preferential magnetic field direction. In N20, the filament is only dense enough in a few regions such that the magnetic field crosses it perpendicularly.

4.1.3 B vs $\rho$

Fig. 4 shows the mean magnetic field strength,

[TABLE]

along with the range containing 95 per cent of the magnetic field strengths, in the density range from $\rho=5\times 10^{-21}$ to $10^{-12}$ g cm*-3*; recall that the initial cloud density is $\rho_{0}=1.22\times 10^{-19}$ g cm*-3*. These curves are the average profiles of each snapshot just prior to each sink formation event before $1.45$ t${}_{\text{ff}}$. Although the calculations follow protostellar collapse to $\rho=10^{-5}$ g cm*-3* before sink particle insertion, the numerical resolution is much lower than in our isolated studies of star formation and the artificial resistivity is likely to weaken the field within the first hydrostatic cores ($\gtrsim 10^{-12}$ g cm*-3*), thus we do not comment on the magnetic field at higher densities.

At low densities ($\rho\lesssim 10^{-19}$ g cm*-3*), the magnetic field strength remains close to the initial field strength; there is larger scatter around this value for the high-$\mu_{0}$ models, and little scatter for low-$\mu_{0}$ models. At intermediate densities ($10^{-17}\lesssim\rho/($g cm*-3*$)\lesssim 10^{-13}$), the relationship converges for the eight models, with the average magnetic field strength increases approximately as $B\propto\rho^{1/2}$; in this range, the average magnetic field strengths only differ by a factor of a few between the models, which is less than the initial spread. All models also have approximately the same dispersion of field strengths in this range, which is approximately an order of magnitude.

The overall scaling of $B\propto\rho^{1/2}$ in this range, and the relative independence of the average magnetic strength on the initial field strength is consistent with the results of recent moving mesh simulations of magnetised, turbulent molecular clouds of Mocz et al. (2017). They find that $B\propto\rho^{2/3}$ for weakly magnetised clouds, transitioning to $B\propto\rho^{1/2}$ for more strongly magnetised clouds, with the transition occurring at around a Alfvénic Mach number of unity. Our initial conditions (magnetic field strength and density) are similar to their more highly magnetised cases. It is interesting to note that they also find that the typical magnetic field strength at densities $\rho\gtrsim 10^{-18}$ g cm*-3* is independent of their initial magnetic field strength (see their fig. 4). They perform simulations in which the initial magnetic field strength is varied by two orders of magnitude (compared to the range of a factor of $\approx 7$ in our simulations). They do not obtain as great a dispersion in magnetic field strengths at $\rho\gtrsim 10^{-15}$ g cm*-3*, but this may simply be because they stop their calculations when the first region collapses.

The wide range of magnetic field strengths in all the models prior to the hydrostatic core phase (i.e. in the range $10^{-17}\lesssim\rho/($g cm*-3*$)\lesssim 10^{-13}$) suggest a wide range of magnetic initial conditions for star formation. The maximum magnetic field strengths here are slightly higher than those presented in the isolated star formation study of Bate et al. (2014), who used similar mass-to-flux ratios (but slightly stronger absolute magnetic field strengths). Given the demonstration of the magnetic braking catastrophe (Allen et al., 2003) in that study, this suggests that disc formation may be hindered in the current simulations, or at least in the regions permeated by strong magnetic fields as given by the upper limits shown in Fig. 4). Discs and disc formation will be further discussed in Section 4.4.

There is no significant difference in results between the ideal and non-ideal MHD models (compare top to bottom panels of Fig. 4.) In agreement with our previous studies (e.g. Wurster et al., 2016, 2018a, 2018b), this indicates that non-ideal MHD processes only affect high-density ($\rho\gtrsim 10^{-12}$ g cm*-3*) gas.

4.1.4 Non-ideal vs ideal MHD

In the outer regions of the cloud where $\rho\lesssim 10^{-18}$ g cm*-3*, the ionisation fraction is $n_{\text{e}}/n\gtrsim 10^{-6}$. This relatively high ionisation rate is a result of cosmic rays easily being able to ionise the low-density gas; in this gas, the temperature is too cold, $T\sim\mathcal{O}\left(10\right)$ K, and the dust grain population is too small to contribute to the ionisation fraction. The drift velocity (i.e. $\bm{v}_{\text{drift}}=\bm{v}_{\text{neutral}}-\bm{v}_{\text{ion}}$) in the bulk of the cloud is \sim$$10^{-4}-10^{-1} km s*-1*, suggesting that the ions are well-coupled to the neutral particles. Given the high ionisation fractions, the non-ideal processes are too weak to contribute significantly to the evolution of the gas. Therefore, the evolution of the large-scale structure is essentially independent of whether ideal or non-ideal MHD is assumed.

Comparison of the bottom two rows of Fig. 1 confirms this. The quantity of gas with $\rho\gtrsim 10^{-15}$ g cm*-3* (c.f. Fig. 2; ideal models not shown) is also not strongly sensitive to non-ideal MHD effects. At any given time, the mass of gas with $\rho\gtrsim 10^{-15}$ g cm*-3* in an ideal model is slightly lower than its non-ideal counterpart, which is reasonable since magnetic fields resist gravitational collapse, and the non-ideal processes weaken the local magnetic field. This decrease of mass in dense gas, however, occurs mainly on small scales, as discussed in Section 4.2 below.

The ionisation fraction decreases for increasing density within the cloud, reaching fractions of $n_{\text{e}}/n\sim 10^{-9}$ in filaments of $\rho\approx 10^{-16}$ g cm*-3*, and fractions of $n_{\text{e}}/n\sim 10^{-14}$ in protostellar discs. Hence we expect non-ideal MHD to act mainly on small scales.

4.2 Small scale structure

Fig. 5 shows a $0.06\times 0.06$ pc2 region of each model to highlight differences caused by changing the initial magnetic field strength, or by using ideal instead of non-ideal MHD. Each panel shows the same spatial region for comparison. Since each model was initialised with the same velocity field, the differences are a direct result of the initial magnetic field strength and non-ideal MHD processes.

In the hydrodynamic and weak field ($\mu_{0}=10$ and $20$) models, Fig. 5 shows a single dense region in each panel that does not strongly connect to any other large scale structure through dense filaments. Although the clump in the hydro model (right panel) appears filamentary in the figure, it is in fact pancake-shaped. The structure in these models forms from the collapse of a short filament, and rotates around the stars as the entire substructure moves throughout the cloud. Each of these dense regions hosts or has hosted several stars, and there are several filaments emanating from these dense regions.

The small scale structure of the stronger field models ($\mu_{0}=5$ and $3$, left two columns) are more filamentary, with several small clumps embedded within the filament; each individual clump is smaller than in the weakly magnetised models. There are fewer filaments in these models, and the ones that exist tend to be denser and narrower than in their weakly magnetised counterparts.

On this scale, there appears to be a trend from filamentary to broken filaments to clumpy as the initial magnetic field strength decreases. As suggested by the large-scale structure, the magnetic field is approximately perpendicular to the dense filaments.

Discs are visible around many of the stars, and these will be discussed in Section 4.4 below.

4.2.1 Non-ideal vs ideal MHD

The differences caused by the inclusion of non-ideal effects are primarily found on small scales. As an example, Figs. 6 and 7 show the evolution of a small region of the gas column density and the density-weighted line-of-sight averaged magnetic field strength (i.e. $\left<|B|\right>=\int|B|\rho\text{d}y^{\prime}/\int\rho\text{d}y^{\prime}$), comparing N05 (top rows) to I05 (bottom rows). In this region, the gas structure is indistinguishable between the two models at $t=1.16$ t${}_{\text{ff}}$. At 1.41 t${}_{\text{ff}}$, the primary filament is slightly more clumpy in I05, and by 1.49 t${}_{\text{ff}}$, the structure of the primary filament has noticeable differences, as indicated by where the stars form in the two models. This again confirms that non-ideal MHD affects the density on small scales, and that there is a greater divergence between an ideal/non-ideal pair as the simulations continue to evolve. Throughout the evolution, there is slightly more dense gas in the non-ideal models compared to their ideal MHD counterparts (c.f. Fig. 2 and associated discussion), although this is not necessarily clear from the image.

Fig. 7 shows that the magnetic field is stronger in the dense regions around the stars (i.e. the discs) compared to the background. Except at 1.49 t${}_{\text{ff}}$, there is no noticeable difference in strength between the ideal and non-ideal MHD models; however, this stronger magnetic field at 1.49 t${}_{\text{ff}}$ corresponds to a region of slightly higher density in I05 (final column of Fig. 6). This is consistent with Fig. 4 which suggest that all models have similar mean magnetic field strengths at a given gas density and an order of magnitude scatter.

4.3 Stellar populations

Sink particles with a radius of 0.5 au are inserted when $\rho_{\text{max}}=10^{-5}$ g cm*-3* is reached. Given our sink radius and insertion criteria, each sink particle represents one star or brown dwarf, thus we use ‘star’ and ‘sink particle’ interchangeably.

4.3.1 Star formation rate

Fig. 8 shows the total number of stars and the total mass in stars as a function of time. Table 1 lists information about the stellar and disc populations at $t=1.45$ t${}_{\text{ff}}$ and at the end of the respective simulations, $t_{\text{final}}$.

By $t=1.45$ t${}_{\text{ff}}$, there are 8–19 stars, depending on the simulation. The hydro model has the highest number of stars, which is to be expected since it lacks magnetic support against gravitational collapse. However, for the magnetised models, the number of stars does not appear to correlate with the initial magnetic field strength or whether or not we include non-ideal MHD. The magnetised model with the highest number of stars is N05, which has a relatively strong initial magnetic field. Caution is required when comparing stellar populations at a given time, since, for example, N05 and I05 are undergoing an epoch of rapid star formation over the $\approx 0.1$ t${}_{\text{ff}}$ before they are stopped.

The integrated star formation rate (i.e., the total mass in stars as a function of time; bottom row of Fig. 8) depends on the initial magnetic field strength in a (nearly) systematic way. The hydro model always has the most mass in stars, followed by the models with $\mu_{0}=20$, $10$ and $5$, respectively. Thus, amongst these models, we find that stronger magnetic fields slow the star formation rate, which is in agreement with Price & Bate (2008, 2009) and the analysis of the dense gas presented above. We further find that the star formation rate does not depend on whether ideal or non-ideal MHD is assumed.

The models with the strongest initial magnetic field strengths (N03 and I03) do not follow the above trend. At any given time $t\gtrsim 1.1$ t${}_{\text{ff}}$, the integrated star formation rate is greater in these models than the other magnetised models and only slightly lower than Hyd. The strong magnetic field qualitatively changes the collapse, with a quicker collapse along the magnetic field lines giving a larger quantity of dense gas (recall Fig. 2). Despite the strong magnetic fields, this larger quantity of dense gas fuels accretion onto the stars. This deviation from the systematic trend is in contradiction with Price & Bate (2008, 2009), however, it may depend on the initial velocity field.

Many stars remain near the clump in which they were born, however, several stars are ejected into the low-density medium due to interactions with other stellar systems; these stars are typically ejected with low ($v\lesssim 5$ km s*-1*) velocity, and we observe only four high-velocity ($v\gtrsim 5$ km s*-1*) stars by $t_{\text{final}}$. Not all stars accrete at similar rates, and by $t_{\text{final}}$, several stars in each model have stopped accreting gas. Most of these non-accreting stars have been ejected into the low-density medium and are amongst the lowest mass stars in our suite. Despite the inclusion of magnetic fields and our low number of stars, we observe ejected and non-accreting stars, as seen previously in the literature (e.g. Bate et al., 2003; Bate & Bonnell, 2005).

4.3.2 Stellar mass distributions

Fig. 9 shows the number of stars in each mass bin at both $t=1.45$ t${}_{\text{ff}}$ (red hatched) and at $t_{\text{final}}$ (blue outline). Given the small number of stars and that the majority of them are still accreting, it is not practical to calculate an IMF, or explicitly comment on stellar population.

Considerably fewer stars are produced in our simulations compared to the magnetised barotropic calculations of Price & Bate (2008), which produced 17–23 and were evolved to 1.27–1.53 t${}_{\text{ff}}$. This is because introducing radiative feedback reduces small-scale fragmentation (Bate, 2009b). As expected, the numbers are more similar to those produced by the magnetised radiation hydrodynamical calculations of Price & Bate (2009), which produced 3–10 objects and were evolved to 1.36–1.54 t${}_{\text{ff}}$. These earlier magnetised star formation calculations used Euler potentials to model ideal MHD, and the earlier radiative transfer scheme of Whitehouse & Bate (2006) which did not include the diffuse interstellar medium treatment and separate gas and dust temperatures. Despite these differences, the results are qualitatively similar.

4.3.3 Multiplicity

The majority of the stars are born as single stars and then become part of multiple systems when they become gravitationally bound during a close interaction (similar to, e.g., Bate, 2012, 2018; Seifried et al. 2013). Many of these stars are born $\sim\mathcal{O}\left(10^{2}-10^{3}\right)$ au and $\sim\mathcal{O}\left(0.1-10\right)$ kyr apart, implying that most of our binary systems are not primordial.

We locate systems hierarchically up to systems containing four mutually bound stars; this is performed simultaneously with determining disc properties, as described in Section 4.4. First, for each star+circumstellar disc, we determine if it is bound to its closest star+circumstellar disc; if so, then this is a binary pair. Once all the binaries are found, we determine if each binary is bound to its closest star or binary system. For the systems with three stars, we search to determine if the nearest star is single and also bound, for a system of four stars. In our clusters, higher order systems typically include one or two binaries. The number of stellar system of each population is listed in the $N_{\text{systems}}$ column in Table 1. Although there are more systems with only one star than higher order systems, these single-star systems represent $\lesssim 45$ per cent of the total number of stars in any given model.

Excluding mergers and one lone exception, once a binary is formed, it persists for the duration of the simulation. The higher order systems, however, are created and destroyed with time, which is particularly noticeable in N05 and I05, which both form many new stars after $1.45$ t${}_{\text{ff}}$.

This high number of multiplicity is consistent with previous hydrodynamic studies (e.g. Bate et al., 2003; Bate, 2009b, 2012, 2019), suggesting that magnetic fields do not significantly hinder or promote the formation of binaries or higher order systems.

4.3.4 Mergers

When two stars come within 27 R⊙, we considered them to have merged (we increased this from the default of 6 R⊙ for computational efficiency since following inspiralling stars requires a small computational timestep). During a merger, we replace the two sink particles with a single sink with the centre of mass properties of its progenitors. The $N_{\text{mergers}}$ columns of Table 1 lists the number of mergers prior to $1.45$ t${}_{\text{ff}}$ and $t_{\text{final}}$.

There are two classes of mergers in our simulations: stars that ‘collide’ during fly-by and merge, and those that form a binary system with a decaying major axis and ultimately merge. In both I03 and I20, the stars form a binary system before merging, whereas the remainder of the mergers occur during a fly-by (where a fly-by is indistinguishable from a binary with a large major axis and high eccentricity).

In I03, the stars were born with a separation of $\sim$260 au, and quickly became gravitationally bound. Their orbit quickly decayed, and within 14 kyr of the birth of the younger star, the periastron separation decreased to $\sim$27 R⊙ and the stars merged. At the time of merger, their apastron was $\sim$110 R⊙.

In I20, the binary system initially has a stable average separation of $\sim$5 au before a close encounter causes the orbit become more eccentric, but maintaining a similar apastron. However, a second encounter causes the orbit to decay while becoming more eccentric until the stars ultimately merge at periastron. This is a merger of two stars of similar masses. Then, $\sim$2.5 kyr later, a low mass star ‘collides’ with this star in a fly-by and merges. Thus, between the two mergers, the primary star rapidly increases its mass, and at $1.45$ t${}_{\text{ff}}$ is the second most massive star in our suite (the most massive star is the star that underwent a merger in I03).

In reality, due to the large merger radius, many of these systems may have formed spectroscopic binaries rather than merger products.

4.4 Protostellar disc properties

We track three separate classifications of discs: circumstellar discs (i.e. a disc around a single star), circumbinary discs (i.e. a common disc around a binary star system), and circumsystem discs (i.e. discs around a system of either three or four stars). Thus, it is possible for several discs to be associated with each star. The number of each of these classifications of discs is listed in the $N_{\text{discs}}$ column of Table 1, where the total number of circumstellar and circumbinary discs is listed in parentheses and the number not in parentheses lists the number of discs associated with the highest order systems from the $N_{\text{systems}}$ column.

To calculate the mass of the discs, we follow the prescription of Bate (2018). For each star, we first sort all the particles by distance. We check the closest particle to determine if it is bound to the sink, has an eccentricity of $e<0.3$ and has density111Unlike Bate (2018), we include the density threshold to prevent our discs from including any low-density filamentary material. $\rho\geq 10^{-14}$ g cm*-3*; if so, it is added to the sink+circumstellar disc system. Each consecutive particle is checked and added to the system if the criteria are met. We do this for all the particles within 2000 au of the star or until another star is encountered. To determine higher order discs, we determine if stellar systems (star+circumstellar discs) are mutually bound, as described in Section 4.3.3. If so, then the above process is repeated, but using the bound pair rather than star+circumstellar disc. We record the circumbinary disc independently from the circumstellar discs for later analysis. This process is repeated with each new system up to a maximum of four stars per system. The mass of each disc is the total mass of the gas that has been added to the disc. The radius of each disc, $R_{\text{disc}}$, is the radius that includes 63.2 per cent of the disc mass.

4.4.1 Formation history

Stars tend to form in isolation (recall Section 4.3.3), thus the initial discs are small, circumstellar discs, assuming a primordial disc forms at all. However, the bulk flow of the gas promotes interactions between the stellar systems. As such, the systems are frequently promoted to or demoted from higher order systems due to these interactions. As stars approach and become bound, their discs usually merge forming a larger circumsystem disc while often retaining small circumstellar discs. Although disc growth is a common consequence, interactions with unbound stars (i.e. fly-by interactions) or bound stars on large orbits can cause a reduction in the disc size through tidal stripping, or totally disrupt the disc. Thus, the masses and radii of our discs are rapidly changing, as previously seen in the hydrodynamical simulations by Bate (2018). We thus focus our analysis on the discs instantaneously existing at $1.45$ t${}_{\text{ff}}$, similar to observing discs in a particular star forming region at a particular time.

Fig. 10 shows the column density of two disc-containing regions of each model; each panel has the same spatial scale, and we show the same region for each each ideal/non-ideal pair (left and right panels in each column, respectively). From a visual comparison, we observe the strong influence of the non-ideal MHD processes on the small-scale evolution.

Many of the panels in Fig. 10 show isolated circumstellar/circumsystem discs (e.g. the left-hand panels of the $\mu_{0}={10}$ and $20$ models and Hyd) or bound systems (e.g. the mutually bound discs in the right-hand panels of the $\mu_{0}={3}$ and $5$ models). At this time, these discs are relatively smooth, and a reasonable analysis of the disc properties can be performed.

The discs in the left-hand panels of N03 and N05 are undergoing a violent change as a result of their recent interactions. More dramatically, disc interaction is occurring at this time in the right-hand panels of the $\mu_{0}={20}$ models and Hyd. A smooth disc in I20 became disrupted at $t\approx 1.39$ t${}_{\text{ff}}$ yielding the disrupted system in the figure, while gas has been stripped from one disc and accreted onto the other in N20. Thus, as discussed by Bate (2018), the formation history of discs in a cluster is violent.

Fig. 11 shows the gas density and magnetic field strength and direction in a slice through the centre of the most well-defined circumstellar disc in each model.

These representative discs indicate that there is no trend amongst the models. It is worth explicitly noting that large discs exist in each model, even those clusters with initially strong magnetic field strengths. This suggests that the angular momentum required for disc formation originates from the turbulent velocity of the gas, and that any hindrance of disc formation by the magnetic field (i.e., magnetic braking) is relatively weak; this will be discussed further in Section 4.4.5. The smaller discs are a result of their environment and the interaction with other stellar systems, rather than a dependence on the magnetic field; for example, the discs in Hyd and N03 are orbiting a circumbinary disc which regulates the sizes of both discs. The larger discs have not recently undergone any interaction with another stellar system, which allows their discs to grow. This further suggests that discs formation and evolution is more strongly dependent on the local velocity field than the local magnetic field.

The majority of the discs are not as well-defined as in Fig. 11. Several single stars do not have circumstellar discs; many of these are low mass stars that have been ejected from their birth clump, thus were likely stripped of their primordial disc (if they even had one) and are not in an environment that is gas-rich enough for the disc to reform. Many of the circumsystem discs are being influenced and disrupted by the orbits of their host stars, and several are undergoing interactions where the discs is being rapid augmented or destroyed. For further discussion, see Bate (2018); for the complete disc population in our study, see Appendix B.

The second row of Fig. 11 does not reveal any jets or outflows, which is true for all the discs at $1.45$ t${}_{\text{ff}}$. This is likely a consequence of the constant changes in disc orientation due to close encounters, and our limited numerical resolution preventing disc wind formation.

4.4.2 Disc sizes

Fig. 12 shows the disc radius, mass and disc-to-stellar mass ratios of the highest-order discs at $1.45$ t${}_{\text{ff}}$. The radius is for the highest order disc in each system (i.e. the disc that surrounds all the stars in the bound stellar system; this corresponds to the unbracketed numbers in Table 1), while the mass is a sum of the masses of all the discs in the system; thus, there may not necessarily be a corresponding radius to each mass if the stellar system does not include a circumsystem disc (e.g. I20).

Circumstellar discs of $10\lesssim r/\text{au}\lesssim 80$ and circumsystem discs of $30\lesssim r/\text{au}\lesssim 500$ form in our models with no obvious dependence on our simulation parameters. Several $10\lesssim r/\text{au}\lesssim 30$ circumstellar discs form even in our strongest magnetic field models (N03 and I03), confirming that the examples in Fig. 11 are a good representation.

We must be cautious about resolution effects in our disc analysis. The mass of an SPH particle is $m_{\text{p}}=10^{-5}$ M⊙, thus our discs typically contain $10^{2}-10^{4}$ particles. Despite the small numbers, they are resolved as per the Jeans mass (Bate & Burkert, 1997) and the Toomre-Mass (Nelson, 2006) criteria. However, Nelson (2006) suggests that the scale height must be resolved by at least four smoothing lengths in the midplane of the disc, in order to prevent numerical fragmentation. Our discs do not meet this criteria, but we also do not observe any disc fragmentation.

4.4.3 Disc masses

As with the radius, disc masses show no obvious dependence on simulation parameters. At $1.45$ t${}_{\text{ff}}$, there are 10 massive discs ($M>0.1$ M⊙), only one of which is a circumstellar disc. Many of these massive discs survive to the end of the simulation, however, the survival of massive discs is not guaranteed. Throughout the evolution of the clusters, there are several discs with $M>0.1$ M⊙, which are then partially or totally disrupted by close encounters. This is to be expected since these massive discs typically form in high-density regions into which new stars are gravitationally attracted, and are typically extended discs whose outer regions are easy to strip off.

Slightly more than half of the discs have masses of $M_{\text{disc}}>0.1M_{\text{star}}$, where $M_{\text{star}}$ is the total stellar mass in the system, with the higher order discs typically having smaller ratios. We also find that most binaries are comprised of nearly equal mass stars.

4.4.4 Magnetic field strength and geometry

The structure of the magnetic field of the representative circumstellar discs is shown in the third and fourth row of Fig. 11, and the arrows in the fifth row indicate the direction of the magnetic field. Despite the reasonably smooth density profiles (top row), spiral structures exist in the magnetic field (see discussion in Wurster et al., 2018c), resulting in a range of magnetic field strengths in each disc. The average magnetic field strength of each disc and the range containing 95 per cent of the field strengths is show in the top row of Fig. 13.

The average field strength is $0.005\lesssim B/\text{G}\lesssim 0.5$, with the range in any given disc spanning up to $\sim$1.5 dex. The disc magnetic field strength is independent of the initial field strength, and the non-ideal processes moderate the magnetic field strengths in the discs, resulting in a narrower range of field strengths in the non-ideal suite than the ideal suite.

The middle and bottom rows of Fig. 13 show the fraction of the magnetic field component in the discs that is parallel to the disc axis (middle row) and in the azimuthal direction (bottom row). The average magnetic field components are calculated by

[TABLE]

where $j\in\{\phi,z\}$, and we individually sum over the $N$ positive and $n$ negative values before performing a linear weighted average of the two terms. For the circumbinary and circumsystem discs, we include the contribution from the lower order discs. Unsurprisingly, there is a strong azimuthal component in most of the discs.

Although the initial magnetic field is $-B_{0}\hat{\bm{z}}$, there is no preference for aligned or anti-aligned magnetic fields, thus the Hall effect is expected to be important in approximately half of the discs (Tsukamoto et al., 2015b); however, other effects may be dominant such that the contribution from the Hall effect is minimal (Wurster & Bate, 2019).

When comparing the poloidal to toroidal magnetic field in the region around the disc, we find that the poloidal component is typically dominant, which is consistent with studies of isolated star formation where the initial magnetic field is (anti-)aligned with the rotation axis (e.g. Bate et al., 2014; Wurster et al., 2018a). The larger-scale magnetic field (out to $\sim$800 au) typically reflects the small-scale magnetic field (out to $\sim$200 au), although in a few cases, a horizontal field twists near the disc to become vertical.

4.4.5 Is there a magnetic braking catastrophe?

The presence of large ($r\gtrsim 20$ au), massive ($M\gtrsim 0.01$ M⊙) discs in all of our calculations suggests that their formation and structure primarily depends on the local environment. Disc formation occurs regardless of the initial magnetic field strength of the progenitor cloud and whether or not non-ideal MHD is employed. As shown in Fig. 4, the magnetic field strength in the dense gas ($\rho\gtrsim 10^{-17}$g cm*-3*) is similar for all models, but varies by $\sim$1 dex at each density. Thus, the local magnetic field strength in star forming clumps does not simply reflect that of the initial strength in the cluster. As a result, large rotationally supported discs form even in ideal MHD models with initially strong magnetic fields.

This suggests that the magnetic braking catastrophe is not a problem in realistic environments due to the turbulent and dynamic nature of the environment and that turbulence and interactions are more important than magnetic fields for disc formation and evolution; this is in agreement with Seifried et al. (2013). Our results also suggest a universal initial condition in local star forming regions (or at least a limited set of initial conditions) that may even be independent of the large scale magnetic field of its host environment. Further investigation of this is beyond the scope of this study, but this possibility should be explored in the future.

Several single stars in our suite do not have circumstellar discs, however, investigating their history suggests that this is a result of their dynamical ejection into a gas-poor environment, rather than as a result of the magnetic fields and magnetic braking catastrophe. Thus, it is likely that nearly all stars host discs, even if only briefly.

Although we have learned much from investigating the magnetic braking catastrophe in idealised simulations of the formation of isolated stars, it appears the main problem is the artificial initial conditions that were employed.

5 Discussion

5.1 Previous cluster scale calculations

Magnetised hydrodynamical simulations of star formation in stellar groups and clusters have lagged up to a decade behind purely hydrodynamical simulations. The first such calculations that attempted to resolve individual stellar systems were the ideal MHD simulations of Price & Bate (2008) using a barotropic equation of state, and Price & Bate (2009) that included radiative transfer and a realistic gas equation of state. These SPH simulations employed Euler potentials to model the magnetic fields, limiting the magnetic field geometries that could be modelled. However, they successfully demonstrated the two main effects of magnetic fields on star cluster formation — the formation of different structures in magnetised molecular clouds compared to those modelled purely with hydrodynamics (e.g. anisotropic turbulent motions, low-density striations, and magnetised voids), and the star formation rate is reduced with increasing magnetic field strength due to magnetic support on large scales. Subsequent grid-based calculations have confirmed the reduction of the star formation rate with magnetic fields (e.g. Padoan & Nordlund, 2011; Padoan et al., 2012; Federrath & Klessen, 2012; Myers et al., 2014), but the magnitude of the effect is only at the level of factors of two or three times slower with strong fields compared to the star formation rate obtained using pure hydrodynamics.

Regarding the effects on stellar masses, the situation is less clear. The calculations of Price & Bate (2008, 2009) produced only small numbers of stars and brown dwarfs and did not find any strong evidence for an effect of magnetic fields on their characteristic mass. The calculations of Myers et al. (2014) produced $\approx 90$ objects each, and they found that strong fields may increase the characteristic mass by a factor of two or three. However, the more recent calculations of Cunningham et al. (2018) which also include protostellar outflows and produce $\approx 20$–$70$ objects each find no evidence for a shift in the characteristic mass for calculations performed with driven turbulence, and a weak trend for a characteristic mass that decreases with increasing magnetic field strength for calculations performed with decaying turbulence. Thus, with the limited statistics currently available, there is no clear dependence of stellar masses on magnetic field strength. Our results do not shed any further light on this question — the numbers of objects formed are too small to detect any weak dependence of the stellar mass function on magnetic fields, should one exist. It seems certain, however, that magnetic fields do not play the dominant role in determining the IMF. Instead, the characteristic stellar mass seems to be set primarily by radiative feedback from protostars (Bate, 2009b, 2012, 2014, 2019; Krumholz, 2011; Krumholz et al., 2012) and by thermodynamic processes (Lee & Hennebelle, 2018b).

The calculations we present here extend the investigation of magnetised star formation in two main ways beyond the studies mentioned above. First, they are the first calculations of the formation of stellar groups that treat the three main non-ideal MHD effects (ambipolar diffusion, the Hall effect, and Ohmic diffusion). We find that the non-ideal processes have little effect on the large-scale structures ($\gtrsim 0.05$ pc). However, including non-ideal MHD does alter structures on smaller scales (i.e., the filaments and cores) and, thus, changes the details of the cloud fragmentation and distribution of protostars. There is no evidence from our current calculations that including non-ideal MHD leads to any changes in the statistical properties of stellar systems, but calculations that produce at least an order of magnitude more objects will be required to investigate this question further. Second, unlike in most of the above MHD calculations, protostellar discs are resolved. The calculations of Price & Bate (2008, 2009) resolved discs as small as $\sim 10$ au, but they used Euler potentials to model the magnetic fields which did not allow magnetic fields to be wound up. Our calculations employ sink particles with accretion radii of only 0.5 au, meaning that binary and multiple protostellar systems are resolved, along with protostellar discs with radii as small as $\approx 10$ au.

5.2 Discs and binaries

The effects of both ideal and non-ideal magnetic fields have been thoroughly investigated in studies of isolated star formation (e.g. Allen et al., 2003; Price & Bate, 2007; Hennebelle & Fromang, 2008; Mellon & Li, 2008; Tomida et al., 2010a, b; Commerçon et al., 2010; Seifried et al., 2011; Tomida et al., 2013; Bate et al., 2014); see Wurster & Li (2018) for a recent review. Early studies found that large, protostellar discs did not form in the presence of strong, ideal magnetic fields — the so-called ‘magnetic braking catastrophe’ (Allen et al., 2003; Galli et al., 2006; Price & Bate, 2007; Hennebelle & Fromang, 2008). These models assumed simple initial conditions, with ordered initial magnetic fields aligned parallel to the rotation axis. However, with misaligned magnetic fields (e.g. Joos et al., 2012; Lewis et al., 2015; Lewis & Bate, 2017), turbulent velocity fields (e.g. Machida & Matsumoto, 2011; Seifried et al., 2012, 2013; Santos-Lima et al., 2013; Joos et al., 2013; Lewis & Bate, 2018) or non-ideal MHD (e.g. Tsukamoto et al., 2015a, b; Wurster et al., 2016; Tsukamoto et al., 2017; Vaytet et al., 2018; Wurster et al., 2018a, b, c), substantial protostellar discs can be formed. The key differences between these simulations and the earlier more idealised models are that the discs are subject to weaker magnetic braking and angular momentum transport due to complex field/rotation geometries and/or weaker magnetic fields (e.g. due to non-ideal effects and/or turbulent or grid-scale magnetic reconnection).

Similarly, past studies of the collapse of isolated magnetised molecular cloud cores have shown that magnetic fields can inhibit fragmentation and multiple star formation (Hosking & Whitworth, 2004; Price & Bate, 2007; Hennebelle & Teyssier, 2008; Commerçon et al., 2010, 2011). However, Wurster et al. (2017b) and Wurster & Bate (2019) showed that whether or not a dense core fragments depends on its initial density (sub-)structure, rotation, and field geometry rather than just on the strength of the magnetic field. Thus, in a turbulent molecular cloud in which dense cores are likely to have considerable substructure and rotation and fields that are misaligned with the net rotation, fragmentation is quite possible even with strong fields.

Seifried et al. (2013) performed calculations of ideally magnetised turbulent molecular cloud cores with masses up to 1000 M⊙, that produce multiple protostellar objects and resolved discs. However, their initial conditions were chosen to be strongly centrally-condensed ($\rho\propto r^{3/2}$) which favours massive star formation, rather than to study the formation of low-mass stellar groups. The overwhelming majority of their stars formed via gravitational collapse of distinct overdense regions, then multiple systems formed via capture rather than disc fragmentation; their discs tended to have masses $0.05\lesssim M/$M⊙$\lesssim 0.1$ and radii of $50\lesssim r/\text{au}\lesssim 150$. In general, similar to the work we present here, their study concluded that discs form as a result of the turbulence within star-forming regions, independent of the the presence of strong magnetic fields, which quickly become disordered.

In our calculations, protostars are surrounded by discs with a variety of sizes, ranging from $\approx 10-500$ au. The multiple systems form from the fragmentation of filaments, rather than disc fragmentation. Similar results are obtained from turbulent radiation hydrodynamical simulations without magnetic fields (Bate, 2012, 2018, 2019). The statistical properties of the multiple systems and discs do not display any obvious trends with initial magnetic field strength, or with whether or not non-ideal MHD processes are included. Thus, it seems that it is the density and velocity structure and dynamical interactions between protostars that dominate the properties of multiple systems and discs, rather than whether or not magnetic fields are present.

6 Summary and conclusion

We have presented a suite of radiation non-ideal magnetohydrodynamic simulations modelling star formation in a 50 M⊙ molecular cloud. Each model was initialised with the same turbulent velocity field, and was threaded by a uniform magnetic field; we tested four magnetic field strengths plus a purely hydrodynamic model. For each magnetised model, we performed both an ideal and non-ideal MHD version, where the latter included Ohmic resistivity, ambipolar diffusion and the Hall effect. Sink particles with radii of 0.5 au were inserted at the opacity limit for fragmentation such that each particle represented one star.

Our key results are as follows:

Magnetic fields influence the large scale structure. The initial mass-to-flux ratio was found to mainly affect the large scale evolution, with models with stronger magnetic fields slower to form high-density clumps due to increased magnetic support. However, our models with the strongest magnetic field defied this trend, with gas rapidly collapsing along the magnetic field lines. With strong magnetic fields, the field lines tend to be perpendicular to dense filaments, while with weak magnetic fields there is a non-trivial parallel component near the filament. Magnetic fields tend to be parallel to low-density filaments. 2. 2.

Initial conditions are quickly erased. Magnetic field strengths in the high-density gas ($\rho\gtrsim 10^{-17}$ g cm*-3*) were found to be independent of the initial mass-to-flux ratio and of the non-ideal MHD processes. At a given high density, the magnetic field strength and range was similar for each model and spanned approximately one order of magnitude. Thus, individual star forming clumps contain a wide range of magnetic field strengths. 3. 3.

Non-ideal MHD acts mainly on small scales. Small scale structures tend to be filamentary for models with strong initial magnetic fields and more clumpy for models with weak initial magnetic fields. Non-ideal MHD, while not affecting the large scale evolution, was found to play a more significant role on $\lesssim 0.05$ pc scales. We found the main effect was to create tighter filaments and to reduce the magnetic field strength in gas in and around the protostellar discs. 4. 4.

Star formation is violent and chaotic. Binary and hierarchical triple/quadruple systems form in our models. Binary stars form primarily by gravitational capture, and nearly all of the binaries survived to the end of the simulation, with their orbits and ellipticities strongly influenced by interactions with other stars. The number of stars does not depend on the initial magnetic field strength, however, there is a general trend of decreasing total stellar mass as the initial magnetic field strength is increased; our strongest field models defied this trend. 5. 5.

There is no magnetic braking catastrophe. Protostellar discs with radii of \sim$$10-80 au form in all of our models, and circumsystem discs have radii up to $\sim$500 au. The magnetic field strength in the protostellar discs is independent of the initial magnetic field strength of the molecular cloud, and the non-ideal processes yield a narrower range of disc field strengths than when assuming ideal MHD. Thus, non-ideal MHD may not be required for disc formation, but it does regulate the magnetic field in discs. The presence of large, massive discs in every model suggests that the magnetic braking catastrophe only arises when modelling the formation of isolated stars from idealised initial conditions, and that in the presence of strong global magnetic fields, turbulence promotes disc formation.

In all of our models, we form stars, stellar systems, and protostellar discs, even in clusters with initially strong magnetic field strengths. Given the weak or non-existent trends in the properties that we investigated, this suggests that observed objects (and possibly trends) can be reproduced, even when including realistic strong magnetic fields. Thus, magnetic fields are not a hinderance to numerical star and protostellar disc formation.

Acknowledgements

We would like to thank the referee, Dr Daniel Seifried, for useful comments that greatly improved the quality of this manuscript. JW and MRB acknowledge support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007- 2013 grant agreement no. 339248). DJP received funding via Australian Research Council grants FT130100034, DP130102078 and DP180104235. Calculations and analyses for this paper were performed on the University of Exeter Supercomputer, Isca, which is part of the University of Exeter High-Performance Computing (HPC) facility, and on the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. The research data supporting this publication will be openly available from the University of Exeter’s institutional repository. Several figures were made using splash (Price, 2007).

Appendix A Resolution study

We perform two high-resolution versions of our models with $\mu_{0}={5}$. For these models, we employ a mass resolution seven times higher than our fiducial resolution models, yielding an $\approx 2$ times improvement in spatial resolution. These models are named N05**h and I05**h for the non-ideal and ideal MHD versions, respectively, and are computed to $t=1.356\approx 1.36$ t${}_{\text{ff}}$ and $t=1.446\approx 1.45$ t${}_{\text{ff}}$, respectively.

A.1 Structure

Fig. 14 shows the global structure of the models in our resolution comparison at the final times of the high resolution models. We find qualitative agreement between the two resolutions, and as expected, the small scale filamentary structure appears better resolved in the high resolution models.

The differences can be better seen on smaller scales, as shown in Fig. 15. This figure shows the better defined structures, with narrower and less wispy filaments in the high resolution models. In several regions, there are denser, more well-defined clumps in the high resolution models compared to similar regions in fiducial resolution counterparts. Given the chaotic nature of the evolution, as time progresses, the differences between the resolutions becomes more apparent. This also suggests that the high resolution clusters evolve slightly faster, since their higher densities results in higher gravitational accelerations.

At $1.36$ t${}_{\text{ff}}$, the qualitative structure is similar between N05 and I05. Although N05**h and I05**h are also similar at this time, there are a few more differences, including two high-density clumps in N05**h and one in I05**h that do not appear in the other model. This suggests that non-ideal effects are more important at higher resolutions, which is reasonable since artificial resistivity is resolution-dependent; there is less resistivity in I05**h compared to I05, while N05**h and N05 should have similar amounts in dense regions since physical resistivity is expected to dominate over artificial resistivity in these regions, and should be resolution-independent.

A.2 Stellar populations

Our star formation prescription does not depend on resolution — we still insert a sink of radius 0.5 au when $\rho_{\text{max}}=10^{-5}$ g cm*-3* is reached. Fig. 16 shows the total number of stars (top panel) and the total mass in stars (bottom panel) as a function of time; Table 2 lists the stellar population properties at $t\approx 1.36$ and $1.45$ t${}_{\text{ff}}$.

Star formation begins later in the high resolution models, with higher-resolution models collapsing slower than their lower-resolution counterparts (as previously discussed in Wurster et al., 2018c). Once the star formation begins in I05**h, the star formation rate is similar to that of I05, and both models have accelerating star formation rates. Thus, although initial star formation may be delayed in the high resolution models, similar numbers of stars ultimately form.

In the non-ideal models, there is a similar delay of the onset of star formation at high resolution. At $1.36$ t${}_{\text{ff}}$, the high resolution model has fewer stars and less total stellar mass but this is simply caused by the later start. Starting at $t\approx 1.36$ t${}_{\text{ff}}$, N05 undergoes a rapid star formation epoch. Given the high-density clumps in Fig. 16, it is likely that N05**h is also about to undergo a rapid star formation epoch.

For all time, the total stellar mass is lower in the high resolution model, which is mainly a result of the delay in the collapse due to better resolved pressure gradients. High-density circumstellar material is also accreted by sink particles at a lower rate due to the lower numerical viscosity (see the Appendix of Bate, 2018).

For the high resolution non-ideal MHD model at $1.36$ t${}_{\text{ff}}$, all stars are single stars, but this is simply because only four stars have formed by this point. In the ideal MHD model at $1.45$ t${}_{\text{ff}}$, both ideal models have a similar number of stars and higher order systems (Table 2).

A.3 Protostellar discs

At $1.36$ t${}_{\text{ff}}$, every star has a circumstellar disc in the high resolution models, and higher order discs form at later times in I05**h. Given the chaotic nature of the simulations, each fiducial resolution disc does not have a high resolution counterpart, thus we are unable to perform a direct disc-to-disc comparison. Fig. 17 shows the gas density and magnetic field strength and direction in a slice through the centre of the the most well-defined circumstellar disc in each model. At both resolutions the magnetic field retains a spiral structure in the discs (third row of Fig. 17), and the azimuthal component remains the largest component.

The high resolution discs generally have smaller radii and less mass than those in the fiducial resolution models. These results are to be expected in the ideal MHD calculations as a consequence of the factor of $\sim$4 decrease in the artificial resistivity by increasing the resolution. The higher artificial resistivity in the fiducial resolution models will artificially promote disc formation (Wurster et al., 2016), increasing their size and decreasing their magnetic field strength. With non-ideal MHD, the non-ideal processes should dominate over artificial resistivity. However, because circumstellar discs tend to grow in both mass and size with time (see also Bate, 2018), the difference in size and mass may simply be due to the delay in the onset of the star formation (i.e., the discs have had less time to grow in the high resolution simulations).

Large protostellar discs still form in both of our high resolution simulations, again suggesting that the disc formation is promoted by the turbulent velocity field and that the magnetic braking catastrophe is a result of idealised initial conditions. However, given that we do see a tentative resolution dependence on disc size and mass, we cannot completely rule out the catastrophe. It is possible that, given enough resolution, strong magnetic fields may again hinder disc formation even in turbulent gas. To study this properly would require even higher resolution simulations, which are currently computationally infeasible.

Appendix B Protostellar disc populations

Figs. 18 and 19 show a 412 $\times$ 412 au2 region around every star in our fiducial resolution simulations at $1.45$ t${}_{\text{ff}}$, and at the times discussed in Appendix A for the high-resolution models; the images are rotated such that the disc is face-on and show the gas column density and density-weighted line-of-sight-averaged magnetic field strength, respectively (i.e. $\left<|B|\right>=\int|B|\rho\text{d}z^{\prime}/\int\rho\text{d}z^{\prime}$). These figures differ slightly from the slice through the centre of the star as in Figs. 11 and 17 to better show how the discs interact with their surroundings, by also showing non-coplanar, nearby objects. Since we plot every star, one may observe the wide variety of disc structures, multiplicities, and the stars that do not form discs.