Disc formation and fragmentation using radiative non-ideal magnetohydrodynamics

TL;DR

This study uses 3D radiative non-ideal MHD simulations to explore how magnetic fields, rotation, and initial conditions influence disc formation and fragmentation during gravitational collapse.

Contribution

It provides a comprehensive analysis of disc formation and fragmentation across a wide parameter space, highlighting the role of non-ideal MHD effects and magnetic field orientation.

Findings

Fragmentation occurs more with faster rotation and weaker magnetic fields.

Non-ideal MHD promotes disc formation and fragmentation in strong magnetic fields.

Bipolar outflows are associated with strong, aligned magnetic fields.

Abstract

We investigate the formation and fragmentation of discs using a suite of three-dimensional smoothed particle radiative magnetohydrodynamics simulations. Our models are initialised as 1M$_\odot$ rotating Bonnor-Ebert spheres that are threaded with a uniform magnetic field. We examine the effect of including ideal and non-ideal magnetic fields, the orientation and strength of the magnetic field, and the initial rotational rate. We follow the gravitational collapse and early evolution of each system until the final classification of the protostellar disc can be determined. Of our 105 models, 41 fragment, 21 form a spiral structure but do not fragment, and another 12 form smooth discs. Fragmentation is more likely to occur for faster initial rotation rates and weaker magnetic fields. For stronger magnetic field strengths, the inclusion of non-ideal MHD promotes disc formation, and several of these models fragment, whereas their ideal MHD counterparts do not. For the models that fragment, there is no correlation between our parameters and where or when the fragmentation occurs. Bipolar outflows are launched in only 17 models, and these models have strong magnetic fields that are initially parallel to the rotation axis. Counter-rotating envelopes form in four slowly-rotating, strong-field models -- including one ideal MHD model -- indicating they form only in a small fraction of the parameter space investigated.

Tables

Table 1. Table 1: Summary of the results at the end of the hydrodynamics and ideal MHD models. The columns are as defined throughout the text and summarised at the beginning of this appendix.

Name	Classification	$t_{\text{disc}}$	$t_{\text{end}}-t_{\text{disc}}$	$d_{\text{frag}}$	$R_{\text{B,disc}}$	$R_{\text{T,disc}}$	$M_{\text{B,disc}}$	$M_{\text{T,disc}}$	$M_{\text{star}}$	Outflow	CRE
		[kyr]	[kyr]	[au]	[au]	[au]	[M⊙]	[M⊙]	[M⊙]
H${\Omega }_5$	Transient	116	0	-	-	-	-	-	-	no	no
H${\Omega }_{25}$	Spiral Arms	119	17.2	-	54.0	75.9	0.138	0.152	0.147	no	no
H${\Omega }_{45}$	Fragmented	125	15.7	137	47.0	149	0.0836	0.112	0.0854	no	no
H${\Omega }_{65}$	Fragmented	136	2.51	50.0	14.5	67.8	0.0356	0.0496	0.00208	no	no
H${\Omega }_{85}$	Fragmented	154	3.22	83.9	17.5	92.4	0.0329	0.0463	0.00362	no	no
I${\Omega }_5{\mu }_3{B}_{\text{-z}}$	Transient	123	0	-	-	-	-	-	-	no	no
I${\Omega }_5{\mu }_5{B}_{\text{-z}}$	Transient	119	0.716	-	-	-	-	-	0.0273	no	no
I${\Omega }_5{\mu }_{10}{B}_{\text{-z}}$	Transient	117	0.358	-	-	-	-	-	0.0226	no	no
I${\Omega }_5{\mu }_{20}{B}_{\text{-z}}$	Transient	117	0.358	-	-	-	-	-	0.0251	no	no
I${\Omega }_{25}{\mu }_3{B}_{\text{-z}}$	Transient	126	0	-	-	-	-	-	-	slow	no
I${\Omega }_{25}{\mu }_5{B}_{\text{-z}}$	Transient	121	5.37	-	-	-	-	-	0.112	slow	no
I${\Omega }_{25}{\mu }_{10}{B}_{\text{-z}}$	Smooth	120	15.8	-	42.0	60.7	0.113	0.122	0.151	no	no
I${\Omega }_{25}{\mu }_{20}{B}_{\text{-z}}$	Spiral Arms	119	16.1	-	40.0	77.3	0.104	0.120	0.153	no	no
I${\Omega }_{45}{\mu }_3{B}_{\text{-z}}$	Transient	131	3.94	-	-	-	-	-	0.0816	slow	no
I${\Omega }_{45}{\mu }_5{B}_{\text{-z}}$	Smooth	127	15.8	-	32.0	57.4	0.0853	0.0949	0.118	slow	no
I${\Omega }_{45}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	125	5.50	51.0	35.0	98.1	0.0498	0.0714	0.0484	no	no
I${\Omega }_{45}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	125	8.80	87.3	43.5	129	0.0693	0.0877	0.0632	no	no
I${\Omega }_{65}{\mu }_3{B}_{\text{-z}}$	Transient	143	4.30	-	-	-	-	-	0.0815	slow	no
I${\Omega }_{65}{\mu }_5{B}_{\text{-z}}$	Spiral Arms	139	15.8	-	36.5	70.4	0.0516	0.0663	0.0932	slow	no
I${\Omega }_{65}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	136	2.86	57.0	28.0	76.6	0.0234	0.0406	0.0244	no	no
I${\Omega }_{65}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	136	3.58	88.0	33.0	102	0.0300	0.0557	0.0280	no	no
I${\Omega }_{85}{\mu }_3{B}_{\text{-z}}$	Smooth	162	15.8	-	12.0	15.4	0.00356	0.00388	0.137	slow	no
I${\Omega }_{85}{\mu }_5{B}_{\text{-z}}$	Fragmented	158	3.22	34.0	35.5	51.3	0.0369	0.0408	0.0236	slow	no
I${\Omega }_{85}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	156	1.79	58.2	25.5	75.6	0.0178	0.037	0.0172	no	no
I${\Omega }_{85}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	155	2.51	72.8	24.5	91.7	0.0181	0.0355	0.0162	no	no
I${\Omega }_5{\mu }_3{B}_{\text{-x}}$	Transient	122	0.716	-	-	-	-	-	0.0404	no	yes
I${\Omega }_5{\mu }_5{B}_{\text{-x}}$	Transient	119	0.358	-	-	-	-	-	0.0264	no	no
I${\Omega }_5{\mu }_{10}{B}_{\text{-x}}$	Transient	117	0.358	-	-	-	-	-	0.0236	no	no
I${\Omega }_5{\mu }_{20}{B}_{\text{-x}}$	Transient	117	0.358	-	-	-	-	-	0.0238	no	no
I${\Omega }_{25}{\mu }_3{B}_{\text{-x}}$	Transient	125	0.358	-	-	-	-	-	0.0240	no	no
I${\Omega }_{25}{\mu }_5{B}_{\text{-x}}$	Transient	121	1.07	-	-	-	-	-	0.0419	no	no
I${\Omega }_{25}{\mu }_{10}{B}_{\text{-x}}$	Smooth	120	15.8	-	32.5	43.0	0.101	0.105	0.193	no	no
I${\Omega }_{25}{\mu }_{20}{B}_{\text{-x}}$	Spiral Arms	119	15.8	-	43.5	79.7	0.103	0.116	0.165	no	no
I${\Omega }_{45}{\mu }_3{B}_{\text{-x}}$	Transient	132	0.716	-	-	-	-	-	0.0274	no	no
I${\Omega }_{45}{\mu }_5{B}_{\text{-x}}$	Smooth	128	15.8	-	34.5	41.0	0.0915	0.0938	0.168	no	no
I${\Omega }_{45}{\mu }_{10}{B}_{\text{-x}}$	Fragmented	126	12.1	150	38.5	158	0.0848	0.0983	0.107	no	no
I${\Omega }_{45}{\mu }_{20}{B}_{\text{-x}}$	Fragmented	125	11.0	70.0	39.5	200	0.0695	0.0899	0.0797	no	no
I${\Omega }_{65}{\mu }_3{B}_{\text{-x}}$	Transient	145	0.716	-	-	-	-	-	0.0237	no	no
I${\Omega }_{65}{\mu }_5{B}_{\text{-x}}$	Spiral Arms	140	15.8	-	39.5	65.4	0.0737	0.0881	0.127	no	no
I${\Omega }_{65}{\mu }_{10}{B}_{\text{-x}}$	Fragmented	137	11.8	79.6	41.0	110	0.0625	0.0868	0.0783	no	no
I${\Omega }_{65}{\mu }_{20}{B}_{\text{-x}}$	Fragmented	136	6.80	36.0	38.5	111	0.0522	0.0578	0.0433	no	no
I${\Omega }_{85}{\mu }_3{B}_{\text{-x}}$	Transient	168	0.358	-	-	-	-	-	0.0224	no	no
I${\Omega }_{85}{\mu }_5{B}_{\text{-x}}$	Spiral Arms	162	15.8	-	44.5	70.0	0.0685	0.0788	0.114	no	no
I${\Omega }_{85}{\mu }_{10}{B}_{\text{-x}}$	Fragmented	158	14.3	71.0	45.0	109	0.0579	0.0729	0.0729	no	no
I${\Omega }_{85}{\mu }_{20}{B}_{\text{-x}}$	Fragmented	155	4.30	55.1	34.5	67.1	0.0314	0.0449	0.0259	no	no

Table 2. Table 2: Summary of the results at the end of the non-ideal MHD models. The columns are as defined throughout the text and summarised at the beginning of this appendix.

Name	Classification	$t_{\text{disc}}$	$t_{\text{end}}-t_{\text{disc}}$	$d_{\text{frag}}$	$R_{\text{B,disc}}$	$R_{\text{T,disc}}$	$M_{\text{B,disc}}$	$M_{\text{T,disc}}$	$M_{\text{star}}$	Outflow	CRE
		[kyr]	[kyr]	[au]	[au]	[au]	[M⊙]	[M⊙]	[M⊙]
