The dependence of protostar formation on the geometry and strength of the initial magnetic field

TL;DR

This study uses simulations to explore how the initial magnetic field's strength and geometry influence protostar formation, outflow structures, and system configurations, revealing significant effects on star and binary formation processes.

Contribution

It provides new insights into how magnetic field parameters affect protostellar collapse, outflows, and system morphology, using detailed magnetohydrodynamics simulations.

Findings

Magnetic field strength influences binary formation likelihood.

Field geometry affects outflow suppression and disc warping.

Outflow speeds do not exceed 8 km/s in the simulations.

Abstract

We report results from twelve simulations of the collapse of a molecular cloud core to form one or more protostars, comprising three field strengths (mass-to-flux ratios, {\mu}, of 5, 10, and 20) and four field geometries (with values of the angle between the field and rotation axes, {\theta}, of 0{\deg}, 20{\deg}, 45{\deg}, and 90{\deg}), using a smoothed particle magnetohydrodynamics method. We find that the values of both parameters have a strong effect on the resultant protostellar system and outflows. This ranges from the formation of binary systems when {\mu} = 20 to strikingly differing outflow structures for differing values of {\theta}, in particular highly suppressed outflows when {\theta} = 90{\deg}. Misaligned magnetic fields can also produce warped pseudo-discs where the outer regions align perpendicular to the magnetic field but the innermost region re-orientates to be perpendicular to the rotation axis. We follow the collapse to sizes comparable to those of first cores and find that none of the outflow speeds exceed 8 km s$^{-1}$. These results may place constraints on both observed protostellar outflows, and also on which molecular cloud cores may eventually form either single stars and binaries: a sufficiently weak magnetic field may allow for disc fragmentation, whilst conversely the greater angular momentum transport of a strong field may inhibit disc fragmentation.

Equations

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla {}^i{v}^i\text{leavevmode}\text{,}$$

$$\frac{\mathrm{d}{v}^i}{\mathrm{d}t}=\frac{1}{\rho }{\nabla }^j{S}^{ij}-{\nabla }^i\varphi \text{leavevmode}\text{,}$$

$$\frac{\mathrm{d}{B}^i}{\mathrm{d}t}=\left(B^j{\nabla }^j\right)v^i-B^i\left({\nabla }^j{v}^j\right)\text{leavevmode}\text{,}$$

$${\nabla }^2\varphi =4\pi G\rho \text{leavevmode}\text{,}$$

$$S^{ij}=-P{\delta }^{ij}+\frac{1}{{\mu }_0}\left(B^i{B}^j-\frac{1}{2}{\delta }^{ij}{B}^2\right)\text{leavevmodenobreak .}$$

$$\frac{\mathrm{d}\text{leavevmode}}{\mathrm{d}t}=\frac{\partial \text{leavevmode}}{\partial t}+v^i{\nabla }^i\text{leavevmode}\text{,}$$

$$h_a=\eta {\left(\frac{m_a}{{\rho }_a}\right)}^{\frac{1}{\nu }}\text{leavevmodenobreak ,}$$

$$\tau =\frac{h}{\sigma c_{\mathrm{c}}}\text{leavevmodenobreak ,}$$

$$\begin{array}{rl}{c}_{\mathrm{c}}^2& =c_{\mathrm{s}}^2+v_{\mathrm{a}}^2=c_{\mathrm{s}}^2+\frac{B^2}{{\mu }_0\rho }\text{leavevmodenobreak .}\end{array}$$

$$\rho$$

$${\rho }_{\mathrm{c},1}{\left(\frac{\rho }{{\rho }_{\mathrm{c},1}}\right)}^{\frac{7}{5}}$$

$${\rho }_{\mathrm{c},1}{\left(\frac{{\rho }_{\mathrm{c},2}}{{\rho }_{\mathrm{c},1}}\right)}^{\frac{7}{5}}{\left(\frac{\rho }{{\rho }_{\mathrm{c},2}}\right)}^{\frac{11}{10}}$$

$$\mu =\frac{{\varpi }_{\mathrm{sphere}}}{{\varpi }_{\mathrm{critical}}}$$

$${\varpi }_{\mathrm{sphere}}=\frac{M_{\mathrm{sphere}}}{\pi r_0^2{B}_0}$$

$${\varpi }_{\mathrm{critical}}=\frac{2c_1}{3}\sqrt{\frac{5}{\pi G{\mu }_0}}.$$

$$B_0^x=B_0\sin \vartheta \text{leavevmodenobreak ,}$$

$$B_0^z=B_0\cos \vartheta \text{leavevmodenobreak ,}$$

$$B_r=\sqrt{B_x^2+B_y^2}$$

$$B_{\mathrm{p}}=\sqrt{B_r^2+B_z^2}\text{,}$$

$$B_{\varphi }=B_{\mathrm{tor}}=\arctan \frac{B_y}{B_z}\text{.}$$

Full text

††pagerange: The dependence of protostar formation on the geometry and strength of the initial magnetic field–References††pubyear: 2015

1 Introduction

