The impact of non-ideal magnetohydrodynamics on binary star formation
MMNRAS , 1–17 (2016) Preprint 9 November 2018 Compiled using MNRAS L A TEX style file v3.0
The impact of non-ideal magnetohydrodynamics on binary starformation
James Wurster , (cid:63) , Daniel J. Price and Matthew R. Bate † School of Physics, University of Exeter, Stocker Rd, Exeter EX4 4QL, UK Monash Centre for Astrophysics and School of Physics and Astronomy, Monash University, Vic 3800, Australia
Submitted: Revised: Accepted:
ABSTRACT
We investigate the effect of non-ideal magnetohydrodynamics (MHD) on the formation of bi-nary stars using a suite of three-dimensional smoothed particle magnetohydrodynamics sim-ulations of the gravitational collapse of one solar mass, rotating, perturbed molecular cloudcores. Alongside the role of Ohmic resistivity, ambipolar diffusion and the Hall effect, wealso examine the effects of magnetic field strength, orientation and amplitude of the densityperturbation. When modelling sub-critical cores, ideal MHD models do not collapse whereasnon-ideal MHD models collapse to form single protostars. In super-critical ideal MHD mod-els, increasing the magnetic field strength or decreasing the initial density perturbation ampli-tude decreases the initial binary separation. Strong magnetic fields initially perpendicular tothe rotation axis suppress the formation of binaries and yield discs with magnetic fields ∼ times stronger than if the magnetic field was initially aligned with the rotation axis. Whennon-ideal MHD is included, the resulting discs are larger and more massive, and the binaryforms on a wider orbit. Small differences in the super-critical cores caused by non-ideal MHDeffects are amplified by the binary interaction near periastron. Overall, the non-ideal effectshave only a small impact on binary formation and early evolution, with the initial conditionsplaying the dominant role. Key words: methods: numerical — magnetic fields — MHD — stars: binaries: general —stars: formation
Most observed stars are members of multiple systems (e.g.Duquennoy & Mayor 1991; Raghavan et al. 2010), so the forma-tion of binaries is crucial to our understanding of the star formationprocess. The most common method for simulating binary formationnumerically to date has been to model the gravitational collapse ofa rotating, uniform sphere given an initial m = 2 density perturba-tion either without (e.g. Boss & Bodenheimer 1979; Boss 1986;Burkert & Bodenheimer 1993; Bate et al. 1995; Machida et al.2008b) or with (e.g. Price & Bate 2007; Hennebelle & Teyssier2008; Boss & Keiser 2013) magnetic fields.Using this approach, Price & Bate (2007) (hereafter PB07)found that strong magnetic fields (mass-to-flux ratios (cid:46) timescritical) could suppress binary formation — leading to the forma-tion of only one star where two would otherwise have formed —influenced more by magnetic pressure than by magnetic tension ormagnetic braking effects. The effect was more pronounced withmagnetic fields initially perpendicular to the rotation axis. A sim- (cid:63) [email protected] † [email protected] ilar study by Hennebelle & Teyssier (2008) considered magneticfields parallel to the axis of rotation and an initial rotation that was ∼ times slower than in PB07. They also found that magnetic fieldsgreatly inhibited binary formation, finding disc fragmentation andbinary formation only for weak initial magnetic fields (mass-to-fluxratios (cid:38) times critical). When the initial density perturbationwas increased, then binaries formed except when the initial cloudwas close to being magnetically critical. They concluded that thesuppression of binary formation was a result of the growth of thetoroidal component of the magnetic field, and not a result of mag-netic braking. Further, large density perturbations were required tofrom binaries given the observed magnetic field strengths. c (cid:13) a r X i v : . [ a s t r o - ph . S R ] D ec Wurster, Price & Bate
However, these studies assumed ideal magnetohydrodynam-ics (MHD), which is a poor approximation for molecular clouds(Mestel & Spitzer 1956), where ionisation fractions are of order n e /n H ∼ − (Nakano & Umebayashi 1986a; Umebayashi &Nakano 1990). Partial ionisation means that non-ideal MHD effects— specifically Ohmic resistivity, the Hall effect, and ambipolar dif-fusion — become important, with the relative importance of eachdepending, amongst other things, on the gas density and magneticfield strength (e.g. Wardle & Ng 1999; Nakano et al. 2002; Tassis &Mouschovias 2007a; Wardle 2007; Pandey & Wardle 2008; Keith& Wardle 2014). The Hall effect also depends on the direction ofthe magnetic field with respect to the axis of rotation (e.g. Braiding& Wardle 2012).Previous studies have examined the effects of non-ideal MHDon the formation of single stars (e.g. Nakano & Umebayashi 1986b;Fiedler & Mouschovias 1993; Ciolek & Mouschovias 1994; Li &Shu 1996; Mouschovias 1996; Mouschovias & Ciolek 1999; Shuet al. 2006; Mellon & Li 2009; Duffin & Pudritz 2009; Dapp &Basu 2010; Machida et al. 2011; Li et al. 2011; Dapp et al. 2012;Tomida et al. 2013, 2015; Tsukamoto et al. 2015a,b; Wurster et al.2016). Tsukamoto et al. (2015b) and Wurster et al. (2016) (here-after WPB16) found that, although the Hall effect was not the dom-inant non-ideal effect in numerical simulations of an isolated form-ing star, it was the controlling factor in disc formation, with a largedisc forming when the initial magnetic field and rotation vectorwere anti-aligned but no disc forming when the initial magneticfield direction was reversed.Here, we evaluate the influence of non-ideal MHD on the for-mation and early evolution of binary systems, following the origi-nal ideal MHD studies of PB07 and Hennebelle & Teyssier (2008).We model 3D non-ideal self-gravitating smoothed particle magne-tohydrodynamics simulations of collapsing, low mass cores, withthe ionisation fractions calculated using the N ICIL library (Wurster2016). We present the numerical formulation in Section 2, the ini-tial conditions in Section 3, the results in Section 4 and the discus-sion and conclusions in Section 5.
We solve the equations of self-gravitating, non-ideal magnetohy-drodynamics in the form d ρ d t = − ρ ∇ · v , (1) d v dt = − ρ ∇ (cid:20)(cid:18) P + 12 B (cid:19) I − BB (cid:21) − ∇ Φ , (2) d B d t = ( B · ∇ ) v − B ( ∇ · v ) + d B d t (cid:12)(cid:12)(cid:12)(cid:12) non-ideal , (3) ∇ Φ = 4 πGρ, (4)where dd t ≡ ∂∂t + v · ∇ is the Lagrangian derivative, ρ is the density, v is the velocity, P is the gas pressure, B is the magnetic field, Φ isthe gravitational potential, I is the identity matrix, and d B d t (cid:12)(cid:12) non-ideal is the non-ideal MHD term which is a sum of the Ohmic resistivity (OR), Hall effect (HE) and ambipolar diffusion (AD) terms,d B d t (cid:12)(cid:12)(cid:12)(cid:12) OR = − ∇ × [ η OR ( ∇ × B )] , (5)d B d t (cid:12)(cid:12)(cid:12)(cid:12) HE = − ∇ × (cid:104) η HE ( ∇ × B ) × ˆ B (cid:105) , (6)d B d t (cid:12)(cid:12)(cid:12)(cid:12) AD = ∇ × (cid:110) η AD (cid:104) ( ∇ × B ) × ˆ B (cid:105) × ˆ B (cid:111) , (7)where η OR , η HE and η AD are the non-ideal MHD coefficients. Weassume units for the magnetic field such that the Alfv´en speed is v A ≡ B/ √ ρ (see Price & Monaghan 2004).We close the equation set using a barotropic equation of state, P = c s,0 ρ ; ρ < ρ c ,c s,0 ρ c ( ρ/ρ c ) / ; ρ c ≤ ρ < ρ d ,c s,0 ρ c ( ρ d /ρ c ) / ( ρ/ρ d ) / ; ρ ≥ ρ d , (8)where c s,0 is the initial isothermal sound speed and ρ c = 10 − and ρ d = 10 − g cm − . These density thresholds are the same as inPrice et al. (2012), Lewis et al. (2015) and WPB16. Although we donot employ full radiation magnetohydrodynamics, the barotropicequation of state is designed to mimic the evolution of the equationof state in molecular clouds (Larson 1969; Masunaga & Inutsuka2000; Machida et al. 2008a). Our value of ρ c is one order of magni-tude lower than the typically chosen value, and is chosen to reducefragmentation in the gas surrounding the protostars in the absenceof radiation feedback.We use Version 1.1 of the N ICIL library (Wurster 2016) tocalculate the non-ideal MHD diffusion coefficients. WPB16 useda precursor to N
ICIL , and the main difference is that this libraryincludes thermal ionisation and more detailed cosmic ray ionisa-tion. The thermal ionisation processes can singly ionise hydrogen,and doubly ionise helium, sodium, magnesium and potassium; themass fractions of the five elements are . , . , . × − , . × − and . × − , respectively (e.g. Asplund et al. 2009;Keith & Wardle 2014). Due to the cool temperatures in this study,we do not expect thermal ionisation to significantly contribute tothe electron population. Cosmic rays have the ability to remove anelectron to create an ion, which may be absorbed by a dust grain.We assume that two species of ions can be created: a heavy ionrepresented by magnesium (Asplund et al. 2009) and a light ionrepresenting hydrogen and helium compounds whose mass is cal-culated from the hydrogen and helium mass fractions (in WPB16we considered only the heavy ion species). For most of the calcula-tions, as in WPB16, we model a single grain species that can absorbthe electrons, and further assume that these grains have radius, bulkdensity, and average electric charge of a g = 0 . µ m, ρ b = 3 g cm − (Pollack et al. 1994) and ¯ Z g < , respectively. Our simulations are performed using the 3D smoothed particlemagnetohydrodynamics (SPMHD) code P
HANTOM (Price & Fed-errath 2010; Lodato & Price 2010), which includes self-gravity.The discretised magnetohydrodynamic equations (see review byPrice 2012) are given in Wurster et al. (2014) and WPB16. Weemploy the constrained hyperbolic divergence cleaning algorithmdescribed by Tricco & Price (2012) and Tricco et al. (2016) to con-trol divergence errors in the magnetic field.We adopt the usual cubic spline kernel, set such that the ratioof the smoothing length to the particle spacing is equivalent to ∼ MNRAS000
HANTOM (Price & Fed-errath 2010; Lodato & Price 2010), which includes self-gravity.The discretised magnetohydrodynamic equations (see review byPrice 2012) are given in Wurster et al. (2014) and WPB16. Weemploy the constrained hyperbolic divergence cleaning algorithmdescribed by Tricco & Price (2012) and Tricco et al. (2016) to con-trol divergence errors in the magnetic field.We adopt the usual cubic spline kernel, set such that the ratioof the smoothing length to the particle spacing is equivalent to ∼ MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation neighbours in three dimensions (Price 2012). We solve Poisson’sequation, ∇ Φ a = 4 πGρ a , following Price & Monaghan (2007) atshort range, and use a k -d tree algorithm similar to that described inGafton & Rosswog (2011) to compute the long range gravitationalinteraction in an efficient manner.Finally, the non-ideal MHD timestep is constrained by d t a 75 ( B = 1090 µ G ) ,slower initial rotations, and the effect of using multiple grain pop-ulations in the non-ideal MHD algorithm.We evolve the simulations that form binaries until at least firstapoastron. In the following, we refer to gas densities ρ > ρ disc,min =10 − g cm − centred on a sink particle as a ‘disc’, although theseare not necessarily long-lived or rotationally supported. The radiusof the disc is defined as where the density drops to 10 per cent ofthe maximum density and the disc mass is defined as the gas masswith ρ > ρ disc,min enclosed within this radius. The star+disc mass isthe sum of the disc and sink particle masses. This definition differsfrom the one used in PB07 and WPB16 to avoid calculating artifi-cially large radii and masses due to tidal bridges between two sinkparticles.When a sink particle is first formed, it represents a first hydro-static core, which exists for − yr (e.g. Tomida et al. 2010;Bate 2011) before collapsing to a stellar core (Larson 1969). Giventhat we follow the binary for . - × yr after the formation ofthe first sink particles, we refer to sink particles as protostars.Our full set of models is summarised in Table A1 of Ap-pendix A, listing the time of first periastron t peri , the initial period MNRAS , 1–17 (2016) Wurster, Price & Bate T calculated using first periastron and first apoastron, and the sep-arations at first periastron R peri and first apoastron R apo . Our first set of simulations assume ideal MHD for comparison withprevious studies. Given that the evolution of the magnetic field inideal MHD is independent of the sign of B , we consider only twoinitial magnetic field geometries: B = − B z and + B x . Fig. 1shows the evolution of the gas column density for the ideal MHDmodel with µ = 5 , A = 0 . and B = − B z , which is represen-tative of our suite of models.As the perturbed cloud collapses, each of the two over-densities collapse into a protostar. Their first periastron and apoas-tron are and au, respectively. Over the first seven periods,the mean periastron and apoastron are and au, with an av-erage period of . t ff , indicating that the binary is dynamicallyevolving. Fig. 2 shows the evolution of the binary separation (topright panel, green curve).As found by PB07 and Hennebelle & Teyssier (2008), boththe strength of the initial perturbation and the mass-to-flux ratio in-fluence the formation of the binaries. Fig. 3 shows the snapshots forour models at t = 1 . and 1.34 t ff . For smaller initial density per-turbations, the time of first periastron, the first periastron and firstapoastron separations all decrease. As the magnetic field strength isdecreased, R peri and R apo both increase since there is less magneticbraking.Not all models yield stable orbits as shown in Fig. 1. Themodel with µ = 10 and A = 0 . , for example, forms two mas-sive, dense discs which fragment prior to first periastron. Each pairof protostars forms a tight binary, and these pairs orbit one anotheron a long period; during their formation, the discs are totally dis-rupted. Near apoastron, they interact with younger protostars whichdisrupts the orbit and makes the system chaotic. Thus, this model,and many other models that have more than two protostars yieldinteractions that hinder a useful comparison.The green lines in Fig. 2 show the evolution of the disc andstar+disc masses, the disc radius, and the mass-weighted plasma β and magnetic field in the disc around one protostar for the idealMHD models with A = 0 . and µ = 10 (left-hand column)and µ = 5 (right-hand column). Over the seven periods, startingat first periastron, the disc radius for the µ = 5 model varies be-tween and au, and its mass varies between . and . M (cid:12) ,where the local minima corresponds to periastron. After the initialrapid growth of the protostar, its subsequent growth is not depen-dent on orbital position, and the fluctuations in the star+disc masscorrespond to fluctuations in the disc mass. The disc is always dom-inated by gas pressure rather than magnetic pressure, with β (cid:38) .As expected from the symmetry of the model, the properties aroundboth protostars follow the same trends until one or both of the discsfragments.As expected, the magnetic field strength is higher in the discsof the µ = 5 models compared to the µ = 10 models. How-ever, they are lower than in the discs produced during the collapseto form an isolated protostar (WPB16), and hence have larger val-ues of plasma β . Thus, magnetic fields are less important in theevolution of the discs in these binary models than in the isolatedprotostar models of WPB16.Weaker magnetic fields produce discs at larger separations,which have a larger gas reservoir than their strong field counter-parts. As shown in WPB16, larger discs form in weaker magnetic fields due to less efficient magnetic braking, independent of the gasreservoir. Thus, these two complementary effects result in largerdiscs in weaker magnetic fields. The largest discs form in the modelwith µ = 10 and A = 0 . while the smallest discs form in themodel with µ = 5 and A = 0 . . We would reach the same con-clusions if the models were instead compared exactly at the time offirst periastron.The weak field model exists on a long-period orbit, and doesnot have a second periastron by the end of the simulation (top leftpanel of Fig. 2, green curve). This allows the discs to essentiallyevolve in isolation, with the radius reaching a steady size of r ∼ au, even though the mass is continually decreasing. The sharp andperiodic decreases observed in the µ = 5 models do not occur. Fig. 4 shows a repeat of the above calculations using an initial mag-netic field perpendicular to the rotation axis (i.e. B = + B x ); aspreviously shown, the gas column density is at t = 1 . and 1.34 t ff .The green lines in Fig. 5 show the evolution of the separation of thetwo protostars and the evolution of the properties of one disc for themodel with µ = 10 and A = 0 . . As with the B = − B z mod-els, after the initial growth of the protostar, the fluctuations in thestar+disc mass correspond to the fluctuations in the disc mass, andthe disc is always supported by gas pressure rather than magneticpressure.For weak magnetic fields ( µ = 10 ; upper subfigure in Fig. 4)initially perpendicular to the axis of rotation, the separation of thebinary at first periastron is larger than its B = − B z counterpart,resulting in less interaction and a shorter period. The correspond-ing disc radii and masses are smaller in the B = + B x models,discussed further in Section 4.3 below.Strong magnetic fields ( µ = 5 ; lower subfigure in Fig. 4)initially perpendicular to the axis of rotation suppress the forma-tion of binaries. For A = 0 . , a binary forms early with first pe-riastron occurring at t ≈ . t ff , compared to t ≈ . t ff for its B = − B z counterpart. The apoastron distance is R apo < au,whereas this is a typical periastron distance for its − B z counter-part. The A = 0 . model forms a binary pair with a semi-majoraxis of 3 au and a common disc, while the A = 0 . model formsa single protostar and disc. For the purposes of our analysis, the A = 0 . model is treated as single protostar.The magnetic field strength is larger and the plasma β issmaller in the discs of the B = + B x models than their B = − B z counterparts, indicating that the magnetic field is more im-portant in the B = + B x models for the evolution of the disc. Comparison of Fig. 3 to Fig. 4 demonstrates that the evolution ofthe magnetic field depends on its initial orientation. At t = 1 . t ff ,the net magnetic field in the discs of the µ = 10 models are ∼ times higher in the B = + B x models than their − B z counterpartswith the same A , despite the same initial strength. Since strongermagnetic fields enhance magnetic braking, the discs in the B =+ B x models are smaller and less massive.To quantify this, Fig. 6 compares the magnetic field strengthin the mid-plane ( z = 0 ) for the ideal MHD models at t = 1 . t ff .The magnetic field strength of the B = + B x models is higherthroughout the mid-plane and in the discs than in their respec-tive B = − B z models. As the vertical collapse proceeds in the MNRAS000 Wurster, Price & Bate T calculated using first periastron and first apoastron, and the sep-arations at first periastron R peri and first apoastron R apo . Our first set of simulations assume ideal MHD for comparison withprevious studies. Given that the evolution of the magnetic field inideal MHD is independent of the sign of B , we consider only twoinitial magnetic field geometries: B = − B z and + B x . Fig. 1shows the evolution of the gas column density for the ideal MHDmodel with µ = 5 , A = 0 . and B = − B z , which is represen-tative of our suite of models.As the perturbed cloud collapses, each of the two over-densities collapse into a protostar. Their first periastron and apoas-tron are and au, respectively. Over the first seven periods,the mean periastron and apoastron are and au, with an av-erage period of . t ff , indicating that the binary is dynamicallyevolving. Fig. 2 shows the evolution of the binary separation (topright panel, green curve).As found by PB07 and Hennebelle & Teyssier (2008), boththe strength of the initial perturbation and the mass-to-flux ratio in-fluence the formation of the binaries. Fig. 3 shows the snapshots forour models at t = 1 . and 1.34 t ff . For smaller initial density per-turbations, the time of first periastron, the first periastron and firstapoastron separations all decrease. As the magnetic field strength isdecreased, R peri and R apo both increase since there is less magneticbraking.Not all models yield stable orbits as shown in Fig. 1. Themodel with µ = 10 and A = 0 . , for example, forms two mas-sive, dense discs which fragment prior to first periastron. Each pairof protostars forms a tight binary, and these pairs orbit one anotheron a long period; during their formation, the discs are totally dis-rupted. Near apoastron, they interact with younger protostars whichdisrupts the orbit and makes the system chaotic. Thus, this model,and many other models that have more than two protostars yieldinteractions that hinder a useful comparison.The green lines in Fig. 2 show the evolution of the disc andstar+disc masses, the disc radius, and the mass-weighted plasma β and magnetic field in the disc around one protostar for the idealMHD models with A = 0 . and µ = 10 (left-hand column)and µ = 5 (right-hand column). Over the seven periods, startingat first periastron, the disc radius for the µ = 5 model varies be-tween and au, and its mass varies between . and . M (cid:12) ,where the local minima corresponds to periastron. After the initialrapid growth of the protostar, its subsequent growth is not depen-dent on orbital position, and the fluctuations in the star+disc masscorrespond to fluctuations in the disc mass. The disc is always dom-inated by gas pressure rather than magnetic pressure, with β (cid:38) .As expected from the symmetry of the model, the properties aroundboth protostars follow the same trends until one or both of the discsfragments.As expected, the magnetic field strength is higher in the discsof the µ = 5 models compared to the µ = 10 models. How-ever, they are lower than in the discs produced during the collapseto form an isolated protostar (WPB16), and hence have larger val-ues of plasma β . Thus, magnetic fields are less important in theevolution of the discs in these binary models than in the isolatedprotostar models of WPB16.Weaker magnetic fields produce discs at larger separations,which have a larger gas reservoir than their strong field counter-parts. As shown in WPB16, larger discs form in weaker magnetic fields due to less efficient magnetic braking, independent of the gasreservoir. Thus, these two complementary effects result in largerdiscs in weaker magnetic fields. The largest discs form in the modelwith µ = 10 and A = 0 . while the smallest discs form in themodel with µ = 5 and A = 0 . . We would reach the same con-clusions if the models were instead compared exactly at the time offirst periastron.The weak field model exists on a long-period orbit, and doesnot have a second periastron by the end of the simulation (top leftpanel of Fig. 2, green curve). This allows the discs to essentiallyevolve in isolation, with the radius reaching a steady size of r ∼ au, even though the mass is continually decreasing. The sharp andperiodic decreases observed in the µ = 5 models do not occur. Fig. 4 shows a repeat of the above calculations using an initial mag-netic field perpendicular to the rotation axis (i.e. B = + B x ); aspreviously shown, the gas column density is at t = 1 . and 1.34 t ff .The green lines in Fig. 5 show the evolution of the separation of thetwo protostars and the evolution of the properties of one disc for themodel with µ = 10 and A = 0 . . As with the B = − B z mod-els, after the initial growth of the protostar, the fluctuations in thestar+disc mass correspond to the fluctuations in the disc mass, andthe disc is always supported by gas pressure rather than magneticpressure.For weak magnetic fields ( µ = 10 ; upper subfigure in Fig. 4)initially perpendicular to the axis of rotation, the separation of thebinary at first periastron is larger than its B = − B z counterpart,resulting in less interaction and a shorter period. The correspond-ing disc radii and masses are smaller in the B = + B x models,discussed further in Section 4.3 below.Strong magnetic fields ( µ = 5 ; lower subfigure in Fig. 4)initially perpendicular to the axis of rotation suppress the forma-tion of binaries. For A = 0 . , a binary forms early with first pe-riastron occurring at t ≈ . t ff , compared to t ≈ . t ff for its B = − B z counterpart. The apoastron distance is R apo < au,whereas this is a typical periastron distance for its − B z counter-part. The A = 0 . model forms a binary pair with a semi-majoraxis of 3 au and a common disc, while the A = 0 . model formsa single protostar and disc. For the purposes of our analysis, the A = 0 . model is treated as single protostar.The magnetic field strength is larger and the plasma β issmaller in the discs of the B = + B x models than their B = − B z counterparts, indicating that the magnetic field is more im-portant in the B = + B x models for the evolution of the disc. Comparison of Fig. 3 to Fig. 4 demonstrates that the evolution ofthe magnetic field depends on its initial orientation. At t = 1 . t ff ,the net magnetic field in the discs of the µ = 10 models are ∼ times higher in the B = + B x models than their − B z counterpartswith the same A , despite the same initial strength. Since strongermagnetic fields enhance magnetic braking, the discs in the B =+ B x models are smaller and less massive.To quantify this, Fig. 6 compares the magnetic field strengthin the mid-plane ( z = 0 ) for the ideal MHD models at t = 1 . t ff .The magnetic field strength of the B = + B x models is higherthroughout the mid-plane and in the discs than in their respec-tive B = − B z models. As the vertical collapse proceeds in the MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation t=1.26t ff Ideal, μ =5B =-B z , A =0.1 t=1.3t ff t=1.34t ff t=1.38t ff t=1.42t ff t=1.46t ff t=1.5t ff t=1.53t ff t=1.57t ff t=1.61t ff t=1.65t ff t=1.69t ff t=1.73t ff t=1.77t ff t=1.81t ff t=1.85t ff t=1.89t ff t=1.93t ff t=1.97t ff t=2.01t ff t=2.05t ff t=2.09t ff t=2.13t ff t=2.17t ff t=2.21t ff t=2.25t ff t=2.29t ff t=2.33t ff t=2.37t ff t=2.41t ff t=2.45t ff t=2.49t ff t=2.53t ff t=2.57t ff l o g c o l u m n d e n s i t y [ g / c m ] t=2.61t ff 400 AU Figure 1. Evolution of the face-on gas column density of the ideal MHD model with µ = 5 , A = 0 . and B = − B z , where µ is the mass-to-flux ratioin units of its critical value and A is the amplitude of the initial m = 2 density perturbation. The frames are at intervals of d t = 0 . t ff , where the free-falltime is t ff = 2 . × yr. Each frame is (1200 AU) , and the grey circles of radius 200 au are included for reference. The small solid white circles representsink particles with the radius of the circle representing the accretion radius of the sink particle. The binary is on a stable elliptical orbit, with first apoastron at R apo = 440 au at t = 1 . t ff . Over the first seven periods, the mean periastron and apoastron are and au, with an average period of d t ≈ . t ff .MNRAS , 1–17 (2016) Wurster, Price & Bate S e p a r a ti on ( a u ) µ =10, A =0.1 Ideal, B =-B z Non-ideal, B =-B z Non-ideal, B =+B z µ =5, A =0.1 0.0 0.1 0.2 0.3 0.4 0.5 M a ss ( M s un ) discstar+disc 0 20 40 60 80 100 120 R d i s c ( a u ) -1 p l a s m a β -4 -3 -2 -1 B ( G ) Time (t ff ) 1.4 1.6 1.8 2.0 2.2 2.4 2.6Time (t ff ) Figure 2. Time evolution of selected models starting from the formation of the sink particles (first representing first hydrostatic cores which collapse toprotostars; hereafter referred to as protostars). Each model is initialised with an initial density perturbation of A = 0 . and B = ± B z . The left- andright-hand column shows three models with µ = 10 and , respectively. Top to bottom : the separation of the two protostars, the disc and star+disc masses,the disc radius, and the mass-weighted plasma β and magnetic field strength in the disc around one protostar. The vertical line at t = 1 . t ff corresponds tothe time of our early analysis. The radius and mass of the disc in the non-ideal MHD models is for the non-fragmented disc. For the ideal MHD model with µ = 5 , the mean period, periastron and apoastron are . t ff , and au, respectively. The oscillations in the non-ideal MHD models with µ = 5 after second periastron are epicycles which are a result of one disc fragmenting and forming a well-behaved tight binary; the plotted binary separation is ofthe two initial protostars, and not to the barycentre of the newly formed tight binary. The local minima and decreases in disc radii correspond to periastron.The protostar continues to accrete mass as the models evolve, while the mass of the discs generally decrease. The increasing separation after t ∼ . t ff inthe non-ideal models with µ = 10 is a result of the primary protostars interacting with younger protostars that modify the orbit of the primaries. The linesterminate at the end of the simulation. B = + B x models, the field is dragged into the mid-plane. Whenthe discs form, the stronger magnetic field is wound into the disc,further enhancing its strength; in the B = + B x models, theazimuthal magnetic field, B φ , is the dominant component. In the B = − B z models, the radial dragging of the magnetic fields en-hances the field strength in the discs compared to the background,but not compared to their B x counterparts.In all models, there is little conversion of horizontal magneticfields into vertical fields or vice versa, hence only a weak verti-cal (horizontal) magnetic field develops in the B = + B x ( − B z )models. For example, on average at t = 1 . t ff , the φ -componentof the magnetic field is ∼ . times stronger than the z -componentfor the model with µ = 5 , A = 0 . and B = + B x , whilethe z -component is . times stronger than the φ -component in itscounterpart model with B = − B x . Previous studies have demonstrated that non-ideal MHD affects theformation and evolution of discs around protostars forming in iso-lation (e.g. Li et al. 2011; Tomida et al. 2015; Tsukamoto et al.2015b,a; WPB16). Further, when the Hall effect is included, thedirection of the magnetic field with respect to the axis of rotationaffects the evolution, with larger discs forming for cases where themagnetic field is anti-aligned with the axis of rotation (Braiding &Wardle 2012; Tsukamoto et al. 2015b; WPB16).We thus perform a suite of non-ideal MHD models with thesame parameters as the ideal MHD models discussed above, exceptthat we also run models where the sign of the magnetic field isreversed. Fig. 7 shows the gas column density at t = 1 . t ff of the B = ± B z models. The effect of non-ideal MHD on the results (compar- MNRAS000 Fig. 7 shows the gas column density at t = 1 . t ff of the B = ± B z models. The effect of non-ideal MHD on the results (compar- MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation t=1.30t ff A =0.2Ideal, B =-B z μ =10 A =0.1 A =0.05t=1.34t ff l o g c o l u m n d e n s i t y [ g / c m ] 400 AUt=1.30t ff A =0.2Ideal, B =-B z μ =5 A =0.1 A =0.05t=1.34t ff l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 3. Face-on gas column density snapshots from six ideal MHD mod-els with B = − B z at two different times (top and bottom row in eachsubfigure). The top subfigure shows results with µ = 10 , while the bot-tom subfigure shows the results with µ = 5 . Decreasing the amplitudeof the initial density perturbation A (left to right) or increasing the mag-netic field strength (i.e. decreasing µ ; top vs. bottom subfigure) decreasesthe first periastron separation as well as the disc mass and radius at firstperiastron. ing rows top to bottom in each subfigure) is small, with t peri differ-ing by less than a percent, though with R peri differing by up to per cent, or a maximum of au. The disc radii and masses differby up to and per cent, respectively. The largest differencesin disc mass between the ideal and non-ideal MHD calculationsoccurs in the calculations with µ = 10 and A = 0 . (left-handcolumn; top subfigure of Fig. 7), which also show the largest differ-ence in R peri . These discs are in the early stages of fragmentation,hence are irregularly shaped, which contributes to this differencein mass. A third protostar forms at t ≈ . t ff , disrupting the hostdiscs. In the ideal MHD model, a fourth protostar is formed at ap-proximately the same time to disrupt the second disc.At these early times, the main differences are caused bychanges in µ and A . This is expected since the density and mag-netic field strengths are only starting to reach the limits wherethe non-ideal effects become important. The left-hand subfigure inFig. 8 show the coefficients for Ohmic, Hall and ambipolar dif-fusion in the discs of the non-ideal models with B = − B z at t = 1 . t ff (we plot the average of the absolute value of the coeffi-cients). Ambipolar diffusion is the dominant effect in the disc, withthe coefficients of the Hall effect and Ohmic resistivity lower by afactor of ten in the disc.For comparison, these values are ∼ dex lower than in thediscs formed in the isolated collapse simulations shown by WPB16.Further, the plasma β is smaller in the isolated collapse models, so t=1.30t ff A =0.2Ideal, B =+B x μ =10 A =0.1 A =0.05t=1.34t ff l o g c o l u m n d e n s i t y [ g / c m ] 400 AUt=1.30t ff A =0.2Ideal, B =+B x μ =5 A =0.1 A =0.05t=1.34t ff l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 4. As in Fig. 3, but with B = + B x . With this magnetic fieldorientation, the initial magnetic field strength is the dominant parameter indetermining the evolution. In the strong magnetic field case (bottom subfig-ure), either no binaries or tight binaries form. we expect that non-ideal MHD will play a more minor role in thebinary case, as we have already shown. The µ = 10 models shown in the upper subfigure of Fig. 7 formlarge discs near first periastron. These discs subsequently fragment,with the non-ideal MHD discs fragmenting before their ideal MHDcounterparts. This fragmentation hinders the analysis of the latetime evolution of the weak field models, so in the rest of this sectionwe focus on the µ = 5 models.The left-hand column of Fig. 9 shows the gas column densityat t = 1 . t ff for the A = 0 . models, which is first apoastron forthese three models. The strong initial density perturbation yields anevolution that is very weakly dependent on the non-ideal MHD pro-cesses. Amongst these three models the first period and first apoas-tron differ by less than and per cent, respectively, comparedto and percent, respectively, for the models with A = 0 . .At this epoch, the difference between these three models are small,with the disc mass and magnetic field strengths at any given radiusdiffering by less than a factor of and , respectively.At second periastron, the discs in the non-ideal MHD mod-els fragment and form more protostars, totally disrupting the discs.This is mainly an artefact of our use of a barotropic equation ofstate to represent the thermodynamics.The middle column of Fig. 9 shows the gas column density ofthe µ = 5 , A = 0 . models at t = 1 . t ff , corresponding to thesecond apoastron in the ideal MHD model. At this epoch, the discsin the non-ideal MHD models are less massive and more extended MNRAS , 1–17 (2016) Wurster, Price & Bate S e p a r a ti on ( a u ) Ideal, B =+B x Non-ideal, B =+B x Non-ideal, B =-B x µ =10, A =0.1 0.0 0.1 0.2 0.3 0.4 0.5 M a ss ( M s un ) discstar+disc 0 20 40 60 80 100 120 R d i s c ( a u ) -1 p l a s m a β -4 -3 -2 -1 B ( G ) Time (t ff ) Figure 5. Time evolution of selected models starting from the formation ofthe protostars as in Fig. 2. Each model is initialised with an initial densityperturbation of A = 0 . , µ = 10 and B = ± B x . Top to bottom :the separation of the two protostars, the disc and star+disc masses, the discradius, and the mass-weighted plasma β and magnetic field strength in thedisc around one protostar. The ideal MHD model has a larger first peri-astron, resulting in a shorter first period than the non-ideal models. Thesubsequent periastron moves the ideal binary onto a long orbit, while thesubsequent orbits of the non-ideal models decrease. than their ideal MHD counterpart (comparing middle and bottompanels to the top panel). This is a result of the different orbital his-tories, which diverge shortly after first periastron at t peri ≈ . t ff .The blue and red lines in the right-hand column of Fig. 2 show theprotostar separation and disc properties for the non-ideal models,which can be directly compared to their ideal MHD counterpart(green line). In the ideal MHD model, there have been two perias-tron approaches by t = 1 . t ff , keeping the discs small and con-centrated around its host protostar; the non-ideal MHD discs havenot interacted with one another again, thus effectively evolved inisolation for the previous d t ≈ . t ff . Moreover, the values of thenon-ideal MHD coefficients rapidly decrease after periastron, mak-ing the later evolution more ideal. The coefficients increase brieflyat second periastron when the close interaction increases the den-sity of the disc. However, the interaction also leads to an increasein the value of plasma β , counteracting any added effect.The right-hand column of Fig. 9 shows the A = 0 . modelsat t = 1 . t ff . Additional clumps of gas form at t ≈ . t ff and r ≈ au from the centre of mass; at this time, the primary protostarsare r ≈ au from the centre of mass. The clumps spiral inwardsand interact with the primary binary starting at t ≈ . t ff , totallydisrupting the primary binary.Throughout this paper, we have compared models at the sameabsolute times. However, this may be an unfair comparison in somecases due to different orbital dynamics. For the above ideal and B =-B z A =0.2Ideal, t=1.34t ff μ =10 A =0.1 A =0.05B =+B x -3-2-1 l o g | B | 400 AUB =-B z A =0.2Ideal, t=1.34t ff μ =5 A =0.1 A =0.05B =+B x -3-2-1 l o g | B | 400 AU Figure 6. Magnetic field strength in the mid-plane ( z = 0 ) at t = 1 . t ff for the ideal MHD models with µ = 10 (top subfigure) and µ = 5 (bottom subfigure). At each initial magnetic field strength, the mid-planemagnetic field is always stronger in the models initialised with B = + B x ,despite the same initial value. non-ideal MHD models with A = 0 . (middle column of Fig 9), t apo ≈ . and . t ff , respectively; the non-ideal MHD modelsevolve very little between t = 1 . and . t ff , so the panels inFig. 9 are representative of both times. The disc mass and radius ofthe ideal MHD model, however, decreases from . M (cid:12) and au to . M (cid:12) and au, respectively, between first and secondapoastron. The intervening periastron passages strip mass from thedisc, concentrating the remaining disc mass closer to its host pro-tostar. This also results in a stronger magnetic field in the innerregions of the disc at t = 1 . t ff . We repeat the above study using a magnetic field initially perpen-dicular to the rotation axis. We consider both B = ± B x since weexpect a B z component to be generated during the evolution. Fig. 10 shows the gas column density of our suite of models with B = ± B x at t = 1 . t ff . The weaker field models ( µ = 10 ; topsubfigure) yield well-separated binaries for all A , with the non-ideal MHD models forming larger and more massive discs. Themagnetic field in the non-ideal MHD discs is approximately con-stant at B ≈ . G, while in the ideal MHD model, it decreasesby a factor of ∼ between the maximum strength and the outeredge of the disc. MNRAS000 Magnetic field strength in the mid-plane ( z = 0 ) at t = 1 . t ff for the ideal MHD models with µ = 10 (top subfigure) and µ = 5 (bottom subfigure). At each initial magnetic field strength, the mid-planemagnetic field is always stronger in the models initialised with B = + B x ,despite the same initial value. non-ideal MHD models with A = 0 . (middle column of Fig 9), t apo ≈ . and . t ff , respectively; the non-ideal MHD modelsevolve very little between t = 1 . and . t ff , so the panels inFig. 9 are representative of both times. The disc mass and radius ofthe ideal MHD model, however, decreases from . M (cid:12) and au to . M (cid:12) and au, respectively, between first and secondapoastron. The intervening periastron passages strip mass from thedisc, concentrating the remaining disc mass closer to its host pro-tostar. This also results in a stronger magnetic field in the innerregions of the disc at t = 1 . t ff . We repeat the above study using a magnetic field initially perpen-dicular to the rotation axis. We consider both B = ± B x since weexpect a B z component to be generated during the evolution. Fig. 10 shows the gas column density of our suite of models with B = ± B x at t = 1 . t ff . The weaker field models ( µ = 10 ; topsubfigure) yield well-separated binaries for all A , with the non-ideal MHD models forming larger and more massive discs. Themagnetic field in the non-ideal MHD discs is approximately con-stant at B ≈ . G, while in the ideal MHD model, it decreasesby a factor of ∼ between the maximum strength and the outeredge of the disc. MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation Ideal, B =-B z A =0.2μ =10t=1.34t ff A =0.1 A =0.05Non-ideal, B =-B z Non-ideal, B =+B z l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Ideal, B =-B z A =0.2μ =5t=1.34t ff A =0.1 A =0.05Non-ideal, B =-B z Non-ideal, B =+B z l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 7. Effect of non-ideal MHD on binary formation. Each panel showsthe face-on column density of a particular model at t = 1 . t ff . Top tobottom in each subfigure : ideal MHD, non-ideal MHD with − B z and with + B z , where the models in the top subfigure are initialised with µ = 10 ,while the bottom subfigure shows models with µ = 5 . The addition ofnon-ideal MHD has only a small effect compared to changing the initialmagnetic field strength (comparing top to bottom subfigures) or perturba-tion amplitude ( A ; comparing columns left to right). By contrast, the µ = 5 and A = 0 . models (left-hand col-umn, bottom subfigure) yield discs that are not significantly differ-ent from one another. The µ = 5 and A = 0 . calculations (mid-dle column, bottom subfigure) forms two protostars by t = 1 . t ff ,which all form a tight binary of semi-major axes ∼ and ∼ aufor the ideal and non-ideal MHD models, respectively. The largersemi-major axis in the non-ideal MHD models results in a largercentral cavity in the circumbinary disc.A single protostar forms in the ideal MHD model with A =0 . , while and protostars form by t = 1 . t ff in the non-idealMHD models with B = ± B x , respectively (right-hand column,bottom subfigure). This plethora of protostars immediately disruptsthe discs, and the remaining evolution is chaotic.As with the B = ± B z models, the dynamics are dominated by µ and A rather than the effect of non-ideal MHD. The right-hand subfigure in Fig. 8 shows the non-ideal MHD coefficients forthe non-ideal models with B = + B x at t = 1 . t ff . As withtheir B = − B z counterparts, ambipolar diffusion is the dominantterm, however, all the coefficients are lower despite the strongermagnetic field in the disc. The initial differences between the weak field models with A =0 . caused by non-ideal effects trigger pronounced differences asthe evolution continues. For example, the larger R peri for the idealmodel results in it reaching second periastron after d t = 0 . t ff ,while the non-ideal models require d t = 0 . t ff to reach second pe-riastron; see the top panel of Fig. 5 and left-hand column of Fig. 11.Subtle differences in mass and radius of the non-ideal MHD modelsnear second periastron causes their future evolution to diverge.Shortly after first periastron, all nine weak field models pro-duce an additional two protostars on orbits external to the primarybinary, and their early evolution is independent of the primary bi-naries. For A = 0 . and . , these external binaries do not in-teract with the primary prior to the end time of t = 2 . t ff , butthey interact near first apoastron in the non-ideal MHD models with A = 0 . .The µ = 5 models with A = 0 . and . retain a binaryuntil the end of the simulation at t = 1 . t ff , with an ellipticalbinary persisting in the former and a single, stable disc persistingin the latter. The middle column of Fig. 11 shows the A = 0 . models at t = 1 . t ff . At this time, the non-ideal MHD modelshave disc masses and radii that are 10 per cent larger and 4 per centsmaller, respectively, than their ideal counterpart. From top to bot-tom in that column, each model has an increasing periastron andapoastron distance, and by t = 1 . t ff , the models have passedthrough periastron , and times, respectively. Non-ideal MHDeffects contribute to these slight differences, but not enough to sig-nificantly change the overall evolution.The strong field models with A = 0 . form a single disc;see the right-hand column of Fig. 11 for gas column densities at t = 1 . t ff . The non-ideal MHD discs are ∼ per cent larger but ∼ per cent less massive as a result of the large central cavity.For all three models, the mass and radius decrease with time. Thegeneral trends amongst the three models are similar between theearly and late epochs. Fig. 12 compares the magnetic field strength in the mid-plane( z = 0 ) for the non-ideal MHD models at t = 1 . t ff ; this figure isdirectly comparable to Fig. 6. Only the B = − B z and + B x mod-els are shown, but the results are similar when the initial magneticfield direction is reversed.Similar to the ideal MHD models, at t = 1 . t ff , the net mag-netic field in the discs are higher in the B = + B x models thantheir − B z counterparts with the same µ and same A . When com-paring a non-ideal MHD model to its ideal MHD counterpart, themagnetic field in the mid-plane, and specifically the disc, is weaker(see also Figs. 2 and 5).As in Joos et al. (2013), we calculate the evolution of the mass-to-flux ratio inside a sphere of fixed radius, R , using µ ( R, t ) = M ( R ) πR (cid:104) B ( R ) (cid:105) (cid:18) M Φ B (cid:19) − crit , (14) MNRAS , 1–17 (2016) Wurster, Price & Bate B =-B z , t=1.34t ff µ =5, A =0.2 η ( c m s - ) η OR | η HE | η AD µ =10, A =0.210 µ =5, A =0.1 η ( c m s - ) µ =10, A =0.110 0 10 20 30 40 50 60 µ =5, A =0.05 η ( c m s - ) Radius (au) 0 10 20 30 40 50 60 µ =10, A =0.05 Radius (au) B =+B x , t=1.34t ff µ =5, A =0.2 η ( c m s - ) η OR | η HE | η AD µ =10, A =0.210 µ =5, A =0.1 η ( c m s - ) µ =10, A =0.110 0 10 20 30 40 50 60 µ =5, A =0.05 η ( c m s - ) Radius (au) 0 10 20 30 40 50 60 µ =10, A =0.05 Radius (au) Figure 8. Average values of Ohmic, ambipolar and Hall diffusion coefficients for the gas in the disc around one protostar at t = 1 . t ff averaged over allgas particles with ρ > ρ disc,min , for the models with B = − B z (left-hand subfigure) and B = + B x (right-hand subfigure). The Hall coefficient is theaverage of its absolute values. The lines switch to cyan at the defined edge of the disc. Ambipolar diffusion is the dominant effect in the disc, although allthree non-ideal coefficients typically differ by less than a factor of ten close to the protostar. These coefficients are smaller than in models that form an isolatedprotostar presented in WPB16. t=1.51t ff μ =5, A=0.2Ideal, B =-B z Non-ideal, B =-B z Non-ideal, B =+B z t=1.66t ff μ =5, A =0.1 t=1.51t ff μ =5, A=0.05 0123 l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 9. Gas column density for nine models with µ = 5 and B = ± B z at a later time. Left-hand column : the models with A = 0 . at t = 1 . t ff , which is first apoastron for all three models; the evolution isonly weakly dependent on the non-ideal effects. Middle column : the modelswith A = 0 . at t = 1 . t ff , which is second apoastron for the ideal MHDmodel; at this time, the ideal MHD model has discs that are more mas-sive and more concentrated near the protostar than in the non-ideal MHDmodels. Right-hand column : the models with A = 0 . at t = 1 . t ff ;additional protostars form at t ≈ . t ff , which interact with the primaryprotostars starting at t ≈ . t ff to disrupt the disc by this time. where M ( R ) is the enclosed mass including any protostars (i.e.sink particles), (cid:104) B ( R ) (cid:105) is the volume-averaged magnetic fieldwithin radius R , and ( M/ Φ B ) crit is the critical mass-to-flux ratiothat is independent of M , R and B . We have previous defined µ ≡ µ ( R = 0 . pc , t = 0) .The mass-to-flux ratio is spatially dependent, thus we plot fourvalues for selected models: Fig. 13 shows µ ( R, t ) for radii of R = 2680 au = 0 . pc (i.e. the initial size of the gas cloud) and au centred on the origin, and Fig. 14 shows µ ( R, t ) for radii of R = 120 and au centred on the first protostar that forms.For each non-ideal model and its ideal MHD counterpart, µ ( R = 2680 au , t ) is similar for all time, while µ ( R = 500 au , t ) begins to diverge at t (cid:38) . t ff . The small differences in µ ( R =500 au , t ) are a result of the weaker magnetic field and more mas-sive discs in the non-ideal models; the larger differences, includingthe sharp drops to µ ( R = 500 au , t ) (cid:28) are a result of the binaryseparation surpassing R = 1000 au; thus, these plots also provideinsight into the orbital properties of the binaries.These results are consistent with those found by Tassis &Mouschovias (2007b) and Joos et al. (2013) who studied the forma-tion of isolated protostars: the mass-to-flux ratio increases aroundthe protostar during its formation, removing memory of its initialvalue. In Tassis & Mouschovias (2007b), their mass-to-flux ratioincreases rapidly as the evolution proceeds and density increases.The increase of the mass-to-flux ratio centred on the first protostarin our models is more gradual after the formation of the protostarsince the sink particle removes the central magnetic field upon for-mation and particle accretion and effectively limits the maximumgas density at ρ ∼ − g cm − .As shown in Fig. 14, the evolution of µ ( R, t ) around the pro-tostar is quantitatively different for each model, reflecting the dif-ferent disc masses and magnetic fields contained within them. Bychanging any one parameter, the mass-to-flux ratio either increasesor decreases, indicating that no one parameter is dominant in deter-mining its evolution immediately around the protostar. The changecaused by including non-ideal MHD is typically smaller than thechange caused by altering another parameter, further suggestingthat non-ideal MHD plays a secondary role in binary formation. In previous studies of collapse to form an isolated first hydrostaticcore, the inclusion of non-ideal MHD was found to permit disc MNRAS000 MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation Ideal, B =+B x A =0.2μ =10t=1.34t ff A =0.1 A =0.05Non-ideal, B =+B x Non-ideal, B =-B x l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Ideal, B =+B x A =0.2μ =5t=1.34t ff A =0.1 A =0.05Non-ideal, B =+B x Non-ideal, B =-B x l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 10. Effect of non-ideal MHD on binary formation when the mag-netic field is initially perpendicular to the axis of rotation. Each panel showsthe face-on column density of a particular model t = 1 . t ff . Top to bot-tom in each subfigure : ideal MHD, non-ideal MHD with B = + B x andwith − B x , where the models in the top (bottom) subfigure are initialisedwith µ = 10 ( µ = 5 ). The non-ideal MHD models with µ = 5 and A = 0 . produce a tight binary with an orbit larger than their ideal MHDcounterpart, yielding a larger central cavity and less massive disc. The non-ideal MHD models with µ = 5 and A = 0 . produce multiple proto-stars, which immediately disrupt the system. For weak magnetic fields, thenon-ideal MHD models yield larger discs than their ideal counterparts. formation depending on the direction of the magnetic field withrespect to the axis of rotation (Tsukamoto et al. 2015b; WPB16).Moreover, the Hall effect was found to influence the formation pro-cess even when it was not the dominant non-ideal MHD effect.In this study, the faster initial rotation and the initial den-sity perturbations result in discs forming in every super-criticalmodel. Fig. 8 suggests that the non-ideal MHD effects may in-fluence mainly the inner regions of the discs. Thus, any changesthat switch the direction of the magnetic field in the disc will beamplified by the Hall effect, which may lead to simulations evolv-ing differently. However, these discs are primarily supported by gaspressure (i.e. β (cid:29) ), so the effect of the global change will depend t=1.90t ff Ideal, B =+B x μ =10, A=0.1Non-ideal, B =+B x =-B x t=1.51t ff μ =5, A=0.2400 AU t=1.51t ff μ =5, A=0.1 0123 l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 11. Gas column density for nine models with B = ± B x at alater times. Left-hand column : models with µ = 10 and A = 0 . at t = 1 . t ff ; the images in this column have frame sizes of (3000 au) sothat the four protostars in the non-ideal MHD models can be seen. Middlecolumn : models with µ = 5 and A = 0 . at t = 1 . t ff . Right-handcolumn : models with µ = 5 and A = 0 . at t = 1 . t ff . In all themodels presented, the ideal MHD models have more concentrated discs thantheir respective non-ideal MHD counterparts. B =-B z A =0.2Non-ideal, t=1.34t ff μ =10 A =0.1 A =0.05B =+B x -3-2-1 l o g | B | 400 AUB =-B z A =0.2Non-ideal, t=1.34t ff μ =5 A =0.1 A =0.05B =+B x -3-2-1 l o g | B | 400 AU Figure 12. As in Fig. 6, the magnetic field strength in the mid-plane ( z = 0 )at t = 1 . t ff , but for the non-ideal MHD models with µ = 10 (top sub-figure) and µ = 5 (bottom subfigure). As with the ideal MHD models, foreach initial magnetic field strength, the mid-plane magnetic field is alwaysstronger in the models initialised with B = + B x , despite the same initialvalue. The magnetic field is weaker in the disc mid-plane when comparedto their ideal MHD counterpart.MNRAS , 1–17 (2016) Wurster, Price & Bate B =-B z µ =5, A =0.2 µ ( R , t ) µ =10, A =0.210 µ =5, A =0.1 µ ( R , t ) µ =10, A =0.110 µ =5, A =0.05 µ ( R , t ) Time (t ff ) 0.5 1.0 1.5 2.0 2.5 µ =10, A =0.05 Time (t ff ) B =+B x µ =5, A =0.2 µ ( R , t ) µ =10, A =0.210 µ =5, A =0.1 µ ( R , t ) µ =10, A =0.110 µ =5, A =0.05 µ ( R , t ) Time (t ff ) 0.5 1.0 1.5 2.0 2.5 µ =10, A =0.05 Time (t ff ) Figure 13. Evolution of the mass-to-flux ratio, µ ( R, t ) , for R = 2680 au = 0 . pc (i.e. the initial size of the gas cloud) and au for the ideal (dashedlines) and non-ideal (solid lines) MHD models. The lines end when one of the models in the panel has reached its end time. The cyan line represent the initialmass-to-flux ratio, µ . The ratio is µ ( R, t ) ∝ M ( R ) / (cid:104) B ( R ) (cid:105) , where M ( R ) is the mass enclosed within a sphere of radius R centred on the origin, and (cid:104) B ( R ) (cid:105) is the volume-averaged magnetic field within the sphere. In each panel, µ ( R = 2680 au , t ) is similar for both models, whereas µ ( R = 500 au , t ) may differ after t (cid:38) . t ff due to difference in disc masses and magnetic fields (small deviations) or due to part or all of the discs leaving the sphere (bigdeviations and sudden drops). B =-B z µ =5, A =0.2 µ ( R , t ) µ =10, A =0.210 µ =5, A =0.1 µ ( R , t ) µ =10, A =0.110 µ =5, A =0.05 µ ( R , t ) Time (t ff ) 1.4 1.6 1.8 2.0 2.2 2.4 µ =10, A =0.05 Time (t ff ) B =+B x µ =5, A =0.2 µ ( R , t ) µ =10, A =0.210 µ =5, A =0.1 µ ( R , t ) µ =10, A =0.110 µ =5, A =0.05 µ ( R , t ) Time (t ff ) 1.4 1.6 1.8 2.0 2.2 2.4 µ =10, A =0.05 Time (t ff ) Figure 14. Evolution of the mass-to-flux ratio, µ ( R, t ) for the ideal (dashed lines) and non-ideal (solid lines) MHD models as in Fig. 13, but for R = 120 and au centred on the first protostar that forms. The differences in µ ( R, t ) between each ideal/non-ideal pair are similar to the differences caused by changingother parameters, suggesting that non-ideal MHD plays a secondary role in binary formation. on the amount of modification by the Hall effect and when it occurs.For example, the initial direction of the magnetic field plays a min-imal role in the evolution of the models with µ = 5 , A ≥ . and B = ± B z , while it triggers a divergence in the evolutionarypaths of µ = 5 , A = 0 . and B = ± B x .In the earlier studies, the Hall effect was found to spin up thedisc when the magnetic field was initially anti-aligned with the axisof rotation. To conserve angular momentum, a counter-rotating en-velope forms (e.g. Tsukamoto et al. 2015b; WPB16). By contrast,the Hall effect does not have a pronounced effect on the rotation ofthe discs in our binary models, and as a result we see no evidencefor counter-rotating envelopes. In the above models, the initial mass-to-flux ratio of µ > meansthat the magnetic field is not strong enough to prevent gravitationalcollapse, thus protostar formation is a foregone conclusion.A gas cloud with a sub-critical mass-to-flux ratio is magnet-ically supported and should not collapse when using ideal MHD.Indeed, our sub-critical models with µ = 0 . do not collapseduring their runtime to t ≈ t ff , and their maximum density neversurpasses ρ ≈ × − g cm − (recall that the initial density is ρ = 7 . × − g cm − and that sink particles are inserted at MNRAS000 Evolution of the mass-to-flux ratio, µ ( R, t ) for the ideal (dashed lines) and non-ideal (solid lines) MHD models as in Fig. 13, but for R = 120 and au centred on the first protostar that forms. The differences in µ ( R, t ) between each ideal/non-ideal pair are similar to the differences caused by changingother parameters, suggesting that non-ideal MHD plays a secondary role in binary formation. on the amount of modification by the Hall effect and when it occurs.For example, the initial direction of the magnetic field plays a min-imal role in the evolution of the models with µ = 5 , A ≥ . and B = ± B z , while it triggers a divergence in the evolutionarypaths of µ = 5 , A = 0 . and B = ± B x .In the earlier studies, the Hall effect was found to spin up thedisc when the magnetic field was initially anti-aligned with the axisof rotation. To conserve angular momentum, a counter-rotating en-velope forms (e.g. Tsukamoto et al. 2015b; WPB16). By contrast,the Hall effect does not have a pronounced effect on the rotation ofthe discs in our binary models, and as a result we see no evidencefor counter-rotating envelopes. In the above models, the initial mass-to-flux ratio of µ > meansthat the magnetic field is not strong enough to prevent gravitationalcollapse, thus protostar formation is a foregone conclusion.A gas cloud with a sub-critical mass-to-flux ratio is magnet-ically supported and should not collapse when using ideal MHD.Indeed, our sub-critical models with µ = 0 . do not collapseduring their runtime to t ≈ t ff , and their maximum density neversurpasses ρ ≈ × − g cm − (recall that the initial density is ρ = 7 . × − g cm − and that sink particles are inserted at MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation -1 B =-B z µ =0.75, A =0.1 µ ( R , t ) -1 0 1 2 3 4 5 6 7 8B =+B x µ =0.75, A =0.1 µ ( R , t ) Time (t ff ) Figure 15. Evolution of the mass-to-flux ratio, µ ( R, t ) , for the ideal(dashed lines) and non-ideal (solid lines) MHD models as in Fig. 13, butfor the sub-critical models with µ = 0 . . The cyan lines are at µ = µ and µ = 1 (i.e. the critical value). The vertical lines represent when theprotostars formed in the non-ideal MHD models. The ideal MHD modelsdo not collapse during their runtime of t ≈ t ff while the non-ideal MHDmodels form protostars at t = 5 . and . t ff , for B = − B z and B x ,respectively. ρ ≈ − g cm − ). As shown by the dashed lines in Fig. 15, after t (cid:38) t ff , the mass-to-flux ratio at R = 2680 , and au areapproximately constant, with µ ( R = 200 au , t ) < . The value of (cid:104) B ( R ) (cid:105) is similar for spheres of both R = 500 and au, butthe former has more enclosed mass, hence the higher mass-to-fluxratio.For the non-ideal MHD models with B = − B z and B x , µ ( R = 200 au , t ) > at t ≈ . and . t ff , respectively, as shownby the solid green lines in Fig. 15. After this time, the central re-gions are no longer magnetically supported, and the clouds collapseto form protostars at t = 5 . and . t ff , respectively. Fig. 16shows the gas column density for the sub-critical models near thetime of protostar formation for the non-ideal MHD models. Onlyone protostar is formed in each model, and no discs form aroundthem.Thus, in our models, the non-ideal MHD effects can diffuseenough magnetic field to allow the central regions of initially sub-critical clouds to collapse to form protostars. The models in our primary suite all use an initial rotation of Ω = 1 . × − rad s − , and binaries form in all but fourof the 36 models. However, (e.g.) Hennebelle & Teyssier (2008)and Machida et al. (2008b) found that initial rotation played an im-portant role in determining the evolution of the system. Althougha full parameter study of the initial rotation is out of the scope ofthis study, we briefly discuss the early evolution of models with theslower initial rotations of Ω = 7 . × − , . × − and . × − rad s − , using µ = 5 , A = 0 . and B = − B z . Idealμ =0.75, B =-B z A =0.1 t=5.66t ff t=5.71t ff t=5.72t ff Non-ideal -202 l o g c o l u m n d e n s i t y [ g / c m ] 150 AU Idealμ =0.75, B =+B x A =0.1 t=5.79t ff t=5.84t ff t=5.85t ff Non-ideal -202 l o g c o l u m n d e n s i t y [ g / c m ] 150 AU Figure 16. Face-on gas column density of the initially sub-critical models( µ = 0 . ) using ideal (top row in each subfigure) and non-ideal (bottomrow in each subfigure) MHD. For both initial magnetic field orientations,the ideal MHD models do not collapse, with their density staying below ρ ≈ × − g cm − . The non-ideal MHD models collapse to formsingle protostars shortly after t = 5 . and . t ff for B = − B z (topsubfigure) and B x (bottom subfigure), respectively. Fig. 17 shows the gas column density of these models at earlytimes.As the initial rotation speed decreases, the initial separationof the binaries decreases, with only a single protostar forming atthe slowest two rotation speeds. The non-ideal MHD effects be-come more important as the initial rotation decreases, with largerand more massive discs forming at slower rotation speeds. This isconsistent with previous studies finding larger discs in non-idealMHD models of isolated protostars than ideal MHD models. How-ever, the non-ideal MHD effects do not change the global morphol-ogy and whether binary or single systems form, which is consistentwith our previous results suggesting that non-ideal MHD has a sec-ondary effect on binary formation and evolution. As this study was in progress, Version 1.2.1 of N ICIL (Wurster2016) was released. This version differs from v1.1 used here bymodelling three grain populations, n − g , n g and n + g , with charges Z = − , , +1 , respectively, rather than a single grain population, n g , with charge ¯ Z < . In v1.2.1, grain number density is con-served, with n g = n − g + n g + n + g , where n g is calculated as inv1.1.To test the effect of the grain model, we run two additionalmodels using v1.2.1: µ = 10 , A = 0 . , B = + B x and µ = 5 , A = 0 . , B = − B z . Fig. 18 shows the gas column density at MNRAS , 1–17 (2016) Wurster, Price & Bate t=1.3t ff IdealB =-B z Ω =1.01x10 -12 s -1 μ =5A =0.1 t=1.31t ff t=1.34t ff t=1.14t ff Ω =7.08x10 -13 s -1 t=1.15t ff t=1.18t ff t=1.04t ff Ω =3.54x10 -13 s -1 t=1.05t ff t=1.08t ff t=1.02t ff Ω =1.77x10 -13 s -1 t=1.03t ff l o g c o l u m n d e n s i t y [ g / c m ] t=1.06t ff 400 AU t=1.3t ff Non-idealB =-B z Ω =1.01x10 -12 s -1 μ =5A =0.1 t=1.31t ff t=1.34t ff t=1.14t ff Ω =7.08x10 -13 s -1 t=1.15t ff t=1.18t ff t=1.04t ff Ω =3.54x10 -13 s -1 t=1.05t ff t=1.08t ff t=1.02t ff Ω =1.77x10 -13 s -1 t=1.03t ff l o g c o l u m n d e n s i t y [ g / c m ] t=1.06t ff 400 AU Figure 17. Face-on gas column density of models using the fiducial ini-tial rotation of Ω = 1 . × − rad s − and the slower rotations of Ω = 7 . × − , . × − and . × − rad s − . All modelsuse µ = 5 , A = 0 . , B = − B z , and ideal (top subfigure) or non-ideal(bottom subfigure) MHD. The panels are chosen such that the protostarsform between the first two columns, and the third column is d t ≈ . t ff af-ter the protostar’s formation. Decreasing the initial rotation speed decreasesthe initial binary separation and, if slow enough, prevents the formation ofbinaries. Non-ideal MHD has a greater influence on the environment of theinitially slower rotating models, forming larger and more massive discs, buthas little effect on the large-scale morphology or the number of protostarsthat form. t=1.34t ff B =+B x , μ =10A =0.1Ideal t=1.40t ff t=1.46t ff Non-ideal, NICIL v1.1Non-ideal, NICIL v1.2.1 0123 l o g c o l u m n d e n s i t y [ g / c m ] 400 AU t=1.34t ff B =-B z , μ =5A =0.01Ideal t=1.40t ff t=1.46t ff Non-ideal, NICIL v1.1Non-ideal, NICIL v1.2.1 0123 l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 18. Face-on gas column density from selected models using idealMHD (top row in each subfigure), NICIL Version 1.1 (middle row in eachsubfigure) and NICIL Version 1.2.1 (bottom row in each subfigure). Ver-sion 1.1 uses a single grain population, n g , with charge ¯ Z < and Ver-sion 1.2.1 models three grain populations, n − g , n g and n + g , with charges Z = − , , +1 , respectively. The models with three grain populationsyield binaries with smaller first periastron separations, and, for the µ = 5 models, first apoastron and first periods that are ∼ . times smaller. selected times, and Fig. 19 shows the radial profile of the grainpopulations and non-ideal MHD coefficients in the disc around oneprotostar at t = 1 . t ff .Modelling three grain populations yields binaries with smallerfirst periastron separations than the single grain model. The firstapoastron separation and first period are similar for the µ = 10 models, but for the µ = 5 models, they are ∼ . times largerwhen using v1.1.Once the discs form, the neutral number density is onlyweakly dependent on the grain model, and is n n ∼ O (10 ) cm − for the duration of the simulation. Averaged over the entire disc, thetotal grain number density differs by (cid:46) between the two grainmodels, with the largest differences occurring at larger radii andat later times during the evolution. At t = 1 . t ff and during the MNRAS000 Face-on gas column density from selected models using idealMHD (top row in each subfigure), NICIL Version 1.1 (middle row in eachsubfigure) and NICIL Version 1.2.1 (bottom row in each subfigure). Ver-sion 1.1 uses a single grain population, n g , with charge ¯ Z < and Ver-sion 1.2.1 models three grain populations, n − g , n g and n + g , with charges Z = − , , +1 , respectively. The models with three grain populationsyield binaries with smaller first periastron separations, and, for the µ = 5 models, first apoastron and first periods that are ∼ . times smaller. selected times, and Fig. 19 shows the radial profile of the grainpopulations and non-ideal MHD coefficients in the disc around oneprotostar at t = 1 . t ff .Modelling three grain populations yields binaries with smallerfirst periastron separations than the single grain model. The firstapoastron separation and first period are similar for the µ = 10 models, but for the µ = 5 models, they are ∼ . times largerwhen using v1.1.Once the discs form, the neutral number density is onlyweakly dependent on the grain model, and is n n ∼ O (10 ) cm − for the duration of the simulation. Averaged over the entire disc, thetotal grain number density differs by (cid:46) between the two grainmodels, with the largest differences occurring at larger radii andat later times during the evolution. At t = 1 . t ff and during the MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation [!h] -2 -1 B =+B x , t=1.34t ff µ =10, A =0.1 n g ( c m - ) n g (1 species)n g |Z g |n g (3 species)n g (Z=-1)n g (Z= 0)n g (Z=+1)10 0 10 20 30 40 η ( c m s - ) Radius (au) η OR | η HE | η AD -2 -1 B =-B z , t=1.34t ff µ =5, A =0.1 n g ( c m - ) n g (1 species)n g |Z g |n g (3 species)n g (Z=-1)n g (Z= 0)n g (Z=+1)10 0 10 20 30 40 η ( c m s - ) Radius (au) η OR | η HE | η AD Figure 19. Top panel in each subfigure : The grain number density in the discaround one protostar at t = 1 . t ff using v1.1 of NICIL (dashed lines) andv1.2.1 (solid lines). For v1.1, there is a single grain population with an aver-age negative charge, thus the dashed blue line is n g | Z g | , which effectivelyrepresents n − g . For v1.2.1, all three grain populations are self-consistentlycalculated, and n g = n − g + n g + n + g . For reference, the neutral grain num-ber density is n n ∼ O (10 ) cm − . While the total grain number density, n g , is only weakly dependent on grain model, the effective, n g | Z g | , andreal, n − g , number densities of negatively charged grains can differ by anorder of magnitude. Bottom panel in each subfigure : As in Fig. 8, averagevalues of Ohmic, ambipolar and Hall diffusion coefficients for the gas. Thegrain model only weakly affects the non-ideal MHD coefficients for weakmagnetic fields, but decreases the values of ambipolar diffusion and Ohmicresistivity in the high density regions of the disc in the strong field models. Non-ideal, μ =10, A =0.1 B =-B z N=3x10 t=1.34t ff B =+B x N=10 l o g c o l u m n d e n s i t y [ g / c m ] 400 AU Figure 20. Gas column density for four non-ideal MHD models with µ =10 and A = 0 . at t = 1 . t ff . The left- (right-) hand column shows themodels with B = − B z ( B = + B x ), and the models in the top andbottom rows are initialised with × and particles, respectively,in the sphere. The higher resolution models produce somewhat larger andmore massive discs that are more susceptible to fragmentation. evolution of the disc, n g > n − g > n + g . To approximate n − g inthe single grain model, we use n − g ≈ n g | Z g | , which varies onlyslightly with both radius and time. This value is ∼ times smallerfor our µ = 10 model, and up to ∼ times smaller for the µ = 5 model.The different grain models affect the calculation of the non-ideal MHD coefficients, η , since we are improving the calculationof n − g and adding a charged species, n + g . When comparing η for thedifferent grain models (bottom panel of each subfigure in Fig. 19) inthe weak field model ( µ = 10 ), ambipolar diffusion is the largestterm for both grain models, followed by the Hall effect. At thistime in the µ = 5 model, all values are similar in the inner disc( R ∼ au, n n ∼ × cm − ), however, the order of theterms differs depending on the grain model; ambipolar diffusion isreduced in strength to be similar to Ohmic resistivity at this radius.This is consistent with Tassis & Mouschovias (2007b) and Kunz& Mouschovias (2010), who find that Ohmic resistivity becomesmore important than ambipolar diffusion at n n ∼ O (10 ) cm − .When considering the average values of η over the entire disc,we find that η AD > | η HE | > η OR , and that n n (cid:46) cm − . Thus,we expect ambipolar diffusion to dominate in our models. The calculations presented above used × particles in the ini-tial sphere. This number satisfies the Jeans criteria (c.f. Section 3),while allowing us to perform a large suite of simulations, even withthe small timesteps required to properly evolve the non-ideal MHDterms. To test the effect of resolution, we ran selected models ini-tialised with particles in the sphere.Fig. 20 shows the non-ideal MHD models at t = 1 . t ff with µ = 10 and A = 0 . using both × and particles in theinitial sphere. The particle models yield discs that are some-what larger and more massive than their lower resolution coun-terparts. This results in greater interaction at first periastron andshorter periods. Although the quantitative values change with res-olution, the trends previously discussed are independent of resolu-tion. The higher resolution discs are more massive and thus moresusceptible to fragmentation. We note this fragmentation would MNRAS , 1–17 (2016) Wurster, Price & Bate likely not be a problem if we included a proper treatment of ra-diation rather than the barotropic equation of state used here. We have presented a suite of simulations studying the effect of non-ideal MHD on the formation and early evolution of binary stars.Our models were initialised as a M (cid:12) rotating, uniform densitysphere, which was given an m = 2 density perturbation with am-plitudes of A = 0 . , . and . . We threaded the sphere with aninitially uniform magnetic field. We tested our suite of simulationsusing both ideal MHD and non-ideal MHD at initial mass-to-fluxratios of 10, 5 and 0.75 times the critical value. The ideal MHDmodels were run using B = − B z and + B x , while the non-idealMHD models were run using B = ± B z and ± B x since the Halleffect depends on the direction of the magnetic field with respectto the axis of rotation. In the models that formed binaries, we fol-lowed the gravitational collapse until at least first apoastron. Allof the simulations were performed using the SPMHD code P HAN - TOM . Our key results are as follows:(i) Sub-critical cores : Using ideal MHD, the sub-critical coresdid not collapse during their runtime of t ≈ t ff , and their max-imum density never surpassed ρ ≈ ρ . When using non-idealMHD, the cores collapsed to form single protostars at t (cid:46) . t ff .(ii) Ideal MHD : B = − B z : Decreasing the amplitude of theinitial density perturbation yields earlier times of first periastron,and smaller separations. Decreasing the magnetic field strength (i.e.increasing µ ) increases first periastron separation and discs sizes.(iii) Ideal MHD : B = + B x : Strong magnetic fields suppressthe formation of binaries, as found by previous authors, with a bi-nary only forming for A = 0 . ; a tight binary with a commondisc forms for A = 0 . and a single protostar and disc forms for A = 0 . . Binaries form in all the weak field models, with largerfirst periastron separations and smaller discs masses and radii thanin their B = − B z counterparts. The magnetic fields in the discare ∼ times stronger than in their B = − B z counterparts, de-spite the same initial strength in the initial cloud.(iv) Non-ideal MHD : B = ± B z : The time of first periastron isnot affected by the inclusion of non-ideal MHD, however, at latertimes, the non-ideal MHD models tend to have longer periods andlarger apoastron separations than the ideal MHD models, as well aslarger and more massive discs. When discs become massive enoughto fragment, the fragmentation occurs in the non-ideal MHD mod-els more easily. The evolution of the − B z and B z models tendsto diverge between first apoastron and second periastron; the dif-ferences are initially small, but are enhanced by the dynamics andsubsequent interactions, rather than influences of non-ideal MHD.(v) Non-ideal MHD : B = ± B x : The non-ideal MHD mod-els that form binaries yield smaller first periastron separations andlarger disc radii compared to their ideal MHD counterparts. Diver-gence in periods and periastron and apoastron separations betweenideal and non-ideal MHD models occurs shortly after first perias-tron, while the divergence between the two non-ideal MHD modelswith the same µ and A occurs later, once local changes havemodified the vertical component of the magnetic field such that theHall effect produces a different evolution in each model.(vi) The Hall effect : Unlike models that form an isolated proto-star, the Hall effect does not have a global impact on the evolutionof the binaries. Rather, local changes to the magnetic field will beenhanced by the Hall effect causing small modification to the evo- lution of the model. These small modifications are further enhancedat periastron, causing the evolutionary paths to slowly diverge.(vii) Rotation speeds : Decreasing the initial rotation speed hin-dered binary formation. At the slowest two speeds tested, onlya single protostar formed. The inclusion of non-ideal MHD didnot affect the global morphology or the number of protostars thatformed.(viii) Grain model : The same qualitative conclusions arereached if the non-ideal MHD algorithm used a single populationwith average charge or three separate grain populations. The firstperiastrons were smaller in the models that used the three grainpopulations.In the formation of binary systems, the initial parameters —amplitude of the initial density perturbation, magnetic field strengthand orientation, and rotation — determine their evolution. The maineffect of non-ideal MHD is to enable the formation of larger andmore massive discs form around the protostars, and produce bina-ries that have larger separations and longer periods. ACKNOWLEDGEMENTS JW and MRB acknowledge support from the European ResearchCouncil under the European Community’s Seventh FrameworkProgramme (FP7/2007- 2013 grant agreement no. 339248). JWalso acknowledges support from the Australian Research Council(ARC) Discovery Projects Grant DP130102078. DJP is funded byARC Future Fellowship FT130100034. This work was supportedby resources on the gSTAR national facility at Swinburne Univer-sity of Technology and by Zen. gSTAR is funded by Swinburne andthe Australian Government’s Education Investment Fund. Severalcalculations for this paper were performed on the University of Ex-eter Supercomputer, a DiRAC Facility jointly funded by STFC, theLarge Facilities Capital Fund of BIS, and the University of Exeter.For the column density figures, we used SPLASH (Price 2007). REFERENCES Alexiades V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Meth.Eng., 12, 31Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481Bate M. R., 2011, MNRAS, 417, 2036Bate M. R., Burkert A., 1997, MNRAS, 288, 1060Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Boss A. P., 1986, ApJS, 62, 519Boss A. P., Bodenheimer P., 1979, ApJ, 234, 289Boss A. P., Keiser S. A., 2013, ApJ, 764, 136Braiding C. R., Wardle M., 2012, MNRAS, 427, 3188Burkert A., Bodenheimer P., 1993, MNRAS, 264, 798Ciolek G. E., Mouschovias T. C., 1994, ApJ, 425, 142Dapp W. B., Basu S., 2010, A&A, 521, L56Dapp W. B., Basu S., Kunz M. W., 2012, A&A, 541, A35Duffin D. F., Pudritz R. E., 2009, ApJ, 706, L46Duquennoy A., Mayor M., 1991, A&A, 248, 485Fiedler R. A., Mouschovias T. C., 1993, ApJ, 415, 680Gafton E., Rosswog S., 2011, MNRAS, 418, 770Hennebelle P., Teyssier R., 2008, A&A, 477, 25Joos M., Hennebelle P., Ciardi A., Fromang S., 2013, A&A, 554, A17Keith S. L., Wardle M., 2014, MNRAS, 440, 89Kunz M. W., Mouschovias T. C., 2010, MNRAS, 408, 322Larson R. B., 1969, MNRAS, 145, 271Lewis B. T., Bate M. R., Price D. J., 2015, MNRAS, 451, 288Li Z.-Y., Shu F. H., 1996, ApJ, 472, 211 MNRAS000 Alexiades V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Meth.Eng., 12, 31Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481Bate M. R., 2011, MNRAS, 417, 2036Bate M. R., Burkert A., 1997, MNRAS, 288, 1060Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Boss A. P., 1986, ApJS, 62, 519Boss A. P., Bodenheimer P., 1979, ApJ, 234, 289Boss A. P., Keiser S. A., 2013, ApJ, 764, 136Braiding C. R., Wardle M., 2012, MNRAS, 427, 3188Burkert A., Bodenheimer P., 1993, MNRAS, 264, 798Ciolek G. E., Mouschovias T. C., 1994, ApJ, 425, 142Dapp W. B., Basu S., 2010, A&A, 521, L56Dapp W. B., Basu S., Kunz M. W., 2012, A&A, 541, A35Duffin D. F., Pudritz R. E., 2009, ApJ, 706, L46Duquennoy A., Mayor M., 1991, A&A, 248, 485Fiedler R. A., Mouschovias T. C., 1993, ApJ, 415, 680Gafton E., Rosswog S., 2011, MNRAS, 418, 770Hennebelle P., Teyssier R., 2008, A&A, 477, 25Joos M., Hennebelle P., Ciardi A., Fromang S., 2013, A&A, 554, A17Keith S. L., Wardle M., 2014, MNRAS, 440, 89Kunz M. W., Mouschovias T. C., 2010, MNRAS, 408, 322Larson R. B., 1969, MNRAS, 145, 271Lewis B. T., Bate M. R., Price D. J., 2015, MNRAS, 451, 288Li Z.-Y., Shu F. H., 1996, ApJ, 472, 211 MNRAS000 , 1–17 (2016) on-ideal MHD and binary formation Li Z.-Y., Krasnopolsky R., Shang H., 2011, ApJ, 738, 180Lodato G., Price D. J., 2010, MNRAS, 405, 1212Machida M. N., Inutsuka S.-i., Matsumoto T., 2008a, ApJ, 676, 1088Machida M. N., Omukai K., Matsumoto T., Inutsuka S.-i., 2008b, ApJ, 677,813Machida M. N., Inutsuka S.-I., Matsumoto T., 2011, PASJ, 63, 555Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350Mellon R. R., Li Z.-Y., 2009, ApJ, 698, 922Mestel L., Spitzer Jr. L., 1956, MNRAS, 116, 503Morris J. P., 1996, PASA, 13, 97Mouschovias T. C., 1996, in Tsinganos K. C., Ferrari A., eds, NATO Ad-vanced Science Institutes (ASI) Series C Vol. 481, NATO AdvancedScience Institutes (ASI) Series C. pp 505–538Mouschovias T. C., Ciolek G. E., 1999, in Lada C. J., Kylafis N. D., eds,NATO Advanced Science Institutes (ASI) Series C Vol. 540, NATOAdvanced Science Institutes (ASI) Series C. p. 305Mouschovias T. C., Spitzer Jr. L., 1976, ApJ, 210, 326Nakano T., Umebayashi T., 1986a, MNRAS, 218, 663Nakano T., Umebayashi T., 1986b, MNRAS, 221, 319Nakano T., Nishi R., Umebayashi T., 2002, ApJ, 573, 199Pandey B. P., Wardle M., 2008, MNRAS, 385, 2269Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., FongW., 1994, ApJ, 421, 615Price D. J., 2007, PASA, 24, 159Price D. J., 2012, Journal of Computational Physics, 231, 759Price D. J., Bate M. R., 2007, MNRAS, 377, 77Price D. J., Federrath C., 2010, MNRAS, 406, 1659Price D. J., Monaghan J. J., 2004, MNRAS, 348, 123Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347Price D. J., Tricco T. S., Bate M. R., 2012, MNRAS, 423, L45Raghavan D., et al., 2010, ApJS, 190, 1Shu F. H., Galli D., Lizano S., Cai M., 2006, ApJ, 647, 382Tassis K., Mouschovias T. C., 2007a, ApJ, 660, 388Tassis K., Mouschovias T. C., 2007b, ApJ, 660, 402Tomida K., Machida M. N., Saigo K., Tomisaka K., Matsumoto T., 2010,ApJ, 725, L239Tomida K., Tomisaka K., Matsumoto T., Hori Y., Okuzumi S., MachidaM. N., Saigo K., 2013, ApJ, 763, 6Tomida K., Okuzumi S., Machida M. N., 2015, ApJ, 801, 117Tricco T. S., Price D. J., 2012, Journal of Computational Physics, 231, 7214Tricco T. S., Price D. J., Bate M. R., 2016, Journal of ComputationalPhysics, 322, 326Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015a,MNRAS, 452, 278Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015b,ApJ, 810, L26Umebayashi T., Nakano T., 1990, MNRAS, 243, 103Wardle M., 2007, Ap&SS, 311, 35Wardle M., Ng C., 1999, MNRAS, 303, 239Wurster J., 2016, PASA, 33, e041Wurster J., Price D., Ayliffe B., 2014, MNRAS, 444, 1104Wurster J., Price D. J., Bate M. R., 2016, MNRAS, 457, 1037 APPENDIX A: MODELS WITH INITIAL CONDITIONSAND SELECTED RESULTS Table A1 summaries the initial parameters of all our models, alongwith the time of first periastron t peri , the initial period T , and thefirst periastron R peri and first apoastron R apo separations. MNRAS , 1–17 (2016) Wurster, Price & Bate µ A B MHD Alternate t peri T R peri R apo parameter ( t ff ) ( t ff ) (au) (au)5 0.2 − B z ideal 1.370 0.28 110 5305 0.1 − B z ideal 1.333 0.25 68 4405 0.05 − B z ideal 1.318 0.097 49 21010 0.2 − B z ideal 1.369 n/a:P 140 n/a:P10 0.1 − B z ideal 1.340 1.7 100 280010 0.05 − B z ideal 1.324 0.67 54 11005 0.2 B x ideal 1.287 0.033 39 1005 0.1 B x ideal 0 0 0 05 0.05 B x ideal 0 0 0 010 0.2 B x ideal 1.320 0.89 210 120010 0.1 B x ideal 1.330 0.39 190 65010 0.05 B x ideal 1.345 0.39 190 5705 0.2 − B z non-ideal 1.372 0.29 95 5005 0.1 − B z non-ideal 1.337 0.52 61 8205 0.05 − B z non-ideal 1.320 0.12 43 2005 0.2 B z non-ideal 1.373 0.24 98 4605 0.1 B z non-ideal 1.337 0.53 63 8005 0.05 B z non-ideal 1.320 0.12 43 20010 0.2 − B z non-ideal 1.374 n/a:P 120 n/a:P10 0.1 − B z non-ideal 1.341 n/a:A 110 n/a:A10 0.05 − B z non-ideal 1.337 n/a:A 62 n/a:A10 0.2 B z non-ideal 1.378 n/a:P 120 n/a:P10 0.2 B z non-ideal 1.340 n/a:A 110 n/a:A10 0.05 B z non-ideal 1.325 n/a:A 58 n/a:A5 0.2 B x non-ideal 1.288 0.036 14 1305 0.1 B x non-ideal 0 0 0 05 0.05 B x non-ideal n/a:P n/a:P n/a:P n/a:P5 0.2 − B x non-ideal 1.288 0.035 14 1205 0.1 − B x non-ideal 0 0 0 05 0.05 − B x non-ideal n/a:P n/a:P n/a:P n/a:P10 0.2 B x non-ideal 1.347 n/a:A 120 n/a:A10 0.1 B x non-ideal 1.332 0.79 150 110010 0.05 B x non-ideal 1.324 0.25 190 47010 0.2 − B x non-ideal 1.347 n/a:A 120 n/a:A10 0.1 − B x non-ideal 1.332 0.79 150 110010 0.05 − B x non-ideal 1.324 0.29 190 4800.75 0.1 − B z ideal sub-critical µ n/a:NC n/a:NC n/a:NC n/a:NC0.75 0.1 B x ideal sub-critical µ n/a:NC n/a:NC n/a:NC n/a:NC0.75 0.1 − B z non-ideal sub-critical µ B x non-ideal sub-critical µ − B z ideal Ω = 7 . × − s − − B z ideal Ω = 3 . × − s − − B z ideal Ω = 1 . × − s − − B z non-ideal Ω = 7 . × − s − − B z non-ideal Ω = 3 . × − s − − B z non-ideal Ω = 1 . × − s − − B z non-ideal NICIL v1.2.1 1.336 0.17 56 35010 0.1 B x non-ideal NICIL v1.2.1 1.332 0.79 130 1100 Table A1. The initial parameters and early results of our suite of models. The first four columns are the initial conditions: the initial mass-to-flux ratio µ inunits of the critical mass-to-flux ratio, the amplitude of the initial density perturbation A , the initial orientation of the magnetic field B , and whether themodel uses ideal or non-ideal MHD. The fifth column lists deviations from the parameters used in the primary suite (see Section 4.8). The remaining columnsare the time of first periastron t peri , the initial period T calculated using first periastron and first apoastron, and the separations at first periastron R peri and firstapoastron R apo . Entires with zeros indicate that only one disc formed, thus T , R peri and R apo do not exists. Entires with n/a:P indicate that one or both discsfragmented near first periastron and became unstable, entires with n/a:A indicate that one or both discs interacted with younger protostars near first apoastronwhich modified the primary binary’s orbit, and entries with n/a:NC indicate that the cloud did not collapse to form protostars; in these three cases, separationsand periods have no useful meaning. MNRAS000