The role of magnetic fields in the formation of multiple massive stars

(Abridged) Context. Most massive stars are located in multiple stellar systems. Magnetic fields are believed to be essential in the accretion and ejection processes around single massive protostars. Aims. Our aim is to unveil the influence of magnetic fields in the formation of multiple massive stars, in particular on the fragmentation modes and properties of the multiple protostellar system. Methods. Using RAMSES, we follow the collapse of a massive pre-stellar core with (non-ideal) radiation-(magneto-)hydrodynamics. We choose a setup which promotes multiple stellar system formation. Results. In the purely hydrodynamical models, we always obtain (at least) binary systems. When more than two stars are present, their gravitational interaction triggers mergers until there are two stars left. The following gas accretion increases their orbital separation and hierarchical fragmentation occurs so that both stars host a comparable disk and stellar system which then form similar disks as well. We identify several modes of fragmentation: Toomre-unstable disk fragmentation, arm-arm collision and arm-filament collision. Disks grow in size until they fragment and become truncated as the newly-formed companion gains mass. When including magnetic fields, the picture evolves: the primary disk produces less fragments, arm-filament collision is absent. Magnetic fields reduce the initial orbital separation but do not affect its further evolution, which is mainly driven by gas accretion. With magnetic fields, the growth of individual disks is regulated even in the absence of fragmentation or truncation. Conclusions. Hierarchical fragmentation is seen in unmagnetized and magnetized models. Magnetic fields, including non-ideal effects, are important because they remove certain fragmentation modes and limit the growth of disks, which is otherwise only limited through fragmentation.


Introduction
Massive stars are very often located in multiple systems (Sana et al. 2012, Chini et al. 2012, Duchêne & Kraus 2013), and the importance of magnetic fields for massive star formation has been growing recently (Girart et al. 2009; for a review, we refer the reader to Tsukamoto et al. 2022), but the exact role of magnetic fields in the formation of multiple massive stars is not clear yet.Disk fragmentation has been shown theoretically (Adams et al. 1989, Shu et al. 1990, Bonnell & Bate 1994) and now observed thanks to high angular resolution facilities such as ALMA (e.g.Ilee et al. 2018, Johnston et al. 2020a, Johnston et al. 2020b).In addition to their role in the origin of massive protostellar outflows predicted by Kölligan & Kuiper (2018), Commerçon et al. (2022), Mignon-Risse et al. (2021a), Oliva & Kuiper (2023b) and detected with unprecedented resolution by Moscadelli et al. (2022) using masers, magnetic fields are expected to limit gas fragmentation and therefore are certainly important for multiple massive star formation (Añez-López et al. 2020).Numerical studies had to tackle both core and disk fragmentation, since their relative contribution to the formation of multiple systems remains to be evaluated, both in the low-mass and high-mass stellar regimes.So far, in the case of high-mass star formation, the interplay between radiation on small scales and magnetic fields have been shown to limit core fragmentation (Commerçon et al. 2011a, Myers et al. 2013) in the ideal magneto-hydrodynamical (MHD) regime, i.e. with a perfect coupling between magnetic fields and the dustgas-mixture.In the case of low-mass pre-stellar cores, this was shown in Commerçon et al. (2010), following the study by Hennebelle & Teyssier (2008) without radiative transfer.At the scale of turbulent molecular clouds, core fragmentation inhibition by magnetic pressure has been shown in a number of studies such as Federrath & Klessen (2012); see also the review by Padoan et al. (2014) and more recently the study of Lebreuilly et al. (2021).Incorporating ambipolar diffusion, a non-ideal MHD effect describing ion-neutral drift, the ability of magnetic fields to reduce disk fragmentation has been further shown in the study of Mignon-Risse et al. (2021b), oriented towards massive star formation.
Observationally, the presence and influence of magnetic fields either in massive star formation or in multiple protostellar systems has been increasingly constrained thanks to polarization measurements.The relative orientation of the field lines with respect to structures of interest (see e.g.Cox et al. 2022) has shown how the gas dynamics is coupled to the magnetic field topology.Furthermore, in a hub-filament structure (a structure which could be crucial in general for massive star formation, see the review by Motte et al. 2018), the magnetic fields strength has been found to be comparable to gravity (Wang et al. 2019, Añez-López et al. 2020).In the low-mass star formation context, the possible interplay between rotation and magnetic fields in gas fragmentation has been exposed by Galametz et al. (2020) with the submillimeter array (SMA).In addition, the magnetic field topology in a circumbinary disk around low-mass protostellar objects has also been revealed by Alves et al. (2018) using ALMA polarization observations and shown to be active at launching outflows (Alves et al. 2017).It suggests a long-lasting influence over the entire protostellar phase.
The present paper builds on recent numerical results.Disk fragmentation in a centrally-condensed system has been recently studied by Oliva & Kuiper (2020) and Mignon-Risse et al. (2023), hereafter Paper I. Those have focused on the properties of gaseous fragments forming in an accretion disk around a massive protostar modelled as a sink cell or sink particle, respectively, for which spiral arms play a prominent role.However, those studies have explored one case in which a disk forms around a central massive protostar and neglected magnetic fields.Hence, we propose to give a complementary point of view by focusing on the formation of multiple stellar systems (described by sink particles) and by studying the consequences of introducing magnetic fields.
Our choice of physical ingredients is well suited for studying protostellar disk fragmentation.First of all, non-ideal MHD is found to impact the disk formation by circumventing the socalled magnetic braking catastrophe (see Wurster & Li 2018 and references therein) and affecting subsequent disk properties that can impact disk fragmentation.In fact, the fragmentation critical length (the so-called Jeans length, Jeans 1902) depends on local pressures: thermal pressure and magnetic pressure mainly.It has been shown that ambipolar diffusion (drift between the ions and the neutrals) produces vertically thermally-supported disks rather than magnetically-supported as in the ideal MHD case (see e.g.Masson et al. 2016 in the low-mass and Commerçon et al. 2022 in the high-mass regime), of smaller size than in the hydrodynamical case.Let us note that another non-ideal MHD effect, Ohmic dissipation produces a thermally-supported midplane, where fragmentation should occur, in the massive protostellar context (Oliva & Kuiper 2023a).Moreover, including all three non-ideal MHD effects (ambipolar diffusion, Ohmic dissipation, Hall effect) in low-mass pre-stellar core collapse calculations, Wurster & Bate (2019) found less disk fragmentation in their magnetized models than in their purely hydrodynamical models.For these reasons, we include MHD with ambipolar diffusion in our magnetized models, as those are already implemented in the RAMSES code (Teyssier 2002, Masson et al. 2012).Finally, radiative transfer and its coupling to hydrodynamics are taken into account in order to accurately compute the gas temperature, essential for an accurate Jeans instability modelling.In that view, we use RAMSES with state-of-the-art physical ingredients: radiative transfer and MHD with ambipolar diffusion.
The combination of these physical ingredients and astronomical unit-scale resolution complements other studies led so far, in particular those carried out on molecular cloud scales and on which stellar feedback by massive stellar clusters is crucial (e.g.Grudić & Hopkins 2018, Ali & Harries 2019, Grudić et al. 2020).Indeed, disk scales -a prerequisite for modelling disk fragmentation -are hardly reached in present-day molecular cloud scale simulations (the finest resolution is, for instance, 1000 AU in Rosen et al. 2021, 100 AU in Mathew & Federrath 2021).In Bate (2018), disk scales are reached but magnetic fields are neglected.In Lebreuilly et al. (2021), disk scales are reached and magnetic fields with ambipolar diffusion are accounted for, but no massive star formation is reported (possibly due to limited statistics); moreover, the large (>100) number of stars and the environmental turbulence significantly complicate the detailed analysis when it comes to disentangling magnetic effects on the evolution of multiple stellar systems from other mechanisms, as mentioned above.By focusing on the small scale physics, we can reveal potential dominant magnetic mechanisms in the formation and evolution of massive stellar systems; if identified, those should be modelled in future large-scale simulation studies for which the properties of stellar systems (mass, multiplicity, separation...) are important.
Finally, this setup's properties are also complementary to other studies of collapse of massive pre-stellar cores.While Myers et al. (2013) and Mignon-Risse et al. (2021b) have studied fragmentation in turbulent clouds, we propose here to consider a non-turbulent massive pre-stellar core with a rotation profile favoring fragmentation (Oliva & Kuiper 2020).The ability of rotation to promote fragmentation has been extensively shown in the literature (e.g.Machida et al. 2005a, Wurster & Bate 2019).Observational results suggested that core rotation could play a role in dragging the magnetic fields (Beuther et al. 2020) and driving fragmentation even in turbulent cores (Girart et al. 2013, Zhang et al. 2014, with SMA).Our consideration is also meant to explore other initial conditions that those studied previously, possibly co-existing in nature, but also making easier the identification of fragmentation modes and of the role played by magnetic fields than in turbulent simulations, as mentioned above.
The paper, which presents results obtained via numerical simulations of massive pre-stellar core collapse, is organized as follows.First, we propose to study the influence of numerical resolution in Sec. 3, on disk fragmentation and stellar multiplicity, to identify the scales of interest for fragmentation and the runs that are converged enough for a deeper physical analysis.Then, based on these results, we study the influence of magnetic fields in Sec. 4. The next section presents the numerical methods.