Molecular clouds, and hence the molecular cloud cores which collapse to form protostars, are magnetised (Crutcher et al., 1993; Crutcher, 2012), as are protostars and the structures surrounding them. In particular, magnetic outflows are proposed as a mechanism to remove angular momentum from protostars — some mechanism to do so is necessary since simulations without magnetic fields (Norman & Wilson, 1978; Bonnell & Bate, 1994) and (semi-)analytic calculations (e.g. Terebey et al., 1984) show that the total angular momentum of cores are too high to form stars without some angular momentum transport mechanism, and Hartmann et al. (1986) find that the rotation rates of T Tauri stars are substantially lower than the energy required to break up the stellar core. If, for certain field strength or geometries, these outflows can be altered or suppressed, this will constrain how star formation proceeds.

Conversely, some cores must fragment to produce binaries, while still removing sufficient angular momentum to form stellar cores which do not break up. By varying the field strength a regime whereby some magnetic braking, plus a magnetic outflow, removes enough angular momentum to allow stellar core formation while still forming a binary may be obtained.

Misalignment between magnetic fields and outflows have been observed on scales of molecular clouds (e.g. Stephens et al., 2014; Hull et al., 2013a, b) down to stellar length-scales, (e.g. Donati et al., 2010). A range of mass-to-flux ratios are observed in molecular clouds, ranging from the observation by Crutcher (2005) that the mean mass-to-flux ratio, expressed in terms of the critical value, for massive star–forming regions is approximately unity, to Beuther et al. (2010) and Giannetti et al. (2013). The latter found a much larger range of mass-to-flux ratios in star forming regions, including both sub- and super-critical clouds, a result reinforced by the Bayesian analysis of Crutcher et al. (2010).

From its invention (Lucy, 1977; Gingold & Monaghan, 1977) the smoothed particle hydrodynamics (SPH) method has been applied to star formation. In Lewis et al. (2015) we used a smoothed particle magnetohydrodynamics (SPMHD) method to model the collapse of a molecular cloud cores to the scale of first hydrostatic cores (Larson, 1969), and the nature of their outflows, with the magnetic field both aligned and misaligned to the rotation axis. We found a significant dependence on the angle between the initial magnetic field and rotation axes (hereafter, $\vartheta$) and the nature of the outflow produced (or, indeed, whether an outflow is produced at all). This result was similar to that obtained by Ciardi & Hennebelle (2010) (who used a Godunov adaptive mesh refinement code rather than SPH), albeit with some differences. In particular, they found that the mass ejection rate reduced as $\vartheta$ was increased, subsequently increasing the protostellar accretion rate, while we observed that for $20\degree\leq\vartheta\leq 60\degree$ the accretion rate decreased.

Lewis et al. (2015) adopted an initial magnetic field strength corresponding to a mass-flux ratio of $\mu=5$. While mass–to–flux ratios close to the critical value are observed, it is instructive to test the effects of weaker fields. Previous work has indicated that the nature of the outflow in an aligned model ($\vartheta=0\degree$) changes with magnetic field strength. This is to be expected — purely hydrodynamic simulations do not produce outflows, so the nature and strength of these outflows must be related to the strength of the magnetic field (see, for example, Price et al., 2003). Bate et al. (2014) found that whilst the velocity of the first-core outflow is essentially unchanged, the width of the outflow (i.e. the opening angle) increases when the initial field is weaker. In addition, stronger and weaker fields will affect the initial collapse of the molecular cloud core (and may prevent it collapsing at all if a sufficient magnetic pressure can be realised).

Therefore, in this work we perform a similar series of calculations to Lewis et al. (2015), but with weaker fields — corresponding to mass-flux ratios of $\mu=10$ and $20$. We present SPMHD simulations of three mass–to–flux ratios and four field geometries. This allows us to probe a range of field strengths. Additionally we can observe how the changing field geometries changes the nature and extent of the outflows observed. We do not include non–ideal MHD effects (see Wardle, 2007; Wurster et al., 2014, 2016; Tsukamoto et al., 2015a, b) and therefore can not probe the regime where $\mu$ is less than or equal to unity.

Section 2 details our numerical method, which is essentially that used in Bate et al. (2014) and section 3 details the initial conditions. Then in sections 4.1, 4.2, 4.3 and 4.4 we present the results of all twelve models and a discussion of the conclusions obtained from them. We then consider the process of accreting material onto the sink particles in section 4.5, before finally comparing our results to some observations in section 4.6.

2 Method

We solve the equations of ideal magnetohydrodynamics with gravity, viz.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with the MHD stress tensor

[TABLE]

In these equations and the rest of the paper, we represent the total derivative by

[TABLE]

the symbols $P$, $\rho$, $v^{i}$, $B^{i}$, and $\phi$ represent the pressure, density, velocity, magnetic and gravitational fields respectively, $\delta^{ij}$ is the Kronecker delta, $\mu_{0}$ is the permeability of free space, and repeated indices imply summation.

We use sphng, a three–dimensional hybrid openmp and mpi smoothed particle radiation magnetohydrodynamics code that originated from a code written by Benz (1990), with subsequent extensive modifications by Bate et al. (1995) and others. We discretise these equations using the SPMHD method derived originally in Price & Monaghan (2004) (see also the review by Price, 2012) and subsequently modified as detailed below. We use the source term correction to the tensile instability proposed by Børve et al. (2001), with the parameter $\chi=1$ (see Lewis et al., 2016). Each particle has an individual smoothing length, $h_{a}$, given by

[TABLE]

where the symbol $a$ denotes individual SPH particles, $m_{a}$ represents the mass of an SPH particle – which is identical for all SPH particles in our model, $\nu=3$ is the number of spatial dimension and $\eta=1.2$ for the cubic spline kernel (Monaghan & Lattanzio, 1985).