N${\Omega }_5{\mu }_3{B}_{\text{-z}}$	Transient	122	0.716	-	-	-	-	-	0.0307	no	yes
N${\Omega }_5{\mu }_5{B}_{\text{-z}}$	Transient	119	0	-	-	-	-	-	-	no	yes
N${\Omega }_5{\mu }_{10}{B}_{\text{-z}}$	Transient	117	0.358	-	-	-	-	-	0.0237	no	no
N${\Omega }_5{\mu }_{20}{B}_{\text{-z}}$	Transient	116	0.716	-	-	-	-	-	0.0254	no	no
N${\Omega }_{25}{\mu }_3{B}_{\text{-z}}$	Smooth	125	15.4	-	32.0	39.4	0.107	0.108	0.190	slow	no
N${\Omega }_{25}{\mu }_5{B}_{\text{-z}}$	Spiral Arms	122	15.8	-	41.0	75.8	0.104	0.117	0.173	no	no
N${\Omega }_{25}{\mu }_{10}{B}_{\text{-z}}$	Spiral Arms	120	15.8	-	50.5	71.2	0.132	0.148	0.154	no	no
N${\Omega }_{25}{\mu }_{20}{B}_{\text{-z}}$	Spiral Arms	119	15.8	-	41.5	79.4	0.0974	0.124	0.145	no	no
N${\Omega }_{45}{\mu }_3{B}_{\text{-z}}$	Smooth*	131	15.7	-	37.0	57.0	0.0923	0.0994	0.145	fast	no
N${\Omega }_{45}{\mu }_5{B}_{\text{-z}}$	Fragmented	127	10.2	73.8	46.0	113	0.0864	0.103	0.0777	no	no
N${\Omega }_{45}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	125	8.25	85.2	36.5	111	0.0547	0.0861	0.0581	no	no
N${\Omega }_{45}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	125	8.80	41.2	45.0	103	0.0713	0.102	0.0587	no	no
N${\Omega }_{65}{\mu }_3{B}_{\text{-z}}$	Spiral Arms	142	14.9	-	33.5	53.1	0.0788	0.0850	0.115	fast	no
N${\Omega }_{65}{\mu }_5{B}_{\text{-z}}$	Fragmented	138	3.27	57.0	25.0	96.2	0.0354	0.0525	0.0236	no	no
N${\Omega }_{65}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	136	3.16	76.4	28.5	94.8	0.0522	0.0692	0.0355	no	no
N${\Omega }_{65}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	135	2.80	49.5	14.5	62.7	0.0357	0.0569	2.68$\times {10}^{-4}$	no	no
N${\Omega }_{85}{\mu }_3{B}_{\text{-z}}$	Fragmented	162	3.07	64.9	27.0	74.7	0.0353	0.0466	0.0207	slow	no
N${\Omega }_{85}{\mu }_5{B}_{\text{-z}}$	Fragmented	168	2.96	68.7	23.5	79.5	0.0446	0.0628	3.32$\times {10}^{-3}$	no	no
N${\Omega }_{85}{\mu }_{10}{B}_{\text{-z}}$	Fragmented	155	2.86	48.8	23.5	84.8	0.0381	0.0564	2.04$\times {10}^{-3}$	no	no
N${\Omega }_{85}{\mu }_{20}{B}_{\text{-z}}$	Fragmented	154	3.16	52.3	25.5	91.9	0.0349	0.0478	2.97$\times {10}^{-3}$	no	no
N${\Omega }_5{\mu }_3{B}_{\text{+z}}$	Transient	122	0.716	-	-	-	-	-	0.0293	fast	no
N${\Omega }_5{\mu }_5{B}_{\text{+z}}$	Transient	119	0.716	-	-	-	-	-	0.0291	no	no
N${\Omega }_5{\mu }_{10}{B}_{\text{+z}}$	Transient	117	0.358	-	-	-	-	-	0.0231	no	no
N${\Omega }_5{\mu }_{20}{B}_{\text{+z}}$	Transient	116	0.716	-	-	-	-	-	0.0261	no	no
N${\Omega }_{25}{\mu }_3{B}_{\text{+z}}$	Transient	125	0	-	-	-	-	-	0	fast	no
N${\Omega }_{25}{\mu }_5{B}_{\text{+z}}$	Smooth*	122	15.4	-	39.0	57.7	0.107	0.118	0.178	no	no
N${\Omega }_{25}{\mu }_{10}{B}_{\text{+z}}$	Fragmented	120	11.8	83.2	38.5	108	0.100	0.111	0.127	no	no
N${\Omega }_{25}{\mu }_{20}{B}_{\text{+z}}$	Spiral Arms	119	15.8	-	42.5	81.4	0.107	0.125	0.142	no	no
N${\Omega }_{45}{\mu }_3{B}_{\text{+z}}$	Smooth	131	15.3	-	29.0	31.8	0.0776	0.0780	0.174	fast	no
N${\Omega }_{45}{\mu }_5{B}_{\text{+z}}$	Fragmented	127	7.94	69.2	36.0	132	0.0628	0.0841	0.0679	no	no
N${\Omega }_{45}{\mu }_{10}{B}_{\text{+z}}$	Fragmented	125	6.09	108	35.0	118	0.0535	0.0721	0.0445	no	no
N${\Omega }_{45}{\mu }_{20}{B}_{\text{+z}}$	Fragmented	125	8.97	71.9	36.0	104	0.0556	0.0896	0.0615	no	no
N${\Omega }_{65}{\mu }_3{B}_{\text{+z}}$	Smooth*	142	15.2	-	31.5	37.9	0.0878	0.0892	0.132	fast	no
N${\Omega }_{65}{\mu }_5{B}_{\text{+z}}$	Fragmented	138	2.86	72.4	25.5	81.1	0.0456	0.0601	8.52$\times {10}^{-3}$	no	no
N${\Omega }_{65}{\mu }_{10}{B}_{\text{+z}}$	Fragmented	136	2.73	55.0	15.0	70.9	0.0361	0.0568	7.78$\times {10}^{-4}$	no	no
N${\Omega }_{65}{\mu }_{20}{B}_{\text{+z}}$	Fragmented	135	2.84	47.4	15.0	73.4	0.0362	0.0583	1.82$\times {10}^{-4}$	no	no
N${\Omega }_{85}{\mu }_3{B}_{\text{+z}}$	Fragmented	162	3.58	51.8	26.0	81.8	0.0355	0.0456	0.0259	slow	no
N${\Omega }_{85}{\mu }_5{B}_{\text{+z}}$	Fragmented	158	2.73	62.6	21.5	77.2	0.0410	0.0607	1.60$\times {10}^{-3}$	no	no
N${\Omega }_{85}{\mu }_{10}{B}_{\text{+z}}$	Fragmented	155	2.80	37.3	15.5	74.0	0.0311	0.0511	2.75$\times {10}^{-4}$	no	no
N${\Omega }_{85}{\mu }_{20}{B}_{\text{+z}}$	Fragmented	154	3.49	55.7	25.5	84.4	0.0372	0.0589	3.08$\times {10}^{-3}$	no	no
N${\Omega }_5{\mu }_3{B}_{\text{-x}}$	Transient	122	0.358	-	-	-	-	-	0.0214	no	yes
N${\Omega }_5{\mu }_5{B}_{\text{-x}}$	Transient	119	0.716	-	-	-	-	-	0.0287	no	no
N${\Omega }_5{\mu }_{10}{B}_{\text{-x}}$	Transient	117	0.358	-	-	-	-	-	0.0217	no	no
N${\Omega }_5{\mu }_{20}{B}_{\text{-x}}$	Transient	117	0.358	-	-	-	-	-	0.0241	no	no
N${\Omega }_{25}{\mu }_3{B}_{\text{-x}}$	Smooth	126	15.8	-	31.5	36.5	0.102	0.103	0.211	no	no
N${\Omega }_{25}{\mu }_5{B}_{\text{-x}}$	Smooth	122	15.8	-	42.0	52.2	0.124	0.128	0.176	no	no
N${\Omega }_{25}{\mu }_{10}{B}_{\text{-x}}$	Spiral Arms	120	15.8	-	50.0	83.4	0.112	0.132	0.162	no	no
N${\Omega }_{25}{\mu }_{20}{B}_{\text{-x}}$	Spiral Arms	119	15.9	-	47.0	98.4	0.112	0.138	0.144	no	no
N${\Omega }_{45}{\mu }_3{B}_{\text{-x}}$	Spiral Arms	132	16.0	-	46.0	59.5	0.108	0.117	0.160	no	no
N${\Omega }_{45}{\mu }_5{B}_{\text{-x}}$	Spiral Arms	128	15.8	-	41.0	84.5	0.0796	0.113	0.126	no	no
N${\Omega }_{45}{\mu }_{10}{B}_{\text{-x}}$	Spiral Arms	125	15.6	-	41.0	83.8	0.0803	0.101	0.0982	no	no
N${\Omega }_{45}{\mu }_{20}{B}_{\text{-x}}$	Spiral Arms	125	16.0	-	48.5	112	0.0867	0.118	0.0887	no	no
N${\Omega }_{65}{\mu }_3{B}_{\text{-x}}$	Spiral Arms	145	15.4	-	39.5	67.9	0.0749	0.102	0.134	no	no
N${\Omega }_{65}{\mu }_5{B}_{\text{-x}}$	Spiral Arms	140	15.6	-	41.5	96.8	0.0651	0.0954	0.0982	no	no
N${\Omega }_{65}{\mu }_{10}{B}_{\text{-x}}$	Fragmented	136	11.9	105	40.0	143	0.0605	0.0822	0.0611	no	no
N${\Omega }_{65}{\mu }_{20}{B}_{\text{-x}}$	Fragmented	135	3.44	41.3	28.0	75.3	0.0443	0.0542	0.0899	no	no
N${\Omega }_{85}{\mu }_3{B}_{\text{-x}}$	Spiral Arms	168	15.5	-	53.0	62.9	0.102	0.110	0.131	no	no
N${\Omega }_{85}{\mu }_5{B}_{\text{-x}}$	Spiral Arms	162	15.4	-	42.5	99.3	0.0545	0.0855	0.0900	no	no
N${\Omega }_{85}{\mu }_{10}{B}_{\text{-x}}$	Fragmented	157	11.6	93.0	34.5	129	0.0325	0.0646	0.0526	no	no
N${\Omega }_{85}{\mu }_{20}{B}_{\text{-x}}$	Fragmented	155	3.47	32.5	23.5	79.6	0.0396	0.0534	0.0591	no	no