Radiation-magnetohydrodynamical model
Simulations are performed with the RAMSES code (Teyssier 2002, Fromang et al. 2006).RAMSES is an adaptive-mesh refinement code (AMR) integrating the equations of radiation-MHD (RMHD).Non-ideal MHD is taken into account in the form of ambipolar diffusion (ion-neutral drift, Masson et al. 2012).For the radiative transfert part, we use an hybrid radiative transfer method (Mignon-Risse et al. 2020).In this method, the M1 method (Levermore 1984, Rosdahl et al. 2013, Rosdahl & Teyssier 2015) is employed to model the propagation and absorption of protostellar radiation, while the flux-limited diffusion (FLD, Levermore & Pomraning 1981, Commerçon et al. 2011b, Commerçon et al. 2014) is used to model dust-gas emitted radiation.The RMHD model consists in this set of equations where, ρ is the gas density, u is the velocity vector, P is the gas thermal pressure, λ is the flux-limiter for the FLD module, E fld is the FLD radiative energy, κ P, is the Planck mean opacity computed at the effective temperature of the primary star, c is the speed of light, F M1 is the M1 radiative flux, F L = (∇ × B) × B is the Lorentz force, φ is the gravitational potential, E T is the total energy which is defined as where is the specific internal energy), E M1 is the M1 radiative energy, B is the magnetic field, E AD is the ambipolar electromotive force, P fld is the FLD radiative pressure (which is not evolved as such but related to the other FLD moments, this is the FLD approximation), κ P,fld is the Planck mean opacity in the FLD module (computed at the local gas temperature), κ R,fld is the Rosseland mean opacity, a R is the radiation constant, P M1 is the M1 radiative pressure and Ė M1 is the injection term of the primary stellar radiation into the M1 module.For conciseness, we do not present the ambipolar electromotive force and the chemical network used to compute the resistivities (Marchand et al. 2016).For more details, we refer the reader to Mignon-Risse et al. (2021b).
Coupling between the M1 and the FLD modules occurs through the term κ P, ρcE M1 in the equation of temporal evolution of the internal energy, which is given by We use the ideal gas relation with the internal specific energy ρ = C v T , where C v denotes the heat capacity at constant volume.

Initial conditions
We use similar initial conditions as Paper I and Oliva & Kuiper (2020), and refer the reader to these papers for a complete setup description.We summarize its main characteristics here.The massive pre-stellar core has a mass M c = 200 M and radius R c = 20 625 AU = 0.1 pc, with density profile ρ(r) ∝ r −3/2 , resulting in a global free-fall time of 37 kyr, but a free-fall time of a few kyr in the innermost regions.Rotation is characterized by an angular frequency Ω(R) ∝ R −3/4 , where R is the cylindrical radius, ensuring a rotational-to-gravitational energy ratio independent of the radius in the midplane of the core, equal to 5%.The initial temperature is 10 K everywhere.Hydrodynamical simulations are performed with the Lax-Friedrich solver and magnetohydrodynamical simulations are performed with the HLLD solver (Miyoshi & Kusano 2005), as in Mignon-Risse et al. ( 2020) and Mignon-Risse et al. (2021b), respectively.
In the magnetohydrodynamical runs, a uniform vertical magnetic field is initialized (parallel to the rotation axis).The mass-to-flux ratio divided by the critical mass-to-flux ratio (Mouschovias & Spitzer 1976) at the core border is µ = 2, corresponding to a magnetic field strength B 0 = 0.43mG.This corresponds to a strongly magnetized pre-stellar core, but since the magnetic field is uniform, the magnetization at the center of the cloud is not as strong as µ = 2 (see Commerçon et al. 2022 andMignon-Risse et al. 2021b).Within a distance of ∼200 AU to the core center, the ratio of angular velocity to magnetic fields strength exceeds the critical value of (ω/B) crit = 3.19×10 −8 c s −1 yr −1 µG −1 with c s the sound speed given in km s −1 (Machida et al. 2005b).
Following from Paper I, a sink with mass 0.01 M is initially present at the center of the cloud for all runs.However, it is not fixed in space and other sink particles are allowed to form (see below).Qualitatively-similar results are obtained without such a central sink initially because of the mechanisms presented in Sec.3.3.This aspect is discussed in Sec.5.1.

Resolution and sink particles
The coarse resolution is level 5 (equivalent to a 32 3 regular grid) and we vary the finest resolution between levels 13 and 16, resulting in a physical resolution of 10 AU to 1.25 AU.Eight runs are considered : four hydrodynamical runs (prefix HYDRO-) and four magnetohydrodynamical runs (prefix MU2-).Suffix "LR" refers to "low resolution" (10 AU), "MR" to "mid-resolution" (5 AU), "HR" to "high resolution" (2.5 AU) and "VHR" to "veryhigh resolution" (1.25 AU).There is no constraint on the sink position neither on the maximal number of sinks to form, unlike the fiducial run of Paper I. As mentioned in Sec.2.2, one sink is initially present at the center of the box.On all AMR levels below the maximum, cells are refined so that the Jeans length is resolved by 12 cells (see Truelove et al. 1997).At the finest level, sink particles (Bleuler & Teyssier 2014) are introduced where a dense clump is found to be bound and Jeans unstable (more details in Commerçon et al. 2022).Sinks only interact gravitationally with the surrounding gas and other sinks, and merge if part of their radius (four cells here) overlap.Accretion onto sinks occurs if gas within the sink cells is above a density threshold given by 1.2 × 10 −13 (dx/10AU) −15/8 g cm −3 (Eq.11 of Hennebelle et al. (2020), see also Paper I), where dx is the finest physical resolution of the simulation.

Influence of spatial resolution
In this section, we investigate the impact of spatial resolution on massive stellar multiplicity, i.e. the number of sink particles, and related properties.

On the presence of accretion structures
Figure 1 shows the column density in the eight runs at the end of the run, ranging from t = 11 kyr for run MU2-VHR(about a quarter of the core free-fall time) to t = 20 kyr for the hydrodynamical runs.To begin with, in the hydro case we report spiral arms and disk-like structures in run HYDRO-LR, with a filament linking two stellar systems.However, some sink particles do not possess disk.Going to higher resolution, runs HYDRO-MR, HYDRO-HR and HYDRO-VHR suggest convergence with respect to those structures, and exhibit also the spiral arms after they grow.
In those runs, all sink particles have a disk.Nevertheless, the individual disks clearly visible in run HYDRO-HR (and HYDRO-VHR) appear in the place of circumbinary disks in run HYDRO-LR.We note that the filaments between the primary sink (already present in the initial conditions) and the secondary sink (i.e. the first sink that forms in the simulation, since the primary is already there) were created concomitantly with the secondary sink particle.We aim to check whether a given spatial resolution allows to capture spiral arm formation.In the hydrodynamical case, a 10 AU resolution is sufficient to capture structures prone to fragmentation, although circumstellar and circumbinary disks can be confused, and a 5 AU to 1.25 AU resolution appears to give a converged qualitative picture.
In all magnetic runs, a pseudo-disk is present.It is characterized by an enhanced density compared with the surrounding medium, a magnetic pressure dominating over thermal pressure, and infall motions.Its size increases with time.In Fig. 1 it is clearly visible for run MU2-MR and MU2-HR with a ∼1000 AU radius and for run MU2-VHR with a ∼700 AU radius.We will not focus on pseudo-disks (Galli & Shu 1993) but rather on rotationally-supported disks in which the centrifugal acceleration ensures radial equilibrium (simply called "disks" for the rest of the paper).As can be seen in the bottom-left panel of Fig. 1, in run MU2-LR there is no rotating disk larger than the sink particle size at t = 10 kyr.Consequently, there is no spiral arm or filament connecting the two sinks or their surrounding accretion structures.In this run, disks appear after about half of the simulation time.In runs MU2-MR, MU2-HR and MU2-VHR, we find individual disks, spiral arms and a filament linking the two disks, as in the hydrodynamical case.In run MU2-MR, individual disks are resolved but over the integration time we do not observe the formation of new spiral arms after a binary has formed.On the opposite, we report that spiral arms keep forming regularly in runs MU2-HR and MU2-VHR, although they are hardly distinguishable on a single time snapshot like Fig. 1.In run MU2-VHR, a third sink particle forms in one of the individual disks after spiral arm collision (reminiscent of Mignon-Risse et al. 2021b).The common accretion disk, where the third sink formed and is still embedded at the end of the run, is the closest structure to what could be called a circumbinary disk in any of these simulations.For the magnetic case, under these conditions we find that a spatial resolution of 2.5 AU is sufficient to resolve the structures prone to fragmentation -which does not mean that convergence is reached, see Sections 3.2 and 3.3.In particular, disks (more easily identifiable by their column density structure than by rotation arguments, see Mignon-Risse et al. 2021b) are smaller than in the hydrodynamical case (see Sec. 4.4;and smaller than in Mignon-Risse et al. 2021b).Consequently, spiral arms and fragmentation induced by spiral arm collision appear at higher resolution than in hydrodynamical runs, and fragmentation might be missed in simulations that do not reach such spatial resolutions requirements.
Let us note that the cavities visible in the bottom-left panel of Fig. 1 (run MU2-LR, also visible later in the other magnetic runs) are reminiscent of those reported in Hennebelle et al. (2020), Commerçon et al. (2022), andMignon-Risse et al. (2021b).These are likely due to the interchange instability (see Appendix B of Mignon-Risse et al. 2021b).