To capture shocks and prevent inter–particle penetration we use artificial viscosity and resistivity terms following Monaghan (1997), Price & Monaghan (2004) and Price & Monaghan (2005). Both artificial dissipation terms are controlled by switches which are spatially and temporally variable, for artificial viscosity we use the Morris & Monaghan (1997) switch with $\alpha_{\text{AV}}\in[0.1,1.0]$ while for artificial resistivity the Tricco & Price (2013) switch with $\alpha_{\text{B}}\in[0.0,1.0]$.

The gravitational force is calculated using the method in Price & Monaghan (2007) with the potential softened using the SPH smoothing kernel, where we use the same binary tree to both find neighbours and to compute the gravitational force.

To maintain the solenoidal constraint on the magnetic field (since $\nabla^{i}B^{i}=0$ only enters the equations of MHD as an initial condition and $\nicefrac{{\mathrm{d}\leavevmode\nobreak\ }}{{\mathrm{d}t}}\nabla^{i}B^{i}$ is not guaranteed to be zero) we use the hyperbolic divergence cleaning derived by Tricco & Price (2012) (based on the earlier Eulerian method of Dedner et al., 2002). We set the damping time-scale,

[TABLE]

such that the cleaning wave is critically damped by choosing $\sigma=0.8$, and set the cleaning wave speed to the maximum magnetosonic wave speed, i.e.

[TABLE]

These calculations were performed using the original cleaning method, as opposed to the modified scheme using $\nicefrac{{\psi}}{{c_{\mathrm{c}}}}$ recently proposed by Tricco et al. (2016).

Each SPH particle is assigned an individual timestep based on the local Courant condition. The equations are discretised in time using a second–order Runge–Kutta–Fehlberg integrator (Fehlberg, 1969). Depending on the exact parameters employed, each calculation ran for between ca. 100 and 400 hours of wall time on two 6-core hyper-threaded CPUs (i.e. 24 execution threads, and therefore between 2,400 and 9,600 hours of computation-core time).

The integrator bug that affected the calculations of Lewis et al. (2015), reported in the Erratum and Addendum to that paper (Lewis et al., 2017), was fixed before the calculations for this paper were performed. Consequently, we also do not use the ‘average $h$’ terms that were used in that previous work.

We add sink particles to the simulation once a critical density of $\rho_{\mathrm{crit}}=10^{-10}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ is exceeded and the usual tests (see Bate et al. (1995) for details) are passed. By doing so we are able to avoid resolving the formation of the stellar core itself and therefore can evolve the simulations further. The accretion radius of the sink particle is set to 1 au, any SPH particle that comes within this distance of the sink and passes the Bate et al. (1995) tests will be removed from the simulation and its mass will be added to the sink particle. Other than through gravity, sink particles do not interact in the simulation, in particular they do not possess a hydrodynamic pressure nor a magnetic field.

3 Initial Conditions

The initial conditions are essentially the same as in previous work, which ranges from the early work of Hosking (2002), and the Euler potential SPMHD method in Price & Bate (2007), to the self-consistent direct induction SPMHD method in Price et al. (2012); Bate et al. (2014); Lewis et al. (2015). Calculations with different set ups have also been performed, for example, Banerjee & Pudritz (2006) who used Bonner-Ebert (Bonnor, 1956; Ebert, 1955) spheres instead of the uniform density sphere used here. An $M_{\mathrm{sphere}}=1\leavevmode\nobreak\ \text{M}_{\odot}$ sphere of cold gas was placed in a warm external medium. The sphere and medium are in pressure equilibrium, however, the sphere is 30 times the density of the medium, therefore the initial temperatures are $\approx 10\leavevmode\nobreak\ \text{K}$ (corresponding to an initial isothermal sound speed of $c_{\mathrm{s}}=2.2\times 10^{4}\leavevmode\nobreak\ \text{cm}\thinspace\text{s}^{-1}$) and $\approx 300\leavevmode\nobreak\ \text{K}$ respectively. The sphere had an initial radius $r_{0}=4\times 10^{16}\leavevmode\nobreak\ \text{cm}$ and therefore an initial density of $\rho_{0}=7.4\times 10^{-18}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$. The sphere was then placed in a cubic periodic box with side lengths $16\times{}10^{16}\leavevmode\nobreak\ \text{cm}$ to prevent numerical artefacts in the magnetic field at the edge of the warm medium.

An effective minimum resolution is provided by the Bate & Burkert (1997) criterion, namely that the Jeans length must be resolved throughout the collapse of the core, which would require ca. 30,000 particles; however, Federrath et al. (2011) find that magnetised turbulence may require a resolution well in excess of this. Bate et al. (2014) find no significant resolution dependence between 1 and 10 million particles; consequently we use 1.5 million SPH particles in the core as a compromise between computational expense and numerical accuracy. The sphere is then surrounded by ca. 700,000 particles to represent the warm medium. Aside from providing a ‘container’ for the sphere and providing a self-consistent magnetic boundary, the warm medium is dynamically uninteresting — if the sphere was instead placed in vacuo similar results would be expected, albeit with a slower initial collapse (which was inadvertently shown in Price & Bate, 2007, who accidentally performed their calculations with a pressure-free external medium).