Equations

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

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

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

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

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

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

$${\frac{\text{d}\mathbf{B}}{\text{d}t}|}_{\text{non-ideal}}=$$

$${\frac{\text{d}u}{\text{d}t}|}_{\text{non-ideal}}$$

$$\mu \equiv \frac{M/{\Phi }_{\text{B}}}{{\left(M/{\Phi }_{\text{B}}\right)}_{\text{crit}}},$$

$$Q=\frac{\kappa c_{\text{s}}}{\pi \Sigma G}$$

$$Q_{\text{m}}=\frac{\kappa \sqrt{c_{\text{s}}^2+v_{\text{A}}^2}}{\pi \Sigma G},$$

$${\kappa }^2=\frac{2\Omega }{r}\frac{\text{d}}{\text{d}r}\left(r^2\Omega \right),$$

$$\frac{H}{R}\lesssim \frac{M_{\text{disc}}}{M_{\text{star}}},$$

$$M_{\text{T}}=\frac{\pi c_{\text{s}}^4}{G^2\Sigma },$$

$${\Sigma }_{\text{max}}=\frac{\pi c_{\text{s}}^4}{G^2{m}_{\text{p}}{N}_{\text{reso}}},$$

$${\Sigma }_{\text{m,max}}=\frac{\pi {\left(c_{\text{s}}^2+v_{\text{A}}^2\right)}^2}{G^2{m}_{\text{p}}{N}_{\text{reso}}}$$

$$H=\frac{\sqrt{c_{\text{s}}^2+v_{\text{A}}^2}}{\Omega }.$$

Full text

††pagerange: Disc formation and fragmentation using radiative non-ideal magnetohydrodynamics–Disc formation and fragmentation using radiative non-ideal magnetohydrodynamics††pubyear: 2016

1 Introduction

Stars are formed from the gravitational collapse of gas clouds, which are observed to have strong magnetic fields (e.g. Heiles & Crutcher, 2005) and low ionisation rates (e.g. Mestel & Spitzer, 1956; Nakano & Umebayashi, 1986; Umebayashi & Nakano, 1990). Dense molecular cloud cores within these clouds are initially rotating with rates of $\beta_{\text{r}}\lesssim 0.15$ with a mean value of $\beta_{\text{r}}\sim 0.02$ (Goodman et al., 1993), where $\beta_{\text{r}}$ is the ratio of rotational to gravitational energy. Since the dense core is initially rotating, a self-gravitating disc is expected to form during the formation of the protostar (e.g. Terebey et al., 1984; Attwood et al., 2009; Bate, 2011; Machida & Matsumoto, 2011), and large gas discs have been inferred to exist around even young (Class 0) protostars (e.g. Dunham et al., 2011; Lindberg et al., 2014; Tobin et al., 2015, 2016).

During the gravitational collapse, a massive first hydrostatic core can become bar unstable to form a gravitationally unstable disc and develop spiral arms (e.g. Bate, 1998; Saigo & Tomisaka, 2006; Saigo et al., 2008). Alternatively, if a self-gravitating disc forms, it may undergo a gravitational instability due to growing non-axisymmetric modes and then form spiral arms (e.g. Papaloizou & Savonije, 1991; Durisen et al., 2007). The disc itself, or more typically the spiral arms, may further be unstable to fragmentation (e.g. Bonnell, 1994; Bonnell & Bate, 1994). There are several criteria that are used to predict the stability of a disc, including the Toomre-Q parameter (Toomre, 1964), the cooling parameter (Gammie, 2001), and the ratio of disc-to-stellar masses (Gammie, 2001).

The study of disc fragmentation is important for better understanding the formation of multiple systems (e.g. Kratter & Matzner, 2006; Nayakshin et al., 2007; Kratter et al., 2010), the formation of brown dwarfs or planets (e.g. Boss, 1997, 1998, 2001; Mayer et al., 2002; Stamatellos et al., 2007), and episodic accretion events (e.g. Vorobyov & Basu, 2005, 2006, 2015). These studies have been performed starting from both gas clouds (e.g. Stamatellos et al., 2011; Forgan & Rice, 2012), and from Keplerian discs (e.g. Mayer et al., 2002; Rice et al., 2003; Stamatellos & Whitworth, 2009; Meru & Bate, 2010, 2012; Vorobyov, 2013; Meru, 2015; Forgan et al., 2017; Mercer & Stamatellos, 2017; Hall et al., 2017). While many studies stop once the fragmentation limit is determined from their initial conditions and physical processes, several studies (e.g. Vorobyov & Basu, 2015; Hall et al., 2017; Forgan et al., 2017) continue to evolve the system to further study the disc and the fragments.

Studies starting from a Keplerian disc typically neglect magnetic fields since discs are expected to be weakly ionised (i.e. with an ionisation fraction of $10^{-12}$; Fromang et al., 2002 and references therein), and contain magnetic dead zones in which the magnetic Reynolds number is lower than some critical value (e.g. Gammie, 1996). However, Forgan et al. (2017) modelled a disc in the presence of ideal magnetic fields and found that the fragments in the magnetised disc were larger and formed at a smaller range of radii than in their purely hydrodynamical counterpart, and that these properties were dependent on the initial magnetic field strength. They concluded that magnetic fields influence the fragmentation limit.

Studies that start from a gas cloud and include magnetic fields to match observations have shown that the resulting disc properties are at least partly dependent on the initial magnetic field strength and geometry (e.g. Allen et al., 2003; Price & Bate, 2007; Hennebelle & Fromang, 2008; Duffin & Pudritz, 2009; Hennebelle & Ciardi, 2009; Commerçon et al., 2010; Seifried et al., 2011; Lewis et al., 2015). These early studies starting from a gas cloud used ideal magnetohydrodynamics (MHD), and suffered from the ‘magnetic braking catastrophe’ (Allen et al., 2003) since the magnetic fields so efficiently extracted angular momentum that large discs failed to form in the presence of strong magnetic fields.

Detailed models of molecular clouds found have ionisation fractions as low as $10^{-14}$ (Nakano & Umebayashi, 1986; Umebayashi & Nakano, 1990), suggesting ideal MHD is a poor approximation. Rather than neglecting magnetic fields and incorrectly assuming a purely hydrodynamic collapse, simulations began including non-ideal MHD to account for ionised and neutral species (e.g. Ciolek & Mouschovias, 1994; Li & Shu, 1996; Mellon & Li, 2009; Duffin & Pudritz, 2009; Dapp & Basu, 2010; Machida et al., 2011; Li et al., 2011; Tomida et al., 2013; Tomida et al., 2015; Tsukamoto et al., 2015b, a; Wurster et al., 2016; Tsukamoto et al., 2017; Wurster et al., 2018a, c). The more recent simulations that included the Hall effect were able to overcome the magnetic braking catastrophe and produce discs comparable with observed discs. These non-ideal MHD studies typically focused on disc formation rather than fragmentation, thus there has yet to be a fragmentation study using non-ideal MHD starting from a molecular cloud core.

In this study, we model the formation and fragmentation of discs using a 3D self-gravitating, smoothed particle, radiative, non-ideal magnetohydrodynamics code. We self-consistently form a disc by allowing a low-mass molecular cloud to collapse, and we follow the evolution until a disc forms and it dissipates, is proved to be stable, or it fragments. Thus, we primarily aim to determine the initial conditions that will lead to a fragmenting disc; we do not follow the evolution of the fragments. In Section 2 we describe our methods and in Section 3 we give our initial conditions. Results are presented and discussed in Sections 4 and 5, respectively; given the large suite, the discussion also includes aspects that are not directly related to fragmentation, but provide additional insights into the effect of varying our free parameters. We conclude in Section 6.

2 Methods

2.1 Radiative non-ideal magnetohydrodynamics

We solve the set of radiative non-ideal MHD equations, given by

[TABLE]

where $\text{d}/\text{d}t\equiv\partial/\partial t+\bm{v}\cdot\bm{\nabla}$ is the Lagrangian derivative, $\rho$ is the density, $\bm{v}$ is the velocity, $P$ the hydrodynamic pressure, $\bm{B}$ is the magnetic field, $\Phi$ is the gravitational potential, $B_{\text{P}}$ is the frequency-integrated Plank function, $E$ is the radiation energy density, $\bm{F}$ is the radiative flux, $\bm{P}$ is the radiation pressure tensor, $\kappa$ is the opacity, $u$ is the specific energy of the gas, $\mathbb{I}$ is the identity matrix, $c$ is the speed of light and $G$ is the gravitational constant; the magnetic field has been normalised such that the Alfvén velocity is defined as $v_{\text{A}}\equiv B/\sqrt{\rho}$ in code units (see Price & Monaghan, 2004). The radiative transfer algorithm is given in Whitehouse et al. (2005) and Whitehouse & Bate (2006) and uses a two-temperature (matter and radiation) flux-limited diffusion approximation and assumes local thermodynamic equilibrium; the opacity is assumed to be independent of frequency and there is no distinction between absorption and total opacities.

The contribution of the non-ideal MHD processes to the magnetic field and the internal energy are

[TABLE]

and

[TABLE]

respectively, where $\bm{J}=\bm{\nabla}\times\bm{B}$ is the current density, and $\eta_{\text{OR}}$, $\eta_{\text{HE}}$ and $\eta_{\text{AD}}$ are the non-ideal MHD coefficients for Ohmic resistivity, the Hall effect and ambipolar diffusion, respectively.

We use version 1.2.1 of the Nicil library (Wurster, 2016) to calculate the non-ideal coefficients, $\eta$. At low temperatures, cosmic rays ionise a heavy ion and a light ion at the canonical rate of $\zeta_{\text{cr}}=10^{-17}$ s*-1* (Spitzer & Tomasko, 1968; Umebayashi & Nakano, 1981). A single dust grain population is modelled as three species: a positively (negatively) charged grain species that has lost (absorbed) an electron and a neutral species. The grains have a radius and bulk density of $a_{\text{g}}=0.1\mu$m and $\rho_{\text{b}}=3$ g cm*-3*, respectively (Pollack et al., 1994).

2.2 Numerical methods