Jeans length
After we have addressed the spatial resolution needed to capture structures of interest, let us see how resolution impacts the Jeans length.The (thermal) Jeans length is defined as where c s = γP/ρ is the local sound speed, with γ = 5/3 (in our calculations) the adiabatic index.The magnetic Jeans length L J,mag is defined in a similar fashion but using c 2 s + v 2 A , where v A = B/ 4πρ is the Alfvén speed, instead of c s (more details in Sec. 4 on the comparison between L J and L J,mag ).For an alternative approach to computing the Jeans length in the MHD case, we refer the reader to the Appendix of Myers et al. (2013), who also discuss the anisotropy of magnetic effects.Note, however, that the difference between the L J,mag they derive and our simple formulation is minor.
Figure 2 shows the minimal Jeans length as a function of the maximal AMR level (a higher AMR level means a finer resolution by a factor of two).Hence, one would like the Jeans length slope to be less steep than a l −0.5 max slope and ideally to become flat, indicating that the Jeans length becomes increasingly well sampled as the resolution increases.
It is visible in Fig. 2 that for certain AMR levels, the smallest L j actually decreases faster than a l −0.5 max slope.It is particularly visible from l max = 13 to l max = 14 (runs HYDRO-LR and HYDRO-MR).Indeed, we find that the maximal density reached increases with the resolution (not shown here for conciseness).Consequently, a higher maximum AMR level does not always result, in those simulations, in a better-sampled Jeans length.Nevertheless, we checked that the minimal Jeans number, that is the number of cells sampling a Jeans length, is generally larger than 8 at the maximum level, so that artificial fragmentation should be avoided (see the discussion in 5.1).There is no   1.In the magnetic case, we compute the thermal (L J , see Eq. 3) and magnetic (L J,mag ) Jeans length.An additional AMR level gives a finer resolution by a factor of two.
clear trend in the magnetic runs for a Jeans length generally decreasing less steeply than l −0.5 max .Nevertheless, it is the case from l max = 14 to l max = 16 for the hydro runs, corresponding to a resolution of 2.5 AU (run HYDRO-HR) to 1.25 AU (run HYDRO-VHR).
In summary, checking the Jeans length sampling at a given resolution only indicates that fragmentation, if reported, is not artificial (Truelove et al. 1997, see also Federrath et al. 2011 and the discussion in Sec.5.1).It does not mean that the Jeans length converges towards a unique value because it depends on the density and temperature, which depend on the maximum resolution.

Multiplicity: number of sink particles
Figure 3 shows the number of sinks as a function of time.We distinguish three phases: initial fragmentation (increasing number of sinks) due to Toomre-unstable disk fragmentation, mergers (number of sinks decreasings towards a value of 2), and secondary fragmentation.In runs HYDRO-VHR and MU2-VHR, only two sinks are formed after the first fragmentation phase so there is no merger phase.
First, it can be seen that the sink multiplicity increases before t = 6 kyr from 1 (the initial sink) to 2 − 5.This is the initial fragmentation epoch, which is highly impacted by the formation of a disk with a density bump on top of it, prone to fragmentation as it is the most Toomre-unstable region (see Paper I).This structure fragments rapidly (in less than one orbital period, Norman & Wilson 1978) into four pieces in runs HYDRO-LR and MU2-LR, likely triggered by the Cartesian grid.As already mentioned in Paper I, these grid effects dominate for smooth initial conditions; other possibilities would be to choose perturbed (e.g.Commerçon et al. 2008, Kuruwita et al. 2017) or turbulent initial conditions (e.g.Gerrard et al. 2019, Kuruwita & Federrath 2019, Mignon-Risse et al. 2021b) to seed the instability.In the other runs, we nevertheless observe the development of non-axisymmetries similar to spiral arms around the primary sink, where additional sinks form.In runs MU2-MR, MU2-HR and MU2-VHR, the disk is affected by magnetic forces and is elongated into a bar, similarly to the structure reported in Machida et al. (2005a) following Machida et al. (2005b).Hence, in those cases, the disk instability has unlikely been seeded by the grid.
During a second phase, the sink multiplicity decreases as sinks merge.The central sink was initially in an unstable equilibrium as it was surrounded by massive companions.Hence, any shift (e.g.caused by numerical errors or any asymmetry) leads this sink to leave the center.We understand the drop in sink multiplicity to be due to gravitational interactions between  sinks.Because of these interactions, the sinks' angular momentum (computed with respect to the origin) changes, favoring the migration of some sinks towards the center, where the primary sink is initially located, until they merge.Eventually, all runs show a binary system at some point between t = 7.5 kyr and t = 10 kyr.
In the hydrodynamical runs, the sink multiplicity increases again during the third phase, that we refer to as secondary fragmentation.At this stage, each sink of the binary system has developed an accretion disk.New sinks form at the extremity of spiral arms or following the collision between a spiral arm and a filament or another spiral arm.The final number of sinks is between 4 (mid-, high-and very-high resolution) and 7 (lowresolution).According to the value of the Jeans number being larger than 8 in the HYDRO-LR run and most fragmentation structures (circumbinary disk, filaments) being captured, fragmentation is physical rather than numerical.Nonetheless, the excess of sinks in run HYDRO-LR can be explained by the lack of convergence regarding the resolution of structures.
In the MU2 runs, the number of sinks at the end of the simulations is between 2 and 3.In run MU2-LR, we have shown above that no clear disk structure forms around the two sinks (bottomleft panel of Fig. 1), hence we do not expect the formation of additional sinks.This system has reached a quiescent state.As shown in run MU2-VHR (bottom-right panel of Fig. 1), secondary fragmentation is not fully prevented in the presence of magnetic fields, as we observe that a third sink particle has formed from spiral arms collision.It occurs in the close neighbourhood of an existing sink, at a distance of 19 AU.This could not have occurred in run MU2-HR with such a small orbital separation, because of the two sinks accretion radii (10 AU each).
Figure 5 shows the orbital separation between the sinks that will eventually become the two most massive sinks as a function of time.We find it to correctly converge with resolution.It will be studied in more details in Sec.4.2.
Overall, we observe hierarchical fragmentation, with each sink possessing its own disk that can lead to additional fragmentation, in the hydrodynamical case as well as in the magnetic case.We find that early fragmentation is compensated by enhanced stellar mergers following gravitational interactions and radial migration.This suggests that high-multiplicity systems born from disk fragmentation likely formed out of hierarchical fragmentation, rather than from a single disk.