We used a barotropic equation of state similar to Machida et al. (2008), without the final $\gamma=\nicefrac{{5}}{{3}}$ step,

[TABLE]

where the critical densities are given by $\rho_{\mathrm{c,1}}=10^{-14}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ and $\rho_{\mathrm{c,2}}=10^{-10}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ (the physical consequences of these choices are discussed in section 4).

The magnetic field is characterised by the dimensionless mass-to-flux ratio, $\mu$, and the angle between the rotation and field axes, $\vartheta$. The mass-to-flux ratio, in terms of the critical value, is defined as

[TABLE]

where the ratio between the magnetic and gravitational forces in a magnetised sphere is given by

[TABLE]

and the critical ratio, i.e. that at which the self-gravitational and magnetic pressure forces will be in equilibrium, is given by

[TABLE]

where the parameter $c_{1}=0.53$ (Mouschovias & Spitzer, 1976). The geometry of the initial field, $B_{0}$ is then defined according to

[TABLE]

[TABLE]

so that when $\vartheta=0\degree$ the field and rotation axes are parallel (i.e. aligned) and, conversely, when $\vartheta=90\degree$ they are perpendicular (i.e. totally misaligned). The rotation axis is aligned with the $z-$axis, so that the magnetic field axis then varies. Owing to the symmetry of the initial field – and the fact the field in the warm medium does not materially evolve over the course of the calculation – the field is also periodic across the boundaries, preventing numerical artifacts at the box edge.

We use mass-to-flux ratios of $\mu=[5,\leavevmode\nobreak\ 10,\leavevmode\nobreak\ 20]$, which give field strengths of $B_{0}=[163,\leavevmode\nobreak\ 81,\leavevmode\nobreak\ 41]$ microgauss and mean ratios of hydrodynamic to magnetic pressure, $\beta=[3.45,\leavevmode\nobreak\ 13.8,\leavevmode\nobreak\ 55.3]$ respectively (i.e. in all cases the simulation begins with the hydrodynamic pressure stronger than the magnetic pressure). For each field strength, we then perform the calculations 4 times, with values of $\vartheta=[0\degree,\leavevmode\nobreak\ 20\degree,\leavevmode\nobreak\ 45\degree,\leavevmode\nobreak\ 90\degree]$ to give a total of 12 models.

Finally, the sphere is placed in solid body rotation around the $z$-axis with $\Omega=1.77\times 10^{-13}\leavevmode\nobreak\ \text{rad}\thinspace\text{s}^{-1}$, which gives a ratio of rotational to gravitational energy of $\beta_{\mathrm{rot}}=0.005$, within the range of values observed by Goodman et al. (1993).

4 Results and Discussion

4.1 Isothermal Collapse

The sphere of cold gas described in section 3 is allowed to collapse. For this discussion, we divide the collapse of the cloud core into two phases: firstly an initial isothermal collapse, where the maximum densities are below $\rho=\rho_{\mathrm{c,1}}=10^{-14}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ and then a second phase where a disc structure and pre-stellar core is formed. We first describe the initial isothermal phase which, aside from slight variations in timings, is comparable for all three field strengths. We then, in sections 4.2, 4.3 and 4.4, consider the three field strengths in more detail, focusing in sections 4.2 and 4.3 on the formation of bipolar outflows and how these vary with the field geometry, and in section 4.4 on how the weakest field strength can form binary and multiple protostellar systems.

The first phase of the collapse extends from $t=0$ to $0.98t_{\mathrm{ff}}$, corresponding to the region in figure 1 where the maximum density increases approximately linearly with time (where the free–fall time, $t_{\mathrm{ff}}=24430\leavevmode\nobreak\ \text{yr}$). The main effect from the initial conditions observed here is that the lower mass–to–flux ratios provide higher magnetic pressures — $P_{\mathrm{mag}}\propto B^{2}$ so the doubling of the field strength between $\mu=10$ and $5$ provides a corresponding quadrupling of the pressure support. This results in the duration first phase of the collapse being slightly extended for the stronger fields. For example, the $\mu=10$ calculations reach a density of $10^{-10}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ $0.01\leavevmode\nobreak\ t_{\mathrm{ff}}$ quicker than the stronger field calculations, delaying the subsequent phase of evolution accordingly. Figure 2 shows the formation of the bipolar outflows, which are launched when a pseudo-disc forms and the density exceeds $\approx 10^{-10}\leavevmode\nobreak\ \text{g}\thinspace\text{cm}^{-3}$ (and consequently the $0.01\leavevmode\nobreak\ t_{\mathrm{ff}}$ delay when $\mu=5$) and which we discuss in sections 4.2 and 4.3.

The field geometry, i.e. $\vartheta$, does not affect the overall rate of the collapse since the magnetic pressure terms, which act to support the cloud against gravity, are isotropic and therefore independent of $\vartheta$. However, the tension component is anisotropic and, in general, acts to stop the fluid moving perpendicular to the magnetic field lines. Consequently this causes the collapsing core to become ellipsoidal. Thus, the $\vartheta=0\degree$ calculations produce an oblate spheroid while conversely the $\vartheta=90\degree$ calculations produce prolate spheroids and $20\degree$ and $45\degree$ produces a tri-axial ellipsoid. Figure 3 shows how the oblateness of the spheroids changes with field strength, with the axis ratios of $0.66$, $0.34$, and $0.24$ for $\mu=5$, $10$, and $20$ respectively. Equivalent values can be obtained for other geometries (although these are tri–axial ellipsoids or prolate spheroids).