To perform our simulations, we use the 3D smoothed particle hydrodynamics (SPH) code sphNG with the inclusion of self-gravity, radiative hydrodynamics and non-ideal MHD; this code is based upon the original version by Benz (1990) and Benz et al. (1990), but has since been substantially modified by Bate et al. (1995) and many additional contributors.

For a review of the discretised MHD equations, see Price (2012). Briefly, we adopt the usual cubic spline kernel, set such that the smoothing length is given by $h=1.2\left(m/\rho\right)^{1/3}$, where $m$ is the mass of an SPH particle; this yields $N_{\text{neigh}}\sim 58$ neighbours in three dimensions. We calculate the gravitational forces following Price & Monaghan (2007) at short range, and use a binary tree to compute the long range gravitational interactions. For magnetic stability, the Børve et al. (2001) source-term approach is used, and artificial resistivity is included to capture the magnetic discontinuities (Price & Monaghan, 2005; Price, 2012); the artificial resistivity parameter is given by $\alpha_{\text{B}}=\max\left(h\left|\nabla\bm{B}\right|/\left|\bm{B}\right|,1\right)$ (Tricco & Price, 2013). We employ the constrained hyperbolic divergence cleaning algorithm of Tricco & Price (2012) and Tricco et al. (2016) to control divergence errors in the magnetic field.

To model radiation transport, we use the flux-limited diffusion method described in Bate et al. (2014), where the method is described in detail in Whitehouse et al. (2005) and Whitehouse & Bate (2006). We use 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. At low temperatures, the mean molecular weight is taken to be $\mu=2.38$, and we use opacity tables from Pollack et al. (1985) and Ferguson et al. (2005).

The code is parallelised using both OpenMP and MPI. It does not include super-timestepping (Alexiades et al., 1996) as used in the non-ideal MHD studies by Wurster et al. (2016, 2017b, 2018b), nor implicit timestepping for Ohmic resistivity as introduced in Wurster et al. (2018a).

3 Initial conditions

We begin with a Bonnor-Ebert sphere (Bonnor, 1956; Ebert, 1955) of radius $R=1.3\times~{}10^{17}$ cm, mass $M=1$ M⊙, temperature $T=8$K, and concentration parameter $\xi=7.45$; this corresponds to a density ratio of 20:1 between the inner and outer regions of the sphere. The ratio of thermal energy to gravitational potential energy is $\alpha=0.50$. The sphere is given an initial solid-body rotation $\bm{\Omega}_{0}=\Omega_{0}\hat{\bm{z}}$, and threaded with a uniform magnetic field $\bm{B}_{0}$.

The sphere is placed in a low-density box of edge length $l=4R$ at a density ratio of 1:30 with the edge of the sphere; the sphere and low-density medium are in pressure equilibrium which prevents the sphere from artificially expanding into the low-density medium. This two-medium set-up allows us to place boundary conditions at the edge of the box rather than the edge of the sphere. This is especially useful for the magnetic field, which is uniform initially and, therefore, we can use periodic boundary conditions at the edges of the box. We use quasi-periodic boundary conditions at the edge of the box, in which SPH particles interact hydrodynamically ‘across the box’, but not gravitationally. The box size was chosen to prevent any boundary effects from influencing the evolution of the sphere.

Sink particles (Bate et al., 1995) of radius 1 au are unconditionally inserted when the maximum density reaches $\rho_{\text{crit}}=5\times 10^{-10}$ g cm*-3*, and we permit only one sink particle to form per simulation. Given the small timesteps required to evolve non-ideal MHD at high densities, the introduction of sink particles is necessary to follow the evolution of the disc long enough to determine its stability. However, we note that sink particles stabilise small discs against instabilities (Machida et al., 2014), thus the final classification of models with small discs may be influenced by the sink particle.

Our simulations include $10^{6}$ SPH particles in the sphere and an additional $1.8\times 10^{5}$ particles in the low-density medium. Resolving the Jeans length throughout the collapse requires at least $3\times 10^{4}$ particles per solar mass (Bate & Burkert, 1997), thus the Jeans mass is well resolved at all times. The equal-mass particles in the sphere are initially placed on a regular close-packed lattice, which is then deformed to produce the Bonnor-Ebert sphere described above; the SPH particles in the warm medium have the same mass as the particles in the sphere and are also placed on a regular close-packed lattice.

We characterise the initial rotation by the orbital rotation at the outer radius of the cloud, $\Omega_{\text{orbit}}=2.45\times 10^{-13}$ rad s*-1*. The five rotation speeds we test are $\Omega_{0}=$ 0.05, 0.25, 0.45, 0.65 and 0.85 $\Omega_{\text{orbit}}$, which correspond to ratios of rotational energy to gravitational potential energy of $\beta_{\text{r}}=4.4\times 10^{-4}$, 0.011, 0.035, 0.074 and 0.13. We are thus primarily exploring the higher end of the distribution of $\beta_{\text{r}}$-values (Goodman et al., 1993), but the faster rotators are more likely to fragment than slower rotators.

Our suite of models consists of purely hydrodynamical models, ideal and non-ideal MHD models. For the magnetised models, we characterise the magnetic field in terms of the normalised mass-to-flux parameter

[TABLE]

where $M/\Phi_{\text{B}}\equiv M/\left(\pi R^{2}B\right)$ is the mass-to-flux ratio and $\left(M/\Phi_{\text{B}}\right)_{\text{crit}}=c_{1}/(3\pi)\sqrt{5/G}$ is the critical value in CGS units where magnetic fields prevent gravitational collapse altogether; here, $\Phi_{\text{B}}$ is the magnetic flux threading the surface of the (spherical) cloud at radius $R$ assuming a uniform magnetic field of strength $B$, and $c_{1}\simeq 0.53$ is a parameter numerically determined by Mouschovias & Spitzer (1976). In this study, we test initial values of $\mu_{0}=$ 3, 5, 10 and 20, which correspond to magnetic field strengths of $B_{0}=$ 25.6, 15.4, 7.69 and 3.85 $\mu$G, respectively.

In our ideal MHD models, we test two magnetic field orientations: $\bm{B}_{0}=\bm{B}_{\text{-z}}\equiv-B_{0}\hat{\bm{z}}$ (i.e. the magnetic field is initially parallel to the rotation axis), and $\bm{B}_{0}=\bm{B}_{\text{-x}}\equiv-B_{0}\hat{\bm{x}}$ (i.e. the magnetic field is initially perpendicular to the rotation axis). In non-ideal MHD, the Hall effect is dependent on $\bm{\Omega}\cdot\bm{B}$ (Braiding & Wardle, 2012), thus we test $\bm{B}_{0}=\bm{B}_{\text{+z}}\equiv B_{0}\hat{\bm{z}}$ in addition to $\bm{B}_{\text{-z}}$ and $\bm{B}_{\text{-x}}$. We only test one orientation perpendicular to the axis of rotation (i.e. $\bm{B}_{\text{-x}}$) since $\bm{\Omega}\cdot\bm{B}_{\pm\text{x}}=0$, suggesting similar results should be obtained for $\bm{B}_{0}=\bm{B}_{\pm\text{x}}$. However, the binary formation study of Wurster et al. (2017b) found slightly different results due to the structures that formed as the systems evolved.

Due to the large suite and computational limitations, the simulations are run until the disc classification can be determined or until $t_{\text{end}}-t_{\text{disc}}\approx 16$ kyr (i.e. $\approx$16 kyr after formation of the disc). Models with small, dense, strongly magnetised discs typically have the shortest end-time relative to the disc formation time (aside from those models whose discs dissipate) due to the short timestep required to resolve the processes; this timestep is decreased even more in the models with the non-ideal MHD processes since they require an even shorter timestep (e.g. Mac Low et al., 1995; Choi et al., 2009; Wurster et al., 2014).

4 Results

In this section, we present the results of our suite of 105 models, and in Section 5 we discuss their implications. Details and properties of the outcome of each model (i.e., classification, disc formation time, simulation end-time, disc radius, disc mass, stellar mass, outflow and envelope properties) are listed in Tables  1 and 2 in Appendix A.

Our magnetised models use the naming convention of a$\Omega_{b}\mu_{c}B_{\text{\emph{d}}}$, where $a=$ I (N) for ideal (non-ideal) MHD, $b$ represents 100$\times$ the initial angular rotation in terms of $\Omega_{\text{orbit}}$, $c$ represent the initial mass-to-flux ratio in units of the critical mass-to-flux ratio $\mu_{0}$, and $d$ represents the orientation of the initial magnetic field ($\pm z$ or $-x$); our hydrodynamic models are named H$\Omega_{b}$. An asterisk, *, in place of a variable indicates every model with the remaining defined components.

4.1 Identifying and classifying discs

We define the total disc radius, $R_{\text{T,disc}}$, as the radius which includes all the gas that satisfies $\rho>\rho_{\text{thresh}}$, where $\rho_{\text{thresh}}=10^{-13}$ g cm*-3* which is approximately the density at which the collapsing gas becomes adiabatic; this radius will include the spiral arms and the gaps between them. The total mass of the disc, $M_{\text{T,disc}}$, includes all the gas with $\rho>\rho_{\text{thresh}}$.

To define the ‘bulk’ disc, we divide the total disc into tori of width 0.5 au and set the height to include all the gas with $\rho>\rho_{\text{thresh}}$. The bulk radius, $R_{\text{B,disc}}$, is then defined as the outer extent of the outermost tori where 80 per cent of the gas particles in the torus have $\rho>\rho_{\text{thresh}}$. The mass of the bulk disc is the total mass of the gas with $\rho>\rho_{\text{thresh}}$ within the bulk radius. For this disc to ‘exist’, $R_{\text{B,disc}}\geq R_{\text{crit,disc}}$ and $M_{\text{B,disc}}\geq M_{\text{crit,disc}}$, where $R_{\text{crit,disc}}=1$ au is the radius of the sink particle and $M_{\text{crit,disc}}=6\times 10^{-4}$ M⊙. This is the disc that we will typically be referencing. The formation time of the disc is the earliest time when the disc ‘exists’, as per these criteria.

We do not include any velocity or pressure criteria in this definition, thus we cannot be certain that these discs are rotationally supported. Thus, when we refer to discs throughout this study, we are actually referring to ‘disc-like structures.’