Conversion of gas into stars
Figure 6 shows the sum of all sink masses as a function of time.We observe a satisfying convergence in the HYDRO run at the end of the simulation, within ≈20% with respect to the HYDRO-VHR run.We nevertheless observe large variations at low resolution (run HYDRO-LR).Let us recall that the density threshold for accretion obeys a scaling relation (Hennebelle et al. 2020), which appears here to lead to self-consistent results.In the magnetic runs, we find an agreemeent on the total mass within ≈20% as well at the end of the runs.
To sum up, the total mass is nearly independent of spatial resolution, especially at late times.

Influence of magnetic fields
In this section, we investigate the impact of magnetic fields on the fragmentation modes and on the properties of massive multiple stellar systems.

Multiplicity and modes of fragmentation
To begin with, we observe two common fragmentation mechanisms between hydrodynamical and magnetic cases, namely Toomre-unstable disk fragmentation (with a transition through an elongated bar in high-resolution magnetic runs) and arm-arm collision.Moreover, in all runs we observe sink mergers and the presence of a binary system after the early fragmentation epoch (Fig. 3).The difference is found in the secondary fragmentation phase, which does not exist in most of the magnetic runs, except for a sink particle born out of spiral arm collision in run MU2-VHR, indicating that fragmentation is not totally suppressed by magnetic fields.Meanwhile, arm-filament collision, which is found to trigger fragmentation in the hydrodynamical runs, is not observed in any magnetic runs.Indeed, filaments linking binaries (as also observed by, e.g.Sadavoy et al. 2018in IRAS 16293-2422) are more diffuse when magnetic fields are included, as can be seen in Fig. 1, which may be due to additional magnetic pressure (as can be the case for larger-scale filaments, see e.g.Federrath & Klessen 2012, Gutiérrez-Vera et al. 2022) in those regions of smaller density than in the disks, where pressure support is mainly magnetic.By broadening the filaments, magnetic pressure prevents, indirectly, fragmentation induced by arm-filament collision.
As shown in Fig. 2, fragmentation is not directly suppressed by magnetic pressure support -although, as discussed above, it prevents indirectly the fragmentation through arm-filament collision -, which has a minor impact on the Jeans length.Indeed, we find that regions of small Jeans length are as numerous in the hydrodynamical as in the magnetic cases.Similarly, within the magnetic runs, the minimal Jeans length and magnetic Jeans length are roughly equal because thermal pressure largely dominates over magnetic pressure.L J and L J,mag converge towards the same value at high resolution, while in run MU2-LR no thermallysupported disks were formed at that time.When looking at the Jeans length more generally, we find that L J and L J,mag differ where the Jeans length is already large and therefore where the gas is already stable against fragmentation.Hence, magnetic pressure is not responsible for these differences in sink multiplicity.More than that, if the scales on which ambipolar diffusion is efficient are not sufficiently resolved, magnetic field piles-up which leads to overestimating the stabilizing role of magnetic pressure.
Overall, including magnetic fields conserve two modes of fragmentation: Toomre-unstable disk fragmentation and fragmentation associated to arm-arm collision.However, no fragmentation induced by arm-filament collision is reported.

On the variations of the orbital separation in binaries
As shown in Fig 5, the orbital separation between the stars that will become the most massive ones, a(t), is found to increase nearly linearly in all runs on timescales of kyr.Quasi-periodic variations on timescales of the order of a kyr are linked to eccentricity (as in, e.g., Park et al. 2022).Variations on shorter timescales, for example between 15 kyr and 18 kyr in run HYDRO-VHR, are due to orbital motions with other companions, born from secondary fragmentation.In the following we study, on the one hand, the impact of magnetic fields on the orbital separation and, on the other hand, the two possible origins of this outspiral motion: gravitational torques from the gas and gas accretion.
Let us first present the equation describing the orbiral separation evolution of the binary.Under the approximation of equalmass binaries (which is the case most of the time in the hydro runs when two stars are present, and in the magnetic runs as well) and centrifugal equilbrium we have a = 16 j 2 /(GM), where j is the specific angular momentum of the binary -also equal to that of each component of the binary -and M is the total mass of the binary.Then, the temporal variation of the orbital separation can be derived as so that the condition for an increasing orbital separation is (the formulation is similar to, e.g.Tiede et al. 2020) When computing the evolution of the sinks angular momentum and mass over timescales of kyr, on which we observe an increasing separation (Fig 5), we find this condition to be met, which is physically consistent with the outspiral we observe.

Impact of magnetic fields on the separation
Do magnetic fields have any impact on the increasing separation?We observe that the orbital separation is reduced by a factor of two in binaries when magnetic fields are present.Since a(t) is nearly linear, it takes the form a(t) = αt where α is a constant, with α HYDRO ≈ 2α MU2 .Hence, the quantity ȧ/a (where the dot indicates the time derivative) is equal in both cases.In other words, the orbital separation evolution is unaffected by magnetic fields.
Magnetic fields only interact with sinks via their effect on the gas.To address the possible angular momentum removal by magnetic braking, we compute the angular momentum at two distinct epochs, on cubes encompassing the binary system.
On the one hand, we focus on early times and small scales around the center of the cloud.Computing the angular momentum in a cube of side 800 AU centered onto the grid origin at 3.5 kyr (typically the time when the secondary sink forms), we find that the specific angular momentum is reduced by 30% in the magnetic case.Hence, by the time the secondary sink forms, part of the angular momentum has been removed by magnetic braking.This impact of magnetic fields explains the difference in the initial orbital separation of binaries.
On the other hand, we focus on later times, which means larger scales because the orbital separation increases.Computing it in a cube of side 2000 AU centered onto the grid origin at t = 5 kyr and 10 kyr, we do not find significant differences between the hydrodynamical and the magnetic runs.How to explain the lack of efficiency of the magnetic braking (with respect to the center of mass, which is the center of rotation of the system)?This mechanism's strength is proportional to the toroidal and poloidal components of the magnetic field, the former developing from coherent (differential) rotation.Between t = 5 kyr and 10 kyr each star has a disk but no circumbinary disk.The consequence is twofold.First, as illustrated in the top panel of Fig. 7, the magnetic field strength outside the individual disks is smaller (by a factor of a few to a factor of ∼100) because the density there is smaller.Second, portions of the gas rotating coherently around the binary indeed generate a toroidal field; locally, B φ / B ≈1, as shown in the bottom panel of Fig. 7.However, these portions are too small, and the local magnetic field strength is too small for the overall magnetic braking to slowdown the rotating gas on timescales shorter than the time it takes for this same gas to reach any of the two individual disks.In those disks, the magnetic braking with respect to the center of mass is negligible, as expected.Hence, the binary system's angular momentum cannot be transported efficiently by magnetic fields.
Therefore, since ȧ/a is equal in equal in both cases, this suggests that the value of a(t) is inherited from its initial value, which differs in the hydrodynamical and in the magnetic cases because of magnetic braking.Hence, magnetic fields reduce the initial separation but not its further evolution.

Gravitational torques contribution
Now, one should investigate the contribution from gravitational torques and gas accretion onto this increasing separation.To probe the former, we compute the sink specific angular momentum j and the gravitational torques from the gas onto the sinks, noted rg φ , which accounts for the Plummer softening specific to gas-sink interactions.The specific angular momentum glob-ally increases with time -with values of the order of 10 18 cm 2 /s, and slope of 10 6 cm 2 /s 2 .This is expected because the orbital separation increases with time whereas the binary mass, which can only increase, promotes orbital shrinking.Theoretically, the gravitational torque must be negative when corresponding to the disk part lagging behind a sink and positive when it is in front of it (see e.g.Tiede et al. 2020).If the gas distribution was perfectly symmetric with respect to the binary axis, the sum of the two should be zero.Meanwhile, spiral arms, companions (when present) and low-density cavities due to the interchange instability (in the magnetic models) should make it deviate from zero.This is supported by our results: the total gravitational torque rg φ , shown in Fig. 8 for runs HYDRO-HR and MU2-HR, starts showing an oscillatory-like behaviour, at about the orbital period, when spiral arms develop and a companion forms (in the hydro runs), as those orbit around the binary and alternatively contribute to a positive and negative torque.On the opposite, low-density cavities always lag behind the sinks; their presence should translate into a positive torque on the sinks; since we find the torque to oscillate around zero, this contribution appears to be negligible compared to the spiral arm contribution.We note, however, that the oscillation amplitude is larger in the magnetic runs than in the hydro runs.This is likely due to the disks being smaller in the magnetic case: the distance to the entire spiral arm structure is smaller and the associated gravitational torque, which decreases quadratically with the distance, is larger.Overall, the envelope of the oscillations in rg φ has an amplitude of ≈10 5 cm 2 /s 2 .This is one order of magnitude too small to explain the slope of the specific angular momentum evolution.Hence, only gas accretion can explain the increasing orbital separation.