Figure 1 shows how the maximum fluid density evolves as a function of time for all twelve calculations. At $0.98-0.99t_{\mathrm{ff}}$ (depending on the calculation) the density begins to rapidly increase. This marks the change from the isothermal and weakly oblate collapsing core to the second phase. The magnetic forces described above which acted to make the core oblate now cause the formation of a flattened disc-like structure. We use the term pseudo–disc to describe this object, which is not a ‘classical’ accretion disc in Keplerian rotation around a protostellar core. Nevertheless it is rotating albeit with a sub–Keplerian rotation profile, is supported by the gas pressure perpendicular to the pseudo-disc plane against further gravitational collapse parallel to the magnetic field axis, and is accreting onto the protostar (i.e. the radial velocity of the disc relative to the protostar is negative). We discuss this further in sections 4.2 and 4.3.

4.2 $\mu=5$

The strongest initial field strength used in this paper is $\mu=5$. After the initial isothermal collapse phase discussed above, all four field geometries form small pseudo-discs and bipolar outflows, which we show at $t=1.02t_{\mathrm{ff}}$ in figure 4.

Once a significant pseudo-disc forms, magnetic and rotational forces become comparable in magnitude to gravitational forces – driving the formation of effects like the bipolar outflows. Previous simulations of magnetised collapse and outflows used strong fields, aligned with the rotation axis (e.g. Tomisaka, 2002; Matsumoto & Tomisaka, 2004; Machida et al., 2008); here we also vary $\vartheta$ similar to the approach in Ciardi & Hennebelle (2010) and Lewis et al. (2015).

Inclination angles $\vartheta\leq 45\degree$ produce a bipolar outflow. It is convenient to divide these outflows into two regions: an overall bulk outflow and an inner collimated jet. We decompose the magnetic field into cylindrical co-ordinates (for example, as in Parker, 1955), such that the (magnitude of the) radial field is given by

[TABLE]

the magnitude of the poloidal field as

[TABLE]

with the azimuthal (i.e. toroidal) field given by

[TABLE]

Figure 5 shows the ratio of the poloidal to toroidal field in a cross–section slice for t the $\vartheta=0\degree$, $45\degree$, and $90\degree$ calculations. A thin (ca. $5\leavevmode\nobreak\ \text{au}$ wide) region dominated by the poloidal field may be observed, centred on the rotation axis and hence the jet. This poloidal field acts to move material away from the pseudo-disc which is dominated by the toroidal field. The disc and its toroidal field act to ‘wind-up’ the magnetic field and then material and angular momentum are ejected in the outflow; the central region of which is collimated by the toroidal component. A comparable effect, albeit with a somewhat lesser toroidal field contribution, is seen for $\vartheta=45\degree$ (central panel of figure 5).

In sharp contrast, the $\vartheta=90\degree$ calculation has no jet — but does have a small (with a maximum height above the disc of $10\leavevmode\nobreak\ \text{au}$) bulk outflow. This is produced by the small region around the pseudo-disc which is dominated by the poloidal magnetic field (shown in the right–hand panel of figure 5). We observe that above this, again centered on the rotation axis, a thin highly toroidal region exists centered on the rotation axis. However, this is not connected to the pseudo-disc and consequently can not act to collimate a jet — in essence, because the pseudo-disc is rotating in the same plane as the magnetic field lines it is unable to perform the same ‘winding’ effect seen when $\vartheta=0\degree$. The presence of a collimated bipolar outflow therefore clearly implies a combination of sufficent magnetic pressure (exemplified by the poloidal field magnitude) to lift material out of the pseudo–disc combined with a toroidal field which can collimate and drive the jet (cf. Lynden-Bell, 1996, 2001).

Earlier, we observed that changing the values of $\vartheta$ changed the collapsing cloud core from oblate to tri–axial. This effect continues into the later phase of the collapse where the consequent pseudo–disc aligns perpendicular to the magnetic field axis (and hence parallel to the major axis of the antecedent ellipsoid) (Galli & Shu, 1993). However, we find that the innermost region of the disc, with a radius of $<20\leavevmode\nobreak\ \text{au}$, attempts to re–align with the rotation axis, caused by the angular momentum at this scale being sufficient to overcome the magnetic tension force which is aligning the larger scale disc. This produces warped pseudo–discs for $\vartheta>0\degree$. This effect is most marked when $\vartheta=90\degree$ where the inner disc has a radius of $25\leavevmode\nobreak\ \text{au}$, but can also be seen for $20\degree$ and $45\degree$ albeit with a radii of $5$–$10\leavevmode\nobreak\ \text{au}$, and with a clearer connection to the rest of the pseudo–disc.

This inner region drives the central collimated region of the outflow. Although the opening angle of the bulk outflow is very large — the diameter of the whole outflow is on the order of $100\leavevmode\nobreak\ \text{au}$ — the collimated inner region is much smaller and is correlated with the inner disc. As a result the outflow initially aligns much closer to the rotation axis than the magnetic field axis, for example when $\vartheta=20\degree$, the inclination of the jet to the rotation axis is less than $10\degree$. The outflow (and in particular the jet) is removing angular momentum from the pseudo–disc so a close alignment between this and the direction of the angular momentum vector is to be expected.