We define four classifications for our models: transient, a smooth disc, spiral arms without fragmentation, and a fragmented disc. The ‘transient’ classification is given when a disc dissipates, i.e. when $R_{\text{B,disc}}<R_{\text{crit,disc}}$ or $M_{\text{B,disc}}<M_{\text{crit,disc}}$ occurs after a disc has formed. Once a disc dissipates, we end the simulation, thus it is possible that a disc may re-form at a later time, but this is out of the scope of our study. The classifications of ‘smooth’ and ‘spiral arms’ are given to models that maintain a disc and do not fragment; this classification is given by visual inspection.

4.2 Identifying fragments

There are many ways to define a fragment, including searching for minima in the gravitational potential (e.g. Smith et al., 2009), finding density maxima $\sim$3 dex greater than their surroundings (e.g. Meru & Bate, 2010), or using density gradients (e.g. Hall et al., 2017). Hall et al. (2017) found that the gravitational potential method was more restrictive than the density gradient method and found only clumps that survived for long periods of time, whereas the density gradient method found more clumps, including those that later merged or were accreted by the star. We tested several methods of locating fragments and all methods found similar formation times for any given clump, within an uncertainty of $\pm\Delta t\approx 0.36$ kyr.

In our study, we locate fragments using a method analogous to determining when sink particles are to be inserted. We find the densest particle with $\rho>10^{-11}$ g cm*-3* that is not associated with the central core. We then calculate $\bm{\nabla}\cdot\bm{v}$ of this particle using the particles within 2 au of it, and we further calculate the total energy $E$ of these particles. If $\bm{\nabla}\cdot\bm{v}<0$ and $E<0$, then we define the clump as a fragment.

Although the fragmentation time is similar between our various methods, the fragmentation radius has larger variability. The reason is that fragments are typically formed in the outer regions of the discs where the gas has not yet formed a stable orbit. Thus the fragments and their progenitor gas radially migrate, and a slight change in our fragmentation criteria will thus have a direct impact on the fragmentation distance. Although we present fragmentation times and distances, we do so with caution, and caution against a rigorous quantitative analysis due to these large uncertainties.

An analysis of the evolution of the fragments is beyond the scope of this study.

4.3 General trends

Representative examples of the classifications are shown in Fig. 1. A graphical summary of the final classification of each model is given in Fig. 2.

Our suite contains 105 models, of which 31 do not form discs, 12 form smooth discs, 21 form spiral arms but do not fragment, and the remaining 41 fragment.

All models with $\Omega_{0}=0.05\Omega_{\text{orbit}}$ are classified as transient models. This result is independent of all magnetic properties.

In the hydrodynamics models, discs fragment for $\Omega_{0}\geq 0.45\Omega_{\text{orbit}}$. In the magnetised models, discs fragment for fast rotations (large $\Omega_{0}$) and weak magnetic field strengths (large $\mu_{0}$). There is a general transition from fragmented, to structured, to smooth, to transients as the rotation rate is decreased and/or the magnetic field strength is increased.

For ideal MHD, the final classifications are approximately independent of the initial direction of the magnetic field (i.e. $\bm{B}_{\text{-z}}$ vs $\bm{B}_{\text{-x}}$), however, the disc and fragmentation properties are dependent on the direction.

When including the non-ideal MHD processes, larger discs tend to form than in their ideal MHD counterparts. For the $\bm{B}_{\text{-x}}$ models, four models fail to form a disc in non-ideal MHD compared to nine in ideal MHD; however, fewer of the non-ideal MHD discs fragment (four compared to six). For $\bm{B}_{\text{-z}}$, more discs form when using non-ideal MHD, and these discs are more likely to fragment than their ideal MHD counterparts.

For $\Omega_{0}\geq 0.45\Omega_{\text{orbit}}$, the classifications are the same for the non-ideal MHD models with $\bm{B}_{\pm\text{z}}$; the classifications are also the same for the $\bm{B}_{\pm\text{z}}$ models with $\mu_{0}=20$. This implies that the Hall effect is not efficient enough in these regimes to affect the classification; in both regimes, the angular momentum contribution from the Hall effect is small compared to the initial angular momentum due the large initial angular momentum or the initial weak magnetic field, respectively. For slow rotators with strong magnetic fields, discs are more likely to form in the models with $\bm{B}_{\text{-z}}$, indicating that the Hall effect is efficient here. Previous studies that investigates the Hall effect in star formation simulations (Tsukamoto et al., 2015b; Wurster et al., 2016; Tsukamoto et al., 2017; Wurster et al., 2018a, c, d) were initialised in this parameter space.

4.4 Disc formation

There are five models that never form a disc, as per the criteria given in Section 4.1. We examine the simulations every 0.36 kyr, thus it is likely that the disc forms and dissipates within a single timestep in these models. Thus, for these models, we substitute the disc formation time with the sink formation time. With the exception of four models (I$\Omega_{25}\mu_{5}B_{\text{-x}}$, I$\Omega_{25}\mu_{5}B_{\text{-z}}$, I$\Omega_{45}\mu_{3}B_{\text{-z}}$ and I$\Omega_{65}\mu_{3}B_{\text{-z}}$), the discs in all the transient models dissipate within 0.72 kyr of disc formation.

Fig. 3 shows the formation time of the discs. The shapes represent the initial mass-to-flux ratios and the colours represent the initial magnetic field orientation; the initial rotational velocities are off-set from their actual value for clarity.

Discs form later for increasing $\Omega_{0}$ and decreasing $\mu_{0}$ (i.e. increasing $B_{0}$). At a given $\Omega_{0}$ and $\mu_{0}$, discs tend to form at similar times for $\bm{B}_{\pm\text{z}}$ and slightly later for $\bm{B}_{\text{-x}}$.

The disc formation time is approximately independent of the non-ideal MHD processes; during the initial collapse, the gas is not dense enough nor is the magnetic field strong enough for the non-ideal MHD processes to significantly affect the evolution of the system. The discs in the hydrodynamical models tend to form at similar times to the weak field magnetic models.

These trends are expected since magnetic fields and angular momentum both support against collapse, thus increasing these values naturally leads to a later disc formation time. When the magnetic field is initially perpendicular to the axis of rotation (i.e. $\bm{B}_{\text{-x}}$), it resists the vertical collapse. As the vertical collapse proceeds, it drags the magnetic field to the mid-plane, where the magnetic field is further amplified. Both resisting the vertical collapse and the stronger mid-plane magnetic field strength given similar initial conditions results in the models with $\bm{B}_{\text{-x}}$ forming discs after those with $\bm{B}_{\pm\text{z}}$.

In many previous studies (e.g. Price & Bate, 2007; Hennebelle & Fromang, 2008; Duffin & Pudritz, 2009; Hennebelle & Ciardi, 2009; Commerçon et al., 2010; Seifried et al., 2011; Wurster et al., 2016), discs failed to form in models that included ideal MHD and strong magnetic fields, consistent with the magnetic breaking catastrophe (Allen et al., 2003). However, most of our models form discs, including 24 of 40 of our ideal MHD models.

To test for fragmentation, this study intentionally investigates a broad range of initial rotation rates, including the high-end tail as empirically determined by Goodman et al. (1993). Previous studies performed fewer simulations, and were initialised with values more indicative of the mean observed rotation rates $\beta_{\text{r}}\lesssim 0.02$ and observed magnetic field strengths, $\mu_{0}\sim 2-5$. In this small parameter space, we find only one non-transient disc (I$\Omega_{25}\mu_{5}B_{\text{-x}}$) when using ideal MHD. Thus, our results are consistent with the literature in the typically explored parameter space.

With ideal MHD, all of the models initially with more magnetic energy than rotational energy (i.e. $E_{\text{0,mag}}>E_{\text{0,rot}}$) either form transient or small smooth discs, whereas models with $E_{\text{0,mag}}<E_{\text{0,rot}}$ form discs that form spiral arms or fragment. This relation does not hold for non-ideal MHD models.

In summary, increasing the initial rotation has the strongest effect on delaying disc formation, followed by increasing the magnetic field strength, and finally by changing the initial orientation of the magnetic field.

4.5 Disc fragmentation

The top panels in Fig. 4 show the formation time of the first fragment relative to the formation time of the disc, and the bottom panels show the radius at which the fragmentation occurred.

We find that 56 per cent of the models that fragment do so within 4 kyr after disc formation, including all but two of the non-ideal MHD models with $\Omega_{0}\geq 0.65\Omega_{\text{orbit}}$. However, there is no trend of fragmentation time with respect to any of our parameters. Of the models with $\Omega_{0}=0.45\Omega_{\text{orbit}}$ that fragment, they do so $>5$ kyr after disc formation. These discs are initially compact, and remain small in the presence of strong magnetic fields (either low-$\mu_{0}$ or $\bm{B}_{\text{-x}}$). Thus, compared to models with weaker magnetic fields which quickly form large discs, these models require additional time to grow in radius and for the magnetic Toomre-Q parameter, $Q_{\text{m}}$, to decrease into the instability regime where it may then fragment (see Section 5.3 below).

As with the fragmentation time, there is no strong correlation between fragmentation distance and our parameters, although models with lower $\Omega_{0}$ tend to fragment at further distances from the centre. Approximately 78 per cent of the models that fragment do so at initial distances of $r\lesssim 80$ au from the centre of the disc. Prior to fragmentation, the disc becomes unstable and forms spiral arms, and it is in the arms where the fragmentation occurs; in most cases, the fragment forms in the middle of the arm rather than near the bulk disc or near the tip.

The fragmentation distance must be taken with caution since the fragments typically migrate during and after formation, which leads to a large uncertainty about their specific formation distance (e.g. Kley & Nelson, 2012; Baruteau et al., 2014; Meru et al., 2019). Like the fragmentation times, the fragmentation distance is independent of the initial parameters in our suite.

4.6 Disc Properties

The discs are continually evolving, as shown by the time sequence in Fig. 1. The evolution of the bulk disc radii and bulk disc mass are shown by the lines in Figs. 5 and 6, respectively, with the final bulk and total values shown by the points; the final radii and masses are also given in Tables 1 and 2.