Gas accretion
As shown in Eq. 5, gas accretion can, depending on its specific angular momentum and mass, both contribute to increasing or decreasing the separation.However, in this run the general trend is an increasing separation.It is therefore primordial to understand why and to see whether it is purely linked to our setup.
Integrating both sides of Eq. 5 we get d j/ j > dM/2M.Then, we integrate the left-hand side from j 0 to j and the right-hand side from M 0 to M , where j 0 (resp.j is the sink specific angular momentum prior to (resp.after) accretion and M 0 (resp.M ) is the sink mass prior to (resp.after) accretion, yielding j / j 0 > √ M /M 0 .To go further, we express j and M as functions of j 0 , M 0 and the accreted gas angular momentum and mass labelled j gas and M gas , as and M = M gas + M 0 , that we replace in the previous inequality.
After some simple algebra, as long as the gas mass increment is negligible compared to the sink mass (M gas M 0 ) so we neglect second-order terms (M gas /M 0 ) 2 , this leads to the following criterion for an increasing orbital separation: where j 0 = √ aGM from radial equilibrium between the sinks.This is a reformulation of the trend found in the pioneer studies of Bate & Bonnell (1997), Bate (2000).Putting in typical values of separation and binary mass from the simulations, we find that gas accretion increases the separation of the binary if (8) We will use those fiducial values for the following reasoning.A relevant question is the following: does the gas, initially located at a given radius and with a given angular momentum (somewhat arbitrary in those simulations), possess enough angular momentum to fulfill condition 8 and therefore increase the orbital separation of the binary, at the time of interest?For our initial rotation profile, the fiducial value of Eq. 8 corresponds to gas initially located at radii > 10 3 AU.The follow-up question, to ensure the validity of the reasoning, is: does the gas initially located at this radius had enough time to reach the binary ?Yes, based on its associated free-fall time.The scenario of an orbital separation driven by gas accretion in our simulations is thus consistent.
The next natural question would be how this value depends on the initial profile.Although this comparison should be taken with caution because a different rotation profile would likely result in a distinct evolution of the system (e.g.Bate 2000), for different initial rotation profiles such as the slow and fast models presented in Commerçon et al. (2022) the gas should be initially located at ∼10 4 AU to drive the outspiral of such a binary.
To sum up on this subsection dedicated to the orbital separation evolution, we find that gas accretion is mainly responsible for the binary outspiral and we can link such a trend to the initial rotation profile of the core.We report the influence of magnetic fields in only one aspect: they remove the angular momentum in the innermost regions of the cloud at early times, setting a smaller initial orbital separation for binaries than in the hydrodynamical case.However, they play no significant role in the following evolution of the orbital separation.

Conversion of gas into stars
As visible in Fig. 6, in the hydrodynamical case, the total mass of all sinks is ∼30 M near the end of the simulation (t = 17.5 kyr) against ∼15 M in the magnetic case.Nonetheless, the mass growth is similar within 30% up to ∼10 kyr.
In the magnetic case the mass increases with a mean accretion rate slightly smaller than 10 −3 M yr −1 .In the hydrodynamical case, the mean accretion rate is Ṁ≈10 −3 M yr −1 up to the secondary fragmentation epoch (indicated by the circles in Fig. 6).After that, it is Ṁ≈3.5×10 −3 M yr −1 .Hence, the change in the mean accretion rate occurs just before secondary fragmen-tation occurs.This is consistent with the fragmentation being linked to spiral arms.Indeed, they carry angular momentum outwards, allowing the central object to accrete.Though, it does not explain why the accretion rate remains at this higher value.A possibility would be the initial density profile ρ(r) ∝ r −1.5 making more gas available for accretion at later times.However, this would not explain the difference between the hydrodynamical and the magnetic case.Moreover, this trend is not observed in the run with a single sink of Paper I (neither in the study of Oliva & Kuiper (2020) with a similar initial density profile).Hence, we do not attribute the enhanced accretion to the initial density profile.Another possibility is that the presence of new close (∼100 AU) companions also triggers the formation of spiral arms on shorter timescales (the orbital timescale between the two close components), encouraging the accretion.Similarly, the sinks born during secondary fragmentation form their own disk with its own spiral arms.A disk and a sink can perturb the secondary disk either via their spiral arms and via tidal forces (see e.g.Savonije et al. 1994), respectively.As a consequence, the other disk grows spiral arms and accretion would be enhanced.
In the following, we investigate the role of sink multiplicity on the accretion rate using this figure and look for patterns.Figure 9 shows the instantaneous accretion rate as a function of time.In the magnetic case, and in particular for run MU2-HR, we observe accretion peak-like features after 10 kyr that could be periodic.We extract their exact time by selecting the time of accretion rates that are local maxima and that are greater than the mean accretion rate in the simulation.We find that the time between two maxima is between 0.75 kyr and 1.0 kyr.This interval matches the semi-orbital period that can be extracted from the sink motions (also visible in the orbital separation plot, Fig. 5, thanks to the non-zero eccentricity of the system).This is consistent with the m = 2 mode expected to develop when the binaries are not too close with the spiral wave pattern speed being the orbital frequency (Savonije et al. 1994).In order to check for this periodicity, we have computed the Fourier power spectrum of the accretion rate after 10 kyr, as shown in Figure 10.We find an enhanced signal at frequencies between 0.7Ω orb and 2Ω orb .In particular, we find no significant periodicity above 2Ω orb .Nonetheless, it is difficult to draw firm conclusions from this plot, even though it suggests a plausible periodicity in the signal.Indeed, we know that the system's orbital frequency is not constant with time, there may be a non-periodic part of the signal that is not strictly constant either (see Fig. 9) and only a few (≈3) orbital periods are captured after 10 kyr within the simulation time.This shows, however, than a link between multiplicity and accretion rate is plausible, in addition to being justified theoretically.
In the hydrodynamical runs, we find that the timescale between such accretion peaks mainly lies in the 0.1 − 0.4 kyr interval.This is also consistent with the previous scenario, since a larger (> 2) sink multiplicity involves more than one characteristic frequency, and in those runs there are also close binary systems with higher orbital frequencies than in the magnetic runs.Because there are several frequencies involved, the periodicity in the accretion rate are more difficult to extract and much less visible (Fig. 9).
Our simulations do not show a direct influence of magnetic effects on the accretion rate.However, they suggest a possible modulation of the accretion rate at the orbital period due to the gravitational influence of a companion.Nevertheless, the system on hand -more than two components, a variable orbital period with time, a binary embedded in a collapsing core -is not well suited to extract such quasi-periodic features.

Individual disk radius
In this Section, our aim is to estimate the radius of the individual disks surrounding sink particles and the underlying process regulating its size, emphasizing the differences observed between hydrodynamical and magnetic models.Several diagnostics can be used, in principle, to derive the radius of a disk.
First, we consider using kinematic signatures, as is done observationally with velocity-positions diagrams (e.g.Murillo et al. 2013, Tobin et al. 2020).One way to do it in the simulation is to compute the deviation with respect to Keplerianity (as done in Paper I) on a cell-by-cell basis, and to reduce this information to a single scalar by computing the mean or median azimuthal value and looking for a transition radius.However, we find the Keplerianity to vary significantly between the cells as it is greatly affected by the stellar companions.The same limitation is encountered when using the disk criteria presented in Joos et al. (2012) andused, e.g., in Mignon-Risse et al. (2020).Therefore, we choose not to focus on a kinematic definition of the disk, although our results still point to structures dominated by rotation.
Observationally, dust continuum fitting is used as an alternative to kinematic signatures of disks (e.g.Patel et al. 2005).As in Mignon-Risse et al. (2021b), we propose to use the surface density as a criterion to define the disk and thus, the disk radius.Integrating the surface density over of depth of 200 AU perpendicular to the disk plane (this is already sufficiently larger than the typical disk height, which is less than 50 AU), we get a surface density Σ(y, z), a two-dimensional function, as it depends on the spatial coordinates in the disk plane.The largest surface density is found to be around Σ∼10 26 cm −2 , decreasing to Σ∼10 23 cm −2 at cylindrical radii of 400 AU.We reduce Σ(y, z) to a 1D function by computing its azimuthal median, taking one sink's location as the coordinates origin; we obtain Σ median (r cyl ).We use a threshold criterion to identify the disk radius as the cylindrical radius within which the surface density Σ median (r cyl ) exceeds a critical value taken to be 10 25 cm −2 .Using this diagnostic, we find that the disk grows in size, reaching 100 to 150 AU size in hydrodynamical runs before secondary fragmentation.After secondary fragmentation, in which the newly-formed companion starts forming its own disk, the primary disk radius decreases and remains between 60 AU and 100 AU in runs HYDRO-HR and HYDRO-VHR.During this epoch, the companion gains mass rapidly and the disk radius is compatible with disk truncation, e.g. about 0.3 of the orbital separation (Papaloizou & Pringle 1977, Paczynski 1977).Meanwhile, in the magnetic runs, the disk radius remains generally smaller than 100 AU even though no companion forms, as shown in Fig. 11, and is compatible with magnetic regulation (Hennebelle et al. 2016).Varying the threshold value to be 5 × 10 25 cm −2 decreases the obtained disk radius by about 20 AU.
Finally, Commerçon et al. (2022) found that the protostellar disks formed in presence of ambipolar diffusion have β = P/P mag > 1, where P mag = B 2 /(8π) is the magnetic pressure -more recently, it has been shown in simulation incorporating Ohmic dissipation instead of ambipolar diffusion as well (Oliva & Kuiper 2023a).Following the former study, Mignon-Risse et al. (2021b) argued that this could serve as a way to measure the disk radius, especially when the dynamics and multiple stellar components make difficult the comparisons between rotation profiles and Keplerian ones, which is the case here.Here, we indeed found that the closest region to the sinks have a β > 1 and that the transition with β < 1 is located near the surface density transition setting the disk edge.With a spatial resolution smaller than 5 AU we derive a disk radius, based on the plasma beta, of   10.Fourier power spectrum of the accretion rate in run MU2-HR after 10 kyr.The x−axis is the Fourier transform frequency normalized by the approximated orbital frequency Ω = 1/2.5 kyr −1 .Data (instantaneous accretion rates) have been linearly interpolated in order to have a uniformly-spaced sample to perform the fast Fourier transform.The gray band indicates the 0.7 − 2Ω orb interval.The low-frequency signal does not point toward a physical periodicity but to the sample time which is only 6 kyr long.More interestingly, we find no significant peak at frequencies higher than 2Ω orb .This suggests that the accretion rate is modulated by frequencies of the order of the orbital frequency Ω orb .
50 AU to 80 AU (see Fig. 11).This disk radius value is comparable to the surface density estimation using the largest surface density threshold among the two values explored here, and is nearly constant when varying the resolution.
To conclude, we found the column density map and the plasma beta to give a relatively compatible (±20 AU depending on the surface density) and converged (for resolution finer than 5 AU) estimate of the disk radius.Disks are roughly twice smaller in the magnetized case compared to the hydrodynamical case before secondary fragmentation.They are compatible with magnetic regulation, while in the hydrodynamical case they grow until fragmentation occurs and the newly-formed companion truncates the primary disk.Hence, in the hydrodynamical case, this result could be interpreted as the disk size not being regulated (it gathers material reaching its centrifugal limit) or being regulated by fragmentation.

Comparison with previous works
In the simulations presented here, the refinement criterion for the AMR levels below the maximum is a Jeans length resolved by 12 cells.While in the original work of Truelove et al. (1997), a number of 4 cells was mentioned (and widely used since then, see e.g.Kuruwita et al. 2017, Gerrard et al. 2019; other authors opted for 8 cells, see e.g.Masson et al. 2016, Rosen et al. 2016).To achieve convergence with SPH methods, Commerçon et al. (2008) showed that in the hydrodynamical case the Jeans length should be sampled by at least 15 cells.In order to resolve minimal dynamo amplification in self-gravity MHD simulations, Federrath et al. (2011) showed that 30 cells per Jeans length were required.Because the present initial conditions are not turbulent and the present disks do not develop turbulent states, we chose to keep the Jeans length resolution to 12 cells to save computational time.
Bate (2018) studied the morphology of protostellar disks from molecular clouds simulations with SPH methods, neglecting magnetic fields.In our high-resolution hydrodynamical simulations (e.g.HYDRO-VHR), we obtain a quadruple system that is somehow similar to that depicted in their Fig. 4 (sinks 41, 89, 76, 83; see also Sigalotti et al. 2018).They have in common the following properties: this is a point-symmetric system composed of two tight pairs surrounded by gas.Additionally, the two pairs are linked together by gas filaments continuously fed by the spiral arms belonging to protostellar disks.Nevertheless, they report circumbinary disks around the tight pairs while we do not.Instead, we find each star in the quadruple system to have a protostellar disk around it.This could be a matter of angular momentum budget of the surrounding gas, since high angular momentum material is more prone to form circumbinary disks (Bate & Bonnell 1997).Finally, in our simulations, the youngest sinks formed from the interaction between the spiral arm belonging to the oldest sinks' disks and the filament linking them; such fragmentation triggered by arm-filament collision seems to occur in their study as well.We find, however, that the inclusion of magnetic effects suppresses this fragmentation mode as magnetic pressure reduces the density contrast around the filament.
In this study, non-ideal MHD was included in the form of ambipolar diffusion.Hence, let us compare the properties of the multiple stellar systems formed in the magnetized runs with those of the binary systems reported in Mignon-Risse et al. (2021b), for which the included physics was identical.The initial conditions are quantitatively different (core mass, rotation...) but the qualitative difference comes from the inclusion of initial velocity dispersion to mimic turbulence while here all the angular momentum originates from differential rotation.Unlike Mignon-Risse et al. (2021b), where a circumbinary structure is present under particular turbulence level and magnetic field strength, we find no circumbinary structure in the binary systems formed here, except in run MU2-VHR where a sink is still orbiting within the accretion disk it was born in.As in Mignon-Risse et al. (2021b), we do find mass ratios of the order of unity in binary systems.This fits the "high angular momentum core" case, while a "low angular momentum core" is predicted to result in unequal-mass binaries by Bate (1997), Bate & Bonnell (1997).The triple system in run MU2-VHR has masses 3.7 M , 3 M and 1.4 M at the end of the run, with a mass ratio close to one before the third sink formed.
Preliminary works in low-mass star formation with ambipolar diffusion reported disk sizes reduced with respect to the standard hydrodynamical case (Masson et al. 2016, Hennebelle et al. 2016).As in the studies of massive isolated collapse of Commerçon et al. (2022), which forms a single-star system, and Mignon-Risse et al. (2021b), which forms a binary system when initial turbulence dominates over magnetic fields, disks are found to be smaller in the non-ideal MHD case than in the hydrodynamical case, based on the column density maps (and on the plasma beta in the magnetic case).Recently, Lebreuilly et al. (2021) confirmed this trend on large samples in the context of collapses within massive (1000 M ) clumps, aiming at describing the population of protoplanetary disks, although this work was mainly oriented towards low-mass star formation.Gerrard et al. (2019) studied the influence of an initialized turbulent magnetic field on disk fragmentation and outflow launching during the collapse of a low-mass pre-stellar core, under the ideal MHD approximation.They claim that a turbulent magnetic field leads to a more isotropic magnetic pressure distribution than a non-turbulent one, ending up in a reduced magnetic pressure gradient and more disk fragmentation.In fact, it can be expected from their assumptions (ideal MHD) that their disks are supported by magnetic pressure; while the disks formed here, in the presence of ambipolar diffusion, are thermally-supported, so that changes in the magnetic pressure distribution should not have a major influence on the disk fragmentation.
Magnetic fields may also play a role in the orbital separation of binaries.Recently, Harada et al. (2021) reported a numericaland-analytical study of the impact of magnetic braking on massive close binaries, in the ideal MHD framework.They show that the orbital separation is larger in the hydrodynamical case than when magnetic fields are introduced (their Eq. 6 with the non-specific angular momentum is equivalent to Eq. 5 here), because of magnetic braking and outflows.Our findings qualitatively agree with this result (see Fig. 5), as we find an orbital separation reduced by a factor of ∼2 when magnetic fields are present.However, we attribute it to the initial separation being smaller, because of early magnetic braking (when both stars were part of a coherent structure), while the further evolution of ȧ/a is comparable; on the opposite, their results only apply to the urther evolution of ȧ/a, since the initial separation is fixed.Moreover, quantitatively, we also find orbital separations larger than 100 AU in the magnetic case, while they do not.Our results differ because they assume a circumbinary structure able to extract angular momentum from the close binary, while we witness the formation of two separate disks except in the secondaryfragmentation formed binary system of run MU2-VHR.Such differences are partially attributed to the different choices of initial mass and angular momentum budget, which, as we have shown, drive the evolution of the orbital separation in binaries formed from initial fragmentation.Our initial conditions are indeed different (magnetic field strength, mass and angular momentum budget).The MHD version is different as well: non-ideal here, ideal in their work.Ideal MHD is found to overestimate the magnetic field strength in the densest (ρ 10 −15 g cm −3 ) regions (Masson et al. 2016), impacting both disk formation and outflows (Commerçon et al. 2022).Nevertheless, the present study supports the interpretation of Harada et al. (2021) that magnetic fields, in the form of magnetic braking, could play a prominent role in the formation of massive close binaries -at the condition of a coherent rotation structure such as a circumbinary disk -, at least initially.Further in time, however, gas accretion equally widens binaries in the hydrodynamical and in the magnetized cases, in the absence of a circumbinary accretion structure.
The formation of a circumbinary disk structure around a lowmass protobinary system in non-turbulent and turbulent cores has been studied by Kuruwita & Federrath (2019), including ideal MHD.They report orbital shrinkage and attribute it to magnetic braking on the infalling gas.Alternatively, this could be due to the sinks not being at centrifugal equilibrium when they form, since they arise from initial perturbations, shrinking their separation until they reach a stable orbit.To understand whether the final evolution of their orbital separation is consistent with gas accretion (that is, no outspiral), we put into Eq.8 the final mass (∼0.5 M ) and orbital separation (∼20 AU) of the protobinary system in their non-turbulent run.This gives a minimum gas angular momentum of 2 × 10 20 cm 2 s −1 to drive the binary outspiral through gas accretion; meanwhile, their protostellar core has a maximal specific angular momentum of ≈10 18 cm 2 s −1 at the core edge.Hence, their pre-stellar core angular momentum budget is insufficient to drive the binary's outspiral.Because they consider relatively small turbulent velocities (0.02 km s −1 and 0.04 km s −1 ), the angular momentum carried out by turbulence is still insufficient to drive the binary's outspiral and this naturally results in a similar final orbital separation in their non-turbulent and turbulent runs.Finally, in contrast to their turbulent runs, we do not report any circumbinary disk formation other than the one in run MU2-VHR where a companion has previously formed.This could be due to the large (> 100 AU) orbital separation in our runs compared to theirs (∼20 AU).
We have investigated the impact of multiplicity on the accretion rate, whose understanding is missing in the star formation context.Indeed, it is not a free parameter but rather an outcome of the self-consistent collapse, even in the studies dedicated to multiple system formation (e.g.Harada et al. 2021).Nonetheless, it has been investigated for binary compact objects (see e.g.Bowen et al. 2018), which are already introduced in the initial conditions.Periodicities in the accretion rate, such as those we investigate in Sec.4.3, have been reported by Bowen et al. (2018), but those are dominated by the periodic inflow from circumbinary disks.Our present work suggests that no circumbinary disk is needed to modulate the accretion rate at orbital-like frequencies.This is in line with the episodes of enhanced accretion at periastron passage reported in Kuruwita & Federrath (2019).
Finally, let us mention the differences with Paper I in terms of multiplicity and disk radius.First, it was shown in Paper I that the multiplicity outcome of the system is sensitive to numerical choices such as the sink properties (number, possibility to merge) and numerical grid as it can seed instabilities leading to fragmentation.Here we have chosen a numerical setup favoring multiple system formation in order to study the influence of magnetic fields on its properties; the approach is somehow similar to Bate (1997), except that here the sinks formed selfconsistently (within the uncertainty related to numerical choices as explicited above).We have checked that the absence of an initial sink particle also leads to a binary system after a few kyr, on which our analysis is performed, and to qualitatively-similar results, which strenghens the results obtained here.This is due to the several mergers presented in Sec.3.3, as mergers are triggered when the multiplicity is higher than 2 in the innermost region of the cloud); for a higher multiplicity to be reached, hierarchical fragmentation is necessary in the present simulations.This study also illustrates how fragmentation, even with an initial central sink particle, leads to symmetry breaking with respect to the axissymmetric initial conditions; allowing the central sink particle to move is key here not to enforce the central symmetry.The present results on the orbital separation also bring additional clues to the debate regarding the link between the early phases and the multiplicity outcome.Indeed, this link between earlyphase, central fragmentation and the ultimate fate of the system is due to the orbital separation increasing, driving stars apart in the form of a multiple system.We show here that such a trend depends on the initial angular momentum and mass distribution of the cloud because the outspiral is driven by gas accretion.Second, we notice that the single disk in Paper I is larger than the disks formed in the binary systems of the HYDRO-runs.In those runs exhibiting a binary system (and in the magnetic runs as well), the disk radius is smaller than the tidal truncation radius, which is ∼0.3 times the orbital separation for equal-mass binaries (Papaloizou & Pringle 1977, Paczynski 1977).Hence, the disk radius is not attributable to tidal truncation by the companion before the secondary fragmentation epoch.This is most likely due to hierarchical fragmentation, in which the specific angular momentum of the fragments is smaller than that of the parent cloud.Hierarchical fragmentation was already mentioned as a possibility to solve the angular momentum problem in star formation (see e.g.Bodenheimer 1978, Zinnecker 1984).After secondary fragmentation however, as the multiplicity becomes larger than 2, the disks are most likely truncated by tidal forces as they extend from one sink to the companion's disk (particularly visible for the HYDRO-HR and HYDRO-VHR runs in Fig. 1).