4.3 $\mu=10$

The subsequent evolution of the $\mu=10$ calculations follows a similar overall pattern to $\mu=5$ calculations. Those with $\vartheta\leq 45\degree$ produce collimated outflows, which are faster and more closely aligned to the rotation axis for lower $\vartheta$. As we discussed earlier, however, the reduction in the magnetic field strength and the consequently reduced magnetic pressure support causes the core to collapse to higher densities more quickly. A more rapid collapse results in a correspondingly earlier formation of the outflow and jet; whilst the reduced magnetic braking produces different pseudo-discs and outflow velocities compare to $\mu=5$.

In the $\mu=5$ calculations we obtained a pseudo-disc with a radii of $\approx 50\leavevmode\nobreak\ \text{au}$ for $\vartheta<45\degree$ and $<15\leavevmode\nobreak\ \text{au}$ for $\vartheta\geq 45\degree$. Figure 6 (central columns) shows how for $\mu=10$ the same initial values of $\vartheta$ produce much larger pseudo–discs: we obtain radii of $\approx 100\leavevmode\nobreak\ \text{au}$ for $\vartheta\leq 45\degree$ and $25\leavevmode\nobreak\ \text{au}$ for $90\degree$. We see the same correlation between pseudo-disc radius and $\vartheta$, whereby the disc radii reduce as the initial inclination of the field axis to the rotation axis is increased, but for $\mu=10$ these are in general larger than the corresponding $\mu=5$ calculation. The apparently larger discs seen in the left column of figure 4, particularly for $\vartheta=0\degree$ compared to figure 6 is due to a projection effect caused by looking along the plane of the disc, which due to preferential collapse of the the fluid parallel to the field lines has a higher density than the rest of the collapsed core. Our calculation does not include a physical viscosity treatment, therefore the main process for transporting angular momentum — and hence allowing material to move towards the centre of the pseudo-disc — is magnetic braking. The calculations do include artificial viscosity, but this is limited using the Morris & Monaghan (1997) switch and should produce negligible angular momentum transport compared to that produced by the magnetic field. The degree of magnetic braking increases with the magnetic field strength, $|B^{i}|$, and $\mu=10$ implies an initial field which is half the magnitude of the $\mu=5$ calculations, causing the formation of the large discs. As we noted for $\mu=5$, the innermost region of the pseudo-disc (which drives the toroidal field and hence the collimated jet component of the outflow) realigns perpendicular to the rotation axis, although it is correspondingly larger for these calculations.

The most important difference between the two field strengths is that the outflows produced with $\mu=10$ are $60\leavevmode\nobreak\ \%$ faster. This effect is most pronounced for the $\vartheta=0\degree$ calculation. Like the corresponding calculation in section 4.2, this system initially had $|v_{z}|=5\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$ with the characteristic collimated jet. However, the inner jet region detaches and a region of faster flowing material with $|v_{z}|=8\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$ forms with the plasma $\beta$ above and below the pseudo–disc reducing from $10^{-2}$ to $10^{-3}$.

Figure 7 shows the ratio of the poloidal to toroidal field for this calculation. Whilst the $\mu=5$ calculation had a weak poloidal component which was then tightly would by the toroidal field from the disc; as the cross–section in figure 7 shows this calculation has a much stronger poloidal component, covering a larger region, above and below the disc. The disc is still dominated by the torodial ‘winding’ field, but the extra pressure contributed by the poloidal region acts to detach the collimated component of the jet and increase the overall speed. This effect, which manifests as an extra (magnetic) pressure contribution, causes the reduced plasma $\beta$ seen. A similar, although weaker, effect causes the partial disconnections seen when $20\degree\leq\vartheta\leq 45\degree$.

We defined $B_{p}$ in section 4.2 as essentially the magnetic pressure. The cause of this extra field is then simply the more rapidly rotating pseudo–disc, as shown in the radially average $v_{\phi}$ plot in figure 8, due to the initially reduced magnetic braking. In essence a larger disc and a faster outflow velocity are linked because the maximum jet velocity is itself linked to the maximum pseudo–disc rotational velocity, and this effect is continued across a range of misalignments. Additionally, because the initial field is weaker in these calculations, the $\vartheta=45\degree$ calculation — which for $\mu=5$ produced a notably weaker jet than for $\vartheta=20\degree$ — has a more comparable velocity and structure, and a similarly larger inner disc region.

4.4 $\mu=20$

Figure 9 shows the results from the $\mu=20$ calculations. These do not form bipolar outflows instead fragmenting into binary (or multiple) systems, as seen in the central column of figure 9. The precursor to this is the formation of a small (ca. $10\leavevmode\nobreak\ \text{au}$ radius) dense rotating disc: without the magnetic braking (and angular momentum loss from the outflow). These smaller discs, which are closer in nature to true accretion discs (cf. figure 10) than the pseudo-discs discussed earlier, are gravitationally unstable. This can be seen from the distinct spiral arm-like structures.