In general, the bulk disc radii continue to grow with time as gas accretes on to them, with final radii of 30-50 au at $t-t_{\text{disc}}\approx 16$ kyr. The models that fragment necessarily have smaller final radii since they have had less time to grow. Models with an initially weak magnetic field typically undergo an initial rapid accretion phase followed by slower growth, and the growth rate for the $\mu_{0}={20}$ models is only trivially dependent on the orientation of the magnetic field or the initial rotation. The dependence on initial conditions becomes more pronounced for stronger magnetic fields, however there are fewer non-transient models to analyse. Where there are slight differences amongst the growth rates, the non-ideal MHD models form larger bulk discs than their ideal MHD counterparts.

Fig. 5 also plots the final total radius (defined at the outermost radius where a gas particles has $\rho>10^{-13}$ g cm*-3*, thus typically represents the radius of the tip of the spiral arm). The total radii spans a large range, and is typically larger for models with weak magnetic fields or high rotations. The smooth disc models typically have total radii of $r_{\text{T,disc}}/r_{\text{B,disc}}\lesssim 2$, while the models that form spiral arms and/or fragment can have $2\lesssim r_{\text{T,disc}}/r_{\text{B,disc}}\lesssim 6$. The large ratio and large total radius indicate the presence of substantial spiral arms (independent of whether or not they fragment).

Fig. 5 also gives an indication about the location of the fragment. Fragmentation typically occurs at $1.5\lesssim r_{\text{frag}}/r_{\text{B,disc}}\lesssim 4$, and the majority (25 of 41) of the models fragment closer to the tip of the spiral arms rather than the edge of the bulk disc; however, three models (N$\Omega_{45}\mu_{20}B_{\text{-z}}$, I$\Omega_{65}\mu_{20}B_{\text{-x}}$ and I$\Omega_{85}\mu_{5}B_{\text{-z}}$) fragment on the edge of the bulk disc. Thus, for fragmentation to occur, the prior formation of spiral arms is clearly beneficial.

The disc masses increase over their lifetime, indicating that they are continually accreting gas (Fig. 6). These values are dependent on the accretion onto the disc, accretion onto the star from the disc, and the instantaneous calculation of the radius. Thus, these values have noticeable variability, and it is challenging to extract meaningful trends, although there is a slight trend of lower discs masses for models with faster initial rotations. In all cases, the majority of the gas is in the bulk disc (comparing the symbols at the end of each curve), indicating that the there is very little dense gas in the spiral arms. The models that fragment tend to have higher fractions of mass in the spiral arms (30-40 per cent), suggesting that the fragments that form are (at least initially) low mass.

Throughout the lifetime of the disc-like structures, they are rotating at sub-Kelperian speeds. Towards the end of the simulations, the rotational speed has increased to being only a few times slower than the Keplerian speed.

5 Discussion

This section further discusses fragmentation, however, we also discuss some other characteristics and trends that arise as a direct result of varying our parameter space.

5.1 Hydrodynamic versus magnetised models

At any given $\Omega_{0}$, the purely hydrodynamic models should collapse faster than their strongly magnetised counterparts since there is no support from magnetic fields. Indeed, the discs form at similar times to the weakly magnetised models, where the magnetic field does not play an important role in the evolution of the system. Without magnetic fields to transport angular momentum away from the collapsing central region, rotationally supported discs form, with the size increasing from models with $\Omega_{0}=0.25$ to $0.45\Omega_{\text{orbit}}$; these discs are larger than their magnetised counterparts at any given $\Omega_{0}$. Models with $\Omega_{0}=0.65$ and $0.85\Omega_{\text{orbit}}$ quickly become unstable and fragment, thus a direct comparison of disc size is not feasible.

Since the hydrodynamic models form discs that are larger than their strongly magnetised counterparts (assuming the discs grow rather than immediately fragment), it is reasonable to expect that if a hydrodynamic model does not fragment, then neither will its magnetised counterparts. This is true, with the exception of N$\Omega_{25}\mu_{10}B_{\text{+z}}$. The disc in H$\Omega_{25}$ is larger and more extended, with wide spiral arms close to the bulk disc. The bulk disc in N$\Omega_{25}\mu_{10}B_{\text{+z}}$ is condensed, with a narrow, extended arm in near isolation. These narrow spiral arms are typically more susceptible to instability and fragmentation than wider arms, frequently forming over-densities. In some cases, the over-densities dissipate back into the arm, but the more frequent result is fragmentation, as in the case of N$\Omega_{25}\mu_{10}B_{\text{+z}}$.

5.2 Ideal versus non-ideal MHD

5.2.1 Parallel versus perpendicular magnetic field

Generally, the discs in the $\bm{B}_{\text{-x}}$ models tend to be slightly larger than their $\bm{B}_{\pm\text{z}}$ counterparts. When they do fragment, it tends to be later, suggesting that this orientation of magnetic field stabilises against fragmentation, but typically delays rather than prevents it. The mid-plane magnetic field strengths tend to be stronger in the $\bm{B}_{\text{-x}}$ models than their $\bm{B}_{\pm\text{z}}$ counterparts (in agreement with Wurster et al., 2017b).

5.2.2 Parallel magnetic field

As discussed above, discs are less likely to form in strong magnetic fields in the ideal MHD approximation (e.g. Allen et al., 2003), and indeed, only our fastest rotating ideal MHD model with $\mu_{0}={3}$ forms a disc. By including non-ideal MHD with the $\bm{B}_{\text{-z}}$ orientation, four of the $\mu_{0}={3}$ models form discs, and N$\Omega_{85}\mu_{3}B_{\text{-z}}$ even fragments. By reversing the direction of the magnetic field, the Hall effect transports the angular momentum in the gas around the protostar outwards, hindering disc formation. Thus, N$\Omega_{5}\mu_{3}B_{\text{+z}}$ and N$\Omega_{25}\mu_{3}B_{\text{+z}}$ fail to form discs. However, the magnetic dissipation from Ohmic resistivity and ambipolar diffusion permit discs to form in N$\Omega_{45}\mu_{3}B_{\text{+z}}$ and N$\Omega_{65}\mu_{3}B_{\text{+z}}$ unlike in their ideal MHD counterparts.

At $\mu_{0}={5}$ and $\Omega_{0}\geq 0.45\Omega_{\text{orbit}}$, the discs fragment for the non-ideal MHD models but not the ideal MHD models, as a result of physical resistivity allowing larger discs to form. The Hall effect contributes oppositely to the angular momentum in the discs for the models with $\bm{B}_{\text{+z}}$ and $\bm{B}_{\text{-z}}$, however, the azimuthal ion velocity is similar to the bulk azimuthal velocity, showing that Hall effect cannot overcome the fast initial rotation to make significant changes to the evolution; the fragmentation is only slightly delayed in N$\Omega_{45}\mu_{5}B_{\text{-z}}$ compared to N$\Omega_{45}\mu_{5}B_{\text{+z}}$. The Hall effect is strong enough at $\Omega_{0}=0.25\Omega_{\text{orbit}}$ such that the disc in N$\Omega_{25}\mu_{5}B_{\text{-z}}$ is larger than in N$\Omega_{25}\mu_{5}B_{\text{+z}}$, but neither fragment.

In the ideal MHD models, although the magnetic field is initially $\bm{B}_{\text{-z}}$, a strong toroidal component develops in the discs. This occurs for all the ideal MHD models that form discs (except for I$\Omega_{5}\mu_{20}B_{\text{-z}}$). The non-ideal MHD models also develop a toroidal component, but it is generally weaker than in their ideal counterparts.

5.2.3 Perpendicular magnetic field

The discs are typically larger in the non-ideal MHD models than their ideal counterparts due to less magnetic braking. Unlike the ideal MHD models, discs form for $\Omega_{0}\geq 0.25\Omega_{\text{orbit}}$ and $\mu_{0}={3}$.

Models N$\Omega_{10}\mu_{45}B_{\text{-x}}$ and N$\Omega_{20}\mu_{45}B_{\text{-x}}$ do not fragment, unlike their ideal MHD counterparts. All four form narrow and dense arms, and the arms in the non-ideal MHD models are permeated with weaker magnetic fields than in their ideal counterparts, thus spread out and do not fragment. The spiral arms in the ideal MHD models become more well defined as they evolve until they fragment $>~{}10$ kyr after the discs has formed.

5.3 Toomre-Q parameter

The Toomre-Q parameter (Safronov, 1960; Toomre, 1964) is given by

[TABLE]

where $c_{\text{s}}$ is the local sound speed, $\kappa$ is the epicyclic frequency and $\Sigma$ is the surface mass density. In the presence of magnetic fields, the magnetic Toomre-Q parameter is given by

[TABLE]

where $v_{\text{A}}$ is the Alfvén velocity.

The epicyclic frequency is given by

[TABLE]

where $\Omega$ is the angular frequency of the disc. For Keplerian discs, $\kappa\approx\Omega$, and this is the version that is commonly presented in the literature. However, during disc formation, the young discs rotate with sub-Keplerian speeds, thus we calculate the Toomre-Q parameter using $\kappa$ rather than $\Omega$. The fluid becomes unstable for $\kappa^{2}<0$.

Discs are susceptible to fragmentation when $Q\lesssim Q_{\text{crit}}$. For an infinitesimally thin hydrodynamics disc, $Q_{\text{crit}}\sim 1$, while for a 3D hydrodynamics disc $Q_{\text{crit}}\sim 1.5-1.7$ (e.g. Durisen et al. 2007, Helled et al. 2014 and references therein). Using 2D shearing box simulations, Kim & Ostriker (2001) determined that for magnetised self-gravitating discs, $Q_{\text{m,crit}}\sim 1.2-1.4$, where this range is given for the growth of non-axisymmetric perturbations (see also Kim et al. 2003).

The left-hand column of Fig. 7 shows the azimuthally averaged $Q_{\text{m}}$ at six times for the three representative cases shown in Fig. 1 that form discs. Given the asymmetric nature of many of our discs, we also calculate $Q_{\text{m}}$ in wedges of $24^{\circ}$ to search for local minima that may not be detectable in the azimuthally averaged values; see the right-hand column of Fig. 7 for the $Q_{\text{m}}$ wedge that contains the minimum value.

In our suite of models, $Q_{\text{m}}$ yields limited insight into disc fragmentation. As our smooth example disc evolves (top panel of Fig. 7), the minimum $Q_{\text{m}}$, given by $Q_{\text{m,min}}$, slowly decreases. When considering the wedge, the values are slightly lower suggesting a slight asymmetry in the disc, but there is still no indication that this disc will fragment (visually confirmed in the second row of Fig. 1).