Comparison with observations
The initial specific angular momentum profile j ∝ R 1.25 (see Sec. 2.2), where we recall that R is the cylindrical radius, is motivated by theoretical purposes (Meyer et al. 2018, Meyer et al. 2019).It may not be representative of typical cores -if such massive pre-stellar cores do exist (see e.g.Zhang et al. 2020 for candidates and Tan et al. 2014, Motte et al. 2018, for recent reviews).For instance, Pineda et al. 2019 measured j ∝ r 1.8 in low-mass dense cores, which is not so different from the measurements of Goodman et al. (1993) who found a specific total angular momentum J/M ∝ r 1.6 (the slope of j and J/M are the same for monotonic profiles, see Pineda et al. 2019).Meanwhile, it could be important in setting the angular momentum versus mass budget and therefore the evolution of the orbital separation of binaries, as studied in Sec.4.2.Nevertheless, this reasoning relies on a high-mass star formation process similar to that of lowmass stars.Large-scale dynamics, such as hierarchical collapse, should be taken into account, but this is beyond the scope of this paper.
We find that, with and without magnetic fields, primary disks form, fragment, and those fragments eventually possess their own disk as well.This hierarchical picture for fragmentation is supported by the recent observations of Suri et al. (2021) in AFGL 2591.Furthermore, the multiplicity of outflows in AFGL 2591 has been used by Gieser et al. (2019) to infer the presence of fragmented cores.Even though outflows are not the main topic of the present paper, let us note that we find magnetic outflows here as in Mignon-Risse et al. (2021a), launched from individual objects in multiple systems (here represented as sink particles) while the collapse is still on-going.Therefore, it confirms magnetic outflows (namely, magnetic tower flows and magnetocentrifugal outflows, as studied in Mignon-Risse et al. 2021a), as relevant candidates for the protostellar outflows observed in multiple protostellar systems.The comparison of polarized emission maps produced in the outflow region with observations is beyond the scope of this paper.
As already reported in Mignon-Risse et al. (2021b), we find a prominent role of spiral arms in the formation of (massive) multiple stellar systems.Interestingly, Tobin et al. (2016) found evidence that the triple protostellar system L1448 IRS3B is still embedded in the spiral arms it may have formed from.This is consistent with the present picture, if low-and high-mass multiple stellar systems share similarities in their accretion process.
Let us note that we do not report runaway stars in those simulations.In runs HYDRO-(of all resolutions) and MU2-VHR, the multiplicity is strictly larger than two, so N−body interactions could possibly eject stars if the computation was integrated for a longer time.Further work is needed to reconcile massive star formation simulations with the high frequency of massive runaway stars reported by observations (see e.g.Dorigo Jones et al. 2020 and references therein).
A large fraction of massive stars are in close binaries (e.g.Sana et al. 2012).While dynamical interactions have been shown to be a relevant process to form close (<10 AU) binaries (Bate et al. 2002), our run MU2-VHR shows that disk fragmentation of magnetized disks could be a possible pathway.Indeed, we find an initial orbital separation of 19 AU (which remained nearly constant for the short integrated time, namely 16 AU after 1.7 kyr) and a common accretion structure.For a possible orbital shrinking to be followed, we would need finer resolution (here, 5 AU is the sink radius).We leave this to further work.
On the opposite, an excess of "twin" (mass ratio between 0.95 and 1) low-mass binaries with orbital separations >1000 AU have been reported in El-Badry et al. (2019) using the Gaia DR2 catalog.This population is unlikely to arise from close binaries widened by positive torques from a circumbinary disk, since no circumbinary disk that large has ever been observed so far, in agreement with the non-magnetized large-scale numerical simulations of Bate (2018).While our results should be scaleddown to low-mass binaries with care, they support gas accretion as a plausible candidate to help forming such a population, after an initial fragmentation phase which could occur within the innermost region of the pre-stellar core.Indeed, we find gas accretion to widen binaries from a few tens AU to 500 AU in the magnetized case with no sign of decrease and with an explicit criterion for further widening explicited in Eq.8.