The $\vartheta=0\degree$ and $20\degree$ calculations form tight binaries with separations of $<5\leavevmode\nobreak\ \text{au}$. We add a cautionary note that in a calculation including additional physical processes, for example a radiative transfer scheme like that in Whitehouse & Bate (2004, 2006), these may not be formed due to the excess thermal support acting to stabilise the disc; alternatively they may merge together. In contrast, the separation between the two protostars for $\vartheta=90\degree$ is $\approx 13\leavevmode\nobreak\ \text{au}$ at $1.03t_{\mathrm{ff}}$. In all cases, the binary pairs are surrounded by a massive circumbinary disc, with a radius of $50\leavevmode\nobreak\ \text{au}$, while the wider binaries have two distinct dense circumstellar discs (with radii of $<5\leavevmode\nobreak\ \text{au}$) embedded within this.

All four values of $\vartheta$ produce an initial protostar by $t=1.01t_{\mathrm{ff}}$, which is comparable to the time frame for a $\mu=10$ calculation. However, these calculations then go on to produce a second protostar by $1.02t_{\mathrm{ff}}$, forming the binary pairs seen in figure 9. In addition, the $\vartheta=0\degree$ model eventually produces a further sink particle from the circumbinary material, shown in figure 12. The ultimate fate of this third potential protostar is unclear — whilst a stable triple star system is possible, by three-body interaction it could be ejected from the system or the two closer sinks could merge leaving a system similar to the $90\degree$ model albeit without the two discs). We note that this is a similar formation process to that proposed for the ternary system observed by Tobin et al. (2016), albeit with a much smaller separation between the binary pair and the younger third protostar. However, an important difference between all the $\mu=20$ models and the stronger fields is the absence of an outflow — without which angular momentum can not be rapidly removed from the disc. This reduction in angular momentum transport means that fragmentation into a binary (which can store more angular momentum) is the result.

Although we noted earlier that these four calculations do not produce a bipolar outflow, there are still regions with substantial magnetic fields. Consequently a significant magnetic pressure is still realised and this produces a region around the circumstellar (or circumbinary) disc with $\beta_{\mathrm{plasma}}<1$. This causes some material to be moved away from the plane of the disc, although not in the rapid and directed manner of a jet. We see, for example, in figure 11 that although a small poloidal region exists near the protostar, there is no collimating toroidal component present. This is more similar to the field structure seen for $\vartheta=90\degree$ in figures 5 and 7, which similarly have no bipolar outflow and no significant toroidal component near the protostar.

In principle, this may indicate that binary stars may be formed by fragmentation even in non–turbulent clouds, or clouds which have not been perturbed (e.g. by an external impulse as in Pringle, 1989) and are no longer axisymmetric; provided, of course, that the field is sufficiently weak so as to not drive a strong outflow. In practice, this may not be as strong a constraint as it seems — in our quasi–ideal MHD the only method of dissipating magnetic energy is by artificial resistivity, which is intentionally limited to the minimum required for stability. In reality, non–ideal effects (resistivity, ambipolar diffusion, and the Hall effect) will act to reduce or transform the field and may therefore allow an effect similar to that seen here to occur for stronger initial fields. Wurster et al. (2016) propose a consistent SPMHD method that extends the ideal MHD presented here to include all three non–ideal effects, and find that the Hall effect and the relative orientation of the rotation and magnetic field axes can have a significant effect on the consequent evolution of the core. Alternatively, as suggested in Machida et al. (2005) it may be sufficient to just increase the initial angular momentum of the cloud to promote fragmentation at higher field strengths (or some combination of both mechanisms).

Conversely, there should exist a range of sufficiently strong initial fields so that no matter what other physical processes are at play, binary or multiple star systems are difficult to form, and hence solitary stars may be produced. Stronger fields naturally drive stronger jets and outflows or increased magnetic braking which transport angular momentum obviating the need for companion stars to store angular momentum.

4.5 Accretion onto protostars

As well as producing outflows and discs with varying morphologies, each model also results in a different profile of accretion onto the sink particle (which represents a forming protostar). As before, a general trend is observed whereby the $\mu=20$ models and $\vartheta=90\degree$ models differ substantially from stronger fields and shallower angles. In figure 13 we plot maps of the material accreted during the 500 yr after a sink particle is inserted.

Figure 14 shows that the amount of material accreted onto the sink is broadly similar (at $t=1.04t_{\mathrm{ff}}$) across all of the $\mu=5$ and $10$ models, which is born out by the similar outflow morphology seen earlier in section 4. All calculations principally accrete material in the plane of their pseudo–disc, as seen in figure 13. When $\vartheta=0\degree$ this clearly aligns with the equator. Similarly, for weak magnetic fields, where $\mu=20$ (cf. figure 9), the accretion disc forms perpendicular to the rotation axis for all field geometries. The bottom row of figure 13 shows how this results in a strongly equatorial accretion profile. For $\vartheta=45\degree$ in the stronger field calculations a sinusoidal structure is seen in figure 13, which is the expected shape for a warped pseudo–disc like those seen in figures 4 and 6. The $\mu=5$, $\vartheta=90\degree$ has not produced a clear disc and has instead taken on a sigmoidal structure. This causes a hemispherically symmetrical and rotationally symmetrical accretion pattern corresponding to each ‘arm’ of the sigmoid. An imprint of a similar structure is seen for $\mu=10$ from the two spiral arms of the disc. However, unlike in Lewis et al. (2015)  we find that the change to the small scale structure results in a roughly equivalent accretion rate for both $\vartheta=45\degree$ and $90\degree$.

4.6 Comparison to observations