Both N$\Omega_{45}\mu_{20}B_{\text{-x}}$ (middle panels) and N$\Omega_{85}\mu_{20}B_{\text{-z}}$ (bottom panels) have $Q_{\text{m,min}}<1.4$ in the global (azimuthally averaged) and wedge profiles. The global value in the former only briefly drops below the 2D stability limit, suggesting that the over-density quickly diffuses rather than collapses and fragments. The global value in the latter is only marginally unstable at the radius where the disc ultimately fragments. When considering the wedges, both models have several regions that are Toomre-unstable, which suggest that both models have regions that are susceptible to fragmentation. However, only the latter model fragments. Thus, we cannot clearly determine the outcome of a model based upon $Q_{\text{m}}$ alone.

Contrary to the top panel of Fig. 7, many of the stable disc models yield regions of $Q_{\text{m}}<1.4$ but do not fragment. Moreover, most of the discs in our suite have regions that are Toomre-unstable, even those that form smooth discs or discs with spiral arms. Thus, we conclude that all models that fragment do so in regions with $Q_{\text{m}}<1.4$, but not all regions with $Q_{\text{m}}<1.4$ necessarily fragment.

To compare the effect of the magnetic field orientation, Fig. 8 shows $Q_{\text{m}}$ for the magnetised models with $\mu_{0}={5}$ and $\Omega_{0}=0.45\Omega_{\text{orbit}}$ 4 kyr after disc formation; the gas column density of these models at this time is show in Fig. 9.

At this time, these models have formed rotationally supported discs. The azimuthally averaged $Q_{\text{m}}$ suggests that these discs are stable, however, the wedge values suggest there are unstable regions for each model except I$\Omega_{45}\mu_{5}B_{\text{-x}}$. Each of these locations of $Q_{\text{m,min}}<1.4$ corresponds to an over-density clearly visible in Fig. 9, with the exception of H$\Omega_{45}$. These over-densities do not correspond to an increase/decrease in any other property, including gas temperature, velocity, magnetic field strength or velocity, and quickly disperse. Models N$\Omega_{45}\mu_{5}B_{\text{+z}}$, N$\Omega_{45}\mu_{5}B_{\text{-z}}$ and H$\Omega_{45}$ ultimately fragment at $t\approx 8$, 10 and 16 kyr after disc formation, respectively.

Thus, the orientation of the magnetic field clearly influences the formation and evolution of the disc, even early on. Most of the example models form clear over-densities, but these do not fragment, again suggesting that $Q_{\text{m}}$ alone does not determine if a disc will fragment.

5.4 Disc-to-star mass ratios

Massive discs are more likely to fragment than less massive discs, and fragmentation is expected for

[TABLE]

where $H=c_{\text{s}}/\Omega$ is the isothermal scale height (Gammie, 2001).

In isothermal, Keplerian discs, it is expected that $H/R\sim 0.1$. Our discs are not isothermal, with the sound speed varying by a factor of $\sim$10 throughout the disc. Taking the average sound speed and average angular frequency over the disc, we find that $0.1\lesssim H/R\lesssim 0.15$ for $t>t_{\text{disc}}+4$ kyr for the bulk discs, and $0.05\lesssim H/R\lesssim 0.1$ for the total discs. Thus, on average, $H/R\sim 0.1$ is an appropriate approximation in our models, which is in general agreement with the literature.

For this calculation, the stellar mass, $M_{\text{star}}$, is given by the mass of the sink particle of radius 1 au, and the disc mass is that of the bulk disc. After the initial accretion phase, the mass ratio is typically ${M_{\text{disc}}}/{M_{\text{star}}}\sim 1$ (this is not true for the models with $\Omega_{0}=0.05\Omega_{\text{orbit}}$ and the ideal MHD models with $\mu_{0}={3}$). This is reasonable since the protostellar discs have just formed, and the central regions have just collapsed to form the protostar. For the duration of the simulation, this ratio is approximately constant, but we cannot comment upon the long term evolution of this ratio.

In general, ${H}/{R}\sim 0.1{M_{\text{disc}}}/{M_{\text{star}}}$ for the models that form spiral arms or fragment, thus in these models, the condition in Eqn. 13 is satisfied. However, these models and many of the smooth models maintain ${H}/{R}<{M_{\text{disc}}}/{M_{\text{star}}}$ throughout the simulation, thus this relationship is a poor discriminate to determine fragmentation.

5.5 Outflows

Although not the main focus of this study, the large parameter space allows us to briefly investigate outflows. Since our models use 1 au sink particles, these outflows are first core outflows (for a more detailed discussion on first core outflows in non-ideal MHD simulations, see Wurster et al., 2018a, c).

Slow ($v_{\text{r}}<1$ km s*-1*), broad outflows are launched in the ideal MHD models with $\mu_{0}\leq{5}$, $\Omega_{0}\geq 0.25\Omega_{\text{orbit}}$ and $\bm{B}_{\text{-z}}$. Slow outflows are also launched in N$\Omega_{25}\mu_{3}B_{\text{-z}}$ and N\Omega_{85}\mu_{3}B_{\text{\pm z}}, while fast ($v_{\text{r}}\sim 1-5$ km s*-1*) outflows are launched in N$\Omega_{5}\mu_{3}B_{\text{+z}}$, N$\Omega_{25}\mu_{3}B_{\text{+z}}$, N\Omega_{45}\mu_{3}B_{\text{\pm z}}, N\Omega_{65}\mu_{3}B_{\text{\pm z}}. Four of these ideal MHD and two of the non-ideal MHD models have transient classifications. No outflows are launched in our models with $\bm{B}_{\text{-x}}$.

From the ideal MHD models, this suggests that magnetic fields with a strong poloidal component are required in addition to a reasonable amount of initial angular momentum. Ohmic resistivity and ambipolar diffusion weaken the magnetic field enough such that outflows are not launched in the $\mu_{0}={5}$ models that include the non-ideal MHD processes, and that fast outflows are launched in six of the nine non-ideal MHD models that launch outflows. Thus, we find that outflows are dependent primarily on direction and strength of the magnetic field, where the strength is necessarily weakened by the inclusion of the non-ideal MHD processes. We generally find outflow speeds decreasing with increasing $\Omega_{0}$, since larger discs are permeated by a similar magnetic flux as smaller discs which results in less magnetic pinching and weaker field strengths in the larger discs.

This result suggests a resolution of conflicting results in the literature. Both Wurster et al. (2017b) and Kuruwita et al. (2017) modelled the formation of binary stars, however, only the models in Kuruwita et al. (2017) launched outflows. The system in Kuruwita et al. (2017) yielded smaller binary separations and smaller discs than Wurster et al. (2017b), and outflows that carried angular momentum away from the protostars were launched. The calculations of Wurster et al. (2017b) produced large discs and the wide binary separations. To verify that these differences were a result of the initial conditions and not a difference in the algorithm (i.e. SPH vs adaptive mesh refinement), we previously ran low-resolution proof-of-concept models using the algorithms from Wurster et al. (2017b) and the initial conditions from Kuruwita et al. (2017) and found that, indeed, outflows were formed. Thus, large discs in ideal MHD simulations do not appear to launch early outflows.

5.6 Counter-rotating envelopes

It has been previously shown that models that include the Hall effect and initial magnetic field orientations of $\bm{B}_{\text{-z}}$ produce counter-rotating envelopes (e.g. Krasnopolsky et al., 2011; Li et al., 2011; Tsukamoto et al., 2015b; Wurster et al., 2016; Tsukamoto et al., 2017; Wurster et al., 2018b, c). This is to conserve angular momentum as the Hall effect spins up the gas around the protostar. Since the Hall effect hinders disc formation in models with $\bm{B}_{\text{+z}}$, these studies found no counter-rotation in models with $\bm{B}_{\text{+z}}$.

In our entire suite, no counter-rotating envelopes form for $\Omega_{0}\geq 0.25\Omega_{\text{orbit}}$, since the Hall effect is not strong enough to overcome the initial rotation of the envelope. As expected from the previous studies, models I$\Omega_{5}\mu_{*}B_{\text{-z}}$ and N$\Omega_{5}\mu_{*}B_{\text{+z}}$ do not form counter-rotating envelopes; N$\Omega_{5}\mu_{3}B_{\text{-z}}$ forms a strong counter rotating envelope, and N$\Omega_{5}\mu_{5}B_{\text{-z}}$ forms a weak one. See the first three panels of the top row in Fig. 10, which shows the azimuthal velocity $v_{\phi}$ and gas density of the $\Omega_{0}=0.05\Omega_{\text{orbit}}$ models with $\mu_{0}={5}$ near the end of their respective simulations; note that $v_{\phi,0}>0$.

In the above cases, and typically discussed throughout the literature, the cause of the counter-rotating envelope is the Hall effect when the magnetic field and rotation vectors are anti-aligned. However, counter-rotating regions may also form if the magnetic field is initially perpendicular to the rotation axis (i.e. $\bm{B}_{\text{-x}}$).

A gas over-density forms along the rotation axis in I$\Omega_{5}\mu_{3}B_{\text{-x}}$ (fourth column in Fig. 10). The low-density gas in the mid-plane rotates rapidly, while above and below form a slow counter-rotating envelope. In all other ideal MHD models the initial rotation is strong enough to prevent a counter-rotating envelop from forming.

In N$\Omega_{5}\mu_{3}B_{\text{-x}}$, the Hall effect causes the gas to misalign from the rotation axis such that the normal to the dense disc is misaligned by $\sim$40∘. The gas is still infalling along the initial rotation axis, and is counter-rotating along the plane of the disc (i.e. $x\sim-y$). These results are similar to that found in Tsukamoto et al. (2017).

Thus, there is a very small parameter space in which counter-rotating envelopes may form. This required parameter space must include an initially slowly rotating envelope and strong magnetic fields that are either $\bm{B}_{0}=\bm{B}_{\text{-x}}$, or $\bm{B}_{0}=\bm{B}_{\text{-z}}$ if the Hall effect is included.

5.7 Resolution

Our models have been performed using a constant resolution of $10^{6}$ particles in the sphere, thus, we cannot explicitly discuss convergence, however, we will briefly comment on resolution.