Conclusions
In this paper, we have investigated the influence of spatial resolution and of magnetic fields on disk fragmentation and on the subsequent formation of massive multiple stellar systems.
We considered eight runs with no restriction on the number of sink particles to be formed, and one sink already present at the center initially.We studied the properties of the multiple stellar systems, varying the resolution and the presence of magnetic fields.At early stages, a density bump-like structure forms after the adiabatic stage and fragments into several sink particles, breaking the axisymmetry via Toomre instability.The sink initially present interacts gravitationally with the other sinks and leaves the center of the cloud.When more than two sinks are present, their gravitational interactions induce radial migration with respect to the center of the cloud.The sinks gaining (losing) angular momentum migrate to larger (smaller) distance from the center, until mergers occur with the other sinks and the multiplicity reduces to 2, in all runs.The orbital separation in binaries increases as they accrete gas possessing angular momentum with respect to the center of mass of the binary, e.g. the center of the cloud.The separation is decreased by a factor of two when magnetic fields are included, as an heritage of their initial orbital separation reduced via magnetic braking compared to the hydrodynamical case.Hence, magnetic fields are promising to explain the formation of close binaries, if there is a common accretion structure.
The increasing orbital separation is driven by accretion of sufficiently high angular momentum gas.Such a trend in the orbital separation is likely to depend on the initial rotation profile.For initial conditions promoting gas-driven outspiral, any fragmentation occuring in the early phases is crucial for the final fate of the system, as briefly presented in Paper I.
A second generation of stars forms from disk fragmentation around the first generation of stars.In more details, fragmentation can be triggered via arm-arm collision or arm-filament collision; in the magnetic case, arm-filament collision is absent.Our results are consistent with a hierarchical picture of fragmentation, in the hydrodynamical case and in the magnetic case (as in, e.g., Park et al. 2022).We have shown that a resolution of 5 AU gives a converged picture in the hydrodynamical case and 2.5 AU in the magnetic case, i.e. about 1/20 of the disk size, in order to resolve structures prone to fragmentation (filaments and spiral arms).Those values depend on the size of the first generation's disks, which are found here to be about twice smaller in the magnetic runs as compared to the hydrodynamical runs, before secondary fragmentation.We only observe secondary fragmentation in the magnetic case for a resolution of 1.25 AU, suggesting that a lack of fragmentation could be, in this particular case, a lack of resolution.Those results will serve as references for future numerical studies.On the opposite, the disk sizes we find here recommend that any resolution above ∼100 AU simply misses disk fragmentation.
Even at low resolution (10 AU here), the convergence regarding the total sink mass (15%) and the orbital separation of binary systems is correct.In the accretion rates of binaries we find hints of modulations at the orbital frequency.This suggests that their tidal forces help the accretion onto the other object.If confirmed by further work, it means that fragmentation does not limit the mass of stars but rather enhances accretion bursts and possibly increases the efficiency of gas conversion into stars in multiple systems.This is of particular interest, since massive stars are more often in multiple systems than low-mass stars (e.g.Chini et al. 2012).This moderates the argument of the fragmentation-induced starvation scenario (Peters et al. 2010), in which fragmentation reduces the maximal stellar mass by sub-dividing the mass reservoir into several lower-mass stars.
Overall, this work suggests that the influence of magnetic fields in the formation of multiple masse stars -without even focusing on the outflows, which are now understood as mainly magnetic -is important in several aspects.We find that magnetic fields regulate the disk size, which appears to be, otherwise, only regulated through disk fragmentation producing a companion, and tidal truncation by the same companion.Among the possible fragmentation modes (arm-arm collision, arm-filament collision), we find the arm-filament collision mode to be suppressed because the filaments are more diffuse in the presence of magnetic fields.Finally, magnetic braking reduces the orbital separation of binaries by a factor of two because i) it reduces their initial separation and ii) it has no impact on the further evolution of the orbital separation, driven through accretion.In the absence of a common accretion structure, magnetic effects are insufficient to form close binaries.For their strong impact on the disk growth and fragmentation and on the initial orbital separation, we conclude that magnetic fields are important in the formation of multiple massive stars.