Observations of young stellar objects (YSOs) substantially more evolved than a first core (Larson, 1969) are plentiful compared to detections of objects which, although they have formed some sort of core, are much less evolved. Detections of an earlier phase — a dense starless core which is approaching the point of forming a star — are better represented, e.g. those found by Crapsi et al. (2005). Therefore, constraints on the dynamics a first hydrostatic core and its environs would be useful. The insertion of sink particles means we do not follow the collapse to exactly first hydrostatic core proposed by Larson (1969). However, our sink particle has comparable dimensions so the overall structure of the pseudo–disc and any outflow should be similar.

In models where we obtain a collimated (or jet–like) outflow we find that the bulk velocity ranges between $\approx 2$ and $\approx 8\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$ depending on the exact configuration of parameters, with higher values of $\vartheta$ corresponding to generally reduced outflow velocities. This substantially increases the range of potential outflow velocities which may be first cores as opposed to ‘very low luminosity objects’ (VeLLOs, see Myers et al., 2005; Kauffmann et al., 2005). One example is given in Pineda et al. (2011), who find that L1451-mm has a slow outflow with a velocity of about $3\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$. Similarly, Dunham et al. (2011) observe an bipolar outflow of comparable velocity in Per–Bolo 58. These velocities are apparently slightly too slow if we only consider lower values of $\vartheta$, however, an outflow of this velocity would be approximately consistent with our results for $\vartheta=45\degree$.

At the other extreme, Chen et al. (2010) find that L1448 IRS2E (another first core candidate) has a clearly collimated outflow with a velocity of around $25\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$. This is much faster than any of the models presented here, which have an artificial limit of $1\leavevmode\nobreak\ \text{au}$ placed on the collapse. Outflow velocity and the maximum rotation speed (which is controlled by the larger of the size of the protostar or sink particle) are closely linked (Price et al., 2003) and consequently this indicates that whatever object is present in L1448 IRS2E, it has a substantially smaller radius than $1\leavevmode\nobreak\ \text{au}$, and may therefore be more evolved than a first core.

In effect we find that slow bipolar outflows are characteristic of first cores at all but the weakest field strengths. This implies that there is a strong constraint on first core candidates: any outflow, if present, should be between $\approx 2$ and $8\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$ and that faster outflows may be indicative of a more evolved object.

5 Conclusions

Using smoothed particle magnetohydrodynamics we have performed twelve simulations (using three mass–to–flux ratios, and four angles between the field and rotation axes) of the collapse of molecular cloud cores to form protostars. Whilst we observe that the initial collapse of the cloud core is broadly unaffected by the choice of parameters, the subsequent evolution can be quite different. All four $\mu=20$ calculations produce binary systems of some form, with the degree of binary separation strongly dependent on the initial field geometry. At shallow angles, i.e. low values of $\vartheta$ the two protostars form very close together, within $\lesssim 5\text{\leavevmode\nobreak\ au}$; at steeper angles two separate discs form orbiting a common barycentre. This may indicate that binary or multiple star systems may be formed in molecular clouds if the initial magnetic field is very weak. We also find that at this field strength significant first core outflows are absent: instead the system conserves angular momentum by forming a binary.

In contrast, no simulation with $\mu\leq 10$ formed a binary for any field geometry. However, the structure and nature of the resultant pseudo–disc and outflow varies between both values of $\mu$ and $\vartheta$. In addition, the $\mu=10$ models are observed to collapse somewhat more quickly than those with a stronger field. We find that for all values of $\vartheta$, the pseudo–disc formed for $\mu=10$ is significantly larger than for $\mu=5$ due to reduced magnetic braking. In particular, when $\vartheta=90\degree$ the stronger field fails to produce a disc at all, compared to the weaker field where a comparatively large disc. Additionally, the outflows produced may have different speeds, ranging from $2$ to $8\leavevmode\nobreak\ \text{km}\thinspace\text{s}^{-1}$ depending on $\mu$, $\vartheta$, and time.

Additionally, changing the initial field geometry from a configuration where the rotation and field axis are aligned to one where they are misaligned can promote the formation of warped pseudo–discs. Here the outer regions of the discs orient perpendicular to the magnetic field axis (due to reduced magnetic pressure support parallel to the field lines) but the inner region begins to re–orient perpendicular to the rotation axis as angular momentum is concentrated near the sink particle.

These results may have consequences for how binary or single star systems are formed. Clearly, a weaker field, with the attendant lack of both magnetic braking and outflow, is liable to fragment into a binary. Similarly, at the strongest field strengths (with aligned fields), the effect of the magnetic braking plus a large collimated outflow may be sufficient to prevent this fragmentation from occurring. However, a large intermediate regime exists which may allow the formation of binary or single star systems depending on the exact parameters of the initial cloud core; although adding additional physics, for example radiative transfer, may affect this. Further work, following the collapse to a protostellar core itself would be necessary to determine the exact boundary between the two regimes.

Acknowledgments

The authors thank D. J. Price for his comments on the draft manuscript.

BTL acknowledges support from an STFC Studentship and Long Term Attachment grant. This work was also supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 339248). MRB’s visit to Monash was funded by an International Collaboration Award from the Australian Research Council (ARC) under the Discovery Project scheme grant DP130102078.

This work used 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. DiRAC is part of the National E-Infrastructure.

Calculations were also performed on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter.

Rendered plots were produced using the splash (Price, 2007) visualisation programme.