Our mass resolution is $m_{\text{p}}=10^{-6}$ M⊙ per particle, thus there are \sim$$10^{5} particles in our discs by the end of the simulation (recall Fig. 6). Given this number of particles, the Jeans mass is still resolved (Bate & Burkert, 1997, recall Section 3).

For discs, the Toomre-mass,

[TABLE]

must be resolved to prevent numerically induced fragmentation. From this equation, Nelson (2006) calculated the maximum resolvable surface density to be

[TABLE]

which we have modified to

[TABLE]

for our magnetised models. Here, $N_{\text{reso}}$ is the number of particles required to resolve this maximum surface density, which Nelson (2006) empirically determined to be $N_{\text{reso}}\sim 6N_{\text{neigh}}$; given our cubic spline kernel, $N_{\text{reso}}\sim 342$ particles. Throughout our suite, we find that the surface density of the discs is a few orders of magnitude lower than $\Sigma_{\text{m,max}}$, hence our discs meet the Toomre-mass criterion.

Resolving the vertical structure of discs is also important. A poorly resolved vertical structure will underestimate the mid-plane density and hence gas pressure, which will inflate the discs and increase the likelihood of fragmentation. For SPH simulations, Nelson (2006) determined that at least four smoothing lengths $h$ are required per scale-height at the disc mid-plane. Analogous to Section 5.4, we calculate the scale height using

[TABLE]

Using the mid-plane smoothing length, $H/h<4$ for $r\lesssim 10$ au, thus, the inner regions of the discs do not meet this criteria and may be under resolved. This is to be expected given the presence of, and boundary effects caused by, the sink particle (e.g. Machida et al., 2014; Wurster et al., 2017a). An under-resolved inner disc should not affect our general conclusions since the disc is not expected to fragment at such small radii. For $r\gtrsim 10$ au, $H/h>4$ is typically satisfied, thus, we can be confident that our discs are vertically resolved.

5.7.1 Convergence studies in the literature

Although our discs are resolved (via the Jeans mass, Toomre-mass, and scale-height), resolution may still affect our results. Convergence studies, especially of disc formation and fragmentation, have been performed frequently throughout the literature, and these studies have suggested that the fragmentation results are resolution-dependent (e.g. Meru & Bate, 2011a, b, 2012; Forgan et al., 2017; Meyer et al., 2018). While most studies have used parameterised cooling rates, Meyer et al. (2018) performed radiation hydrodynamics simulations and also found that decreasing the resolution prevents fragmentation. While increasing the resolution allows the disc to initially fragment at similar times in their fiducial and high resolution simulations, the future evolutions diverge, with more fragments forming at high resolution than fiducial resolution.

Independent of resolution, all of the hydrodynamic models in Forgan et al. (2017) fragmented. For their magnetised models, increasing their particle number by a factor of two yielded an additional model that fragmented. Resolution also played a role in the disc formation study of Wurster et al. (2016): When modelling ideal MHD, smaller discs formed in higher resolution simulations, while in the models that included the non-ideal MHD processes, the disc masses differed by less than 5 per cent between their two resolutions.

The above studies show that the convergence issue is persistent. This issue arises in models starting from a pre-initialised disc and a molecular cloud core, studies investigating low-mass and high-mass star formation, and is independent of numerical method (i.e. SPH or a grid code).

Since our study meets the resolution criteria discussed above, models radiation hydrodynamics, and we do not model the evolution of the fragments or subsequent fragmentation, we believe that our qualitative results will be unaffected by resolutions. Furthermore, non-ideal MHD models tend to be less sensitive to resolution than their ideal MHD counterparts. However, for definitive quantitative results, a convergence study would be needed.

5.8 Comparison to other studies

Most fragmentation studies, including ours, show that fragmentation does not occur at small disc radii (e.g. Stamatellos & Whitworth, 2008; Boley, 2009; Clarke, 2009; Forgan et al., 2017). At low opacities, however, fragmentation is possible at small radii due to the higher cooling rate of dense gas and the shorter lifetime of the first hydrostatic code (e.g. Meru & Bate, 2010; Bate, 2014, 2019).

Magnetic fields are typically neglected in disc fragmentation studies, however, Forgan et al. (2017) initialise their discs with a toroidal magnetic field. At the longest cooling time, no magnetised disc fragmented (although the hydrodynamics discs did), while at the shorter cooling times, all the magnetised discs fragmented. Their discs were initialised to produce a considerable number of fragments (if unstable), rather than the few fragments that formed in the models we present here. Nonetheless, they find that a large fraction of the hydrodynamical disc fragmented, whereas the magnetised discs fragmented in a narrow band. In our models, there is no significant difference in the fragmentation distance between the hydrodynamic and magnetised models (see Fig. 4). Similar to Forgan et al. (2017), we find the inclusion of magnetic fields stabilises the discs, and either delays or prevents fragmentation.

One significant difference between the two studies is that their hydrodynamical and magnetised discs were initialised with the same parameters (e.g. radius, surface density profile) whereas our discs formed self-consistently, thus our magnetised discs were typically much smaller and less massive than the hydrodynamic discs, and this size difference in our models also likely contributed to the delayed or prevention of disc formation.

Although not explicitly investigating the fragmentation limit, Zhao et al. (2018) investigated disc formation and fragmentation in the presence of magnetic fields, and included Ohmic resistivity and ambipolar diffusion. Similar to here, they concluded that faster rotating discs promoted fragmentation, and that the disc was more likely to fragment in the presence of weak magnetic fields compared to strong fields. They also find a diversity in where the fragments form and how they evolve.

6 Summary and conclusion

We have presented a suite of simulations studying the formation and fragmentation of discs around protostars in the presence of magnetic fields. Our models were initialised as 1M⊙ Bonnor-Ebert spheres, which collapsed to form protostars typically surrounded by massive protostellar discs. We followed the evolution until the final classification – fragmented, spiral arms, smooth, or transient – could be determined. Our suite included ideal MHD, non-ideal MHD and purely hydrodynamical models, where the non-ideal MHD models included Ohmic resistivity, ambipolar diffusion and the Hall effect. We tested five initial rotations $\Omega_{0}$, four initial magnetic field strengths $\mu_{0}$, and three (non-ideal MHD) or two (ideal MHD) orientations of the magnetic field. Our simulations were radiation hydrodynamics simulations that were performed using the SPH code sphng.

Our key results are as follows:

Disc-like structures (herein referred to ‘discs’) formed later for models with faster initial rotations and/or stronger magnetic fields. 2. 2.

Of our 105 models, 41 fragmented, 21 formed spiral structures but did not fragment, 12 formed smooth discs, and 31 did not form discs. Discs were more likely to fragment for faster initial rotations (i.e. $\Omega_{0}\geq 0.45\Omega_{\text{orbit}}$), and for weaker magnetic fields (i.e. $\mu_{0}\geq{10}$). Non-ideal effects became important for strong magnetic fields (i.e. $\mu_{0}\leq{5}$), and many of these discs with $\bm{B}_{0}=\bm{B}_{\pm\text{z}}$ fragmented while their ideal MHD counterparts did not. 3. 3.

For the discs that fragmented, there was no clear correlation between fragmentation time or distance and our initial parameters (rotation rate, magnetic field strength and orientation). 4. 4.

The magnetic Toomre-Q parameter, $Q_{\text{m}}$, cannot be used in isolation in our models to determine if a model will fragment. All models that fragment do so in a region of $Q_{\text{m}}<1.4$, but not all models that entered this regime fragmented. When comparing the ratio of disc-to-stellar masses, the ratio suggested that all of our discs were unstable to fragmentation. This is a result of the young protostar that is still accreting mass from the disc. 5. 5.

Outflows were launched from models with initially strong, $\bm{B}_{0}=\bm{B}_{\pm\text{z}}$ magnetic fields. Fewer models that include the non-ideal MHD processes launch outflows since these processes weaken the magnetic field. 6. 6.

Counter-rotating envelopes form only under specific conditions: an initially slowly rotating envelope with strong magnetic fields that are either $\bm{B}_{0}=\bm{B}_{\text{-x}}$, or $\bm{B}_{0}=\bm{B}_{\text{-z}}$ if the Hall effect is included. 7. 7.

Discs masses up to $\sim$0.1 M⊙were obtained (i.e., up to $\sim$10 per cent of the SPH particles that were initially in the sphere). The Jeans mass and Toomre mass were resolved throughout the calculations, and the scale-height was resolved for $r>10$ au. However, a proper resolution study would be required for a thorough discussion of convergence.

Given our results, we cannot predict when or where a disc will fragment. However, qualitatively, fragmentation is promoted in faster rotating models and in models with weaker magnetic fields.

Acknowledgements

We would like to thank the referee for useful and insightful comments that 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). The calculations for this paper were performed on the University of Exeter Supercomputers, Zen and Isca. The former was a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter; the latter is part of the University of Exeter High-Performance Computing (HPC) facility. Analyses were performed the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. Additional analyses were performed using using 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 column density figures were made using splash (Price, 2007). The research data supporting this publication are openly available from the University of Exeter’s institutional repository, https://ore.exeter.ac.uk/repository.

Appendix A Appendix

Tables 1 and 2 summarise the models at the end of the simulation. The simulations are ended once the final classification can be determined. A few models have small, transient asymmetries; although these features are not well-defined and persistent spiral arms, it means the disc is not completely smooth. For simplicity, we give these models the classification of the most persistent state, and identify them with an asterisk after their classification in the following tables. Along with its classification, each model is listed with its disc formation time, the end-time relative to the disc formation time, and the radius and mass of both the bulk and total discs at this time. The disc properties are calculated as described in Section 4.1; recall that ‘disc’ refers to ‘disc-like structure.’ For the models that fragment, the data we present is from the first output that contains the fragment; for the transient models, the data we present is from the first output in which a disc is not present. If a disc never forms, we present the sink formation time, and the relative end-time is set to zero. Also included is the stellar mass (i.e the mass of the sink particle). The penultimate column indicates if a fast ($v_{\text{r}}\sim 1-5$ km s*-1*) or slow ($v_{\text{r}}<1$ km s*-1*) outflow existed at any time during the simulation, and the final column indicates if a counter-rotating envelope (CRE) existed at any time during the simulation.