Fig. 1 .
Fig. 1.Column density computed over 200 AU along the x−axis (the rotation axis) in the hydro runs (top) and in the magnetic runs (bottom) at low resolution (left panels), mid-resolution (center-left), high resolution (center-right), and very high resolution (right panels), displayed in the x = 0 plane, at the end of the runs.White dots indicate sink particles positions.Gas velocity vectors are overplotted.

Fig. 2 .
Fig.2.Minimal Jeans length as a function of the maximal AMR level l max (a higher AMR level means a finer resolution), computed at t = 10 kyr: each point refers to a simulation of Table1.In the magnetic case, we compute the thermal (L J , see Eq. 3) and magnetic (L J,mag ) Jeans length.An additional AMR level gives a finer resolution by a factor of two.

Fig. 3 .
Fig. 3. Number of sink particles as a function of time in the hydro runs (left panel) and in the magnetic runs (right panel).

Fig. 4 .
Fig.4.Trajectory of sinks #1 and #2 in run HYDRO-HR in the (y, z)plane (perpendicular to the initial rotation axis).The sinks location at t = 5 kyr, t = 10 kyr, and at the time of secondary fragmentation, is indicated by the dots, the triangles, and the empty squares, respectively.

Fig. 6 .
Fig. 6.Sum of all sink masses in the HYDRO runs (left panel) and the MU2 runs (right panel) as a function of time.On the left panel, the filled circles indicate the secondary formation epoch.The blue region shows the value obtained in the highest resolution runs ±20%.

Fig. 7 .
Fig. 7. Map of B/B 0 , where B 0 is the initial, uniform magnetic field strength, and B φ /B in the orbital plane in run MU2-HR at t = 10 kyr.Sink particles are denoted by white circles.

Fig. 8 .
Fig. 8. Gravitational torque rg φ acting on binaries in runs HYDRO-HR and MU2-HR as a function of time.

Fig. 9 .
Fig. 9. Accretion rate in the HYDRO runs (left panel) and the MU2 runs (right panel) as a function of time.On the left panel, the filled circles indicate the secondary formation epoch.

Fig.
Fig.10.Fourier power spectrum of the accretion rate in run MU2-HR after 10 kyr.The x−axis is the Fourier transform frequency normalized by the approximated orbital frequency Ω = 1/2.5 kyr −1 .Data (instantaneous accretion rates) have been linearly interpolated in order to have a uniformly-spaced sample to perform the fast Fourier transform.The gray band indicates the 0.7 − 2Ω orb interval.The low-frequency signal does not point toward a physical periodicity but to the sample time which is only 6 kyr long.More interestingly, we find no significant peak at frequencies higher than 2Ω orb .This suggests that the accretion rate is modulated by frequencies of the order of the orbital frequency Ω orb .

Fig. 11 .
Fig. 11.Disk radius as a function of time, computed from the surface density with a threshold value N = 10 25 cm −2 (black curve), N = 5 × 10 25 cm −2 (gray curve) and the plasma beta (green curve) in run MU2-VHR, and with the plasma beta in run MU2-HR (green dashed curve).

Table 1 .
Finest spatial resolution and magnetization level of the different runs.The cost of each run is given in Appendix A.

Table A .
1. Computational cost (in CPUkhr) and CO 2,e footprint estimate (in kg) of the simulations presented in Table1.