Angular Momentum Transport in Binary Star Formation: The Enhancement of Magnetorotational Instability and Role of Outflows

The formation of binary stars is highly influenced by magnetic fields, which play a crucial role in transporting angular momentum. We conducted 3D numerical simulations of binary star accretion via a circumbinary disk, taking into account a magnetic field perpendicular to the disk and an infalling envelope. Our simulations reproduce the following phenomena: (1) the spiral arms associated with circumstellar disks; (2) the turbulence in the circumbinary disk, induced by magnetorotational instability (MRI); (3) a fast outflow launched from each circumstellar disk; and (4) a slow outflow from the circumbinary disk. The binary models exhibit a higher α-parameter than the corresponding single-star models, indicating that the binary stars enhance the MRI turbulence. Moreover, an infalling envelope also enhances the turbulence, leading to a high α-parameter. While the spiral arms promote radial flow, causing the transfer of mass and angular momentum within the circumbinary disk, the MRI turbulence and outflows are the main drivers of angular momentum transfer to reduce the specific angular momentum of the system.


INTRODUCTION
About half of solar-type stars are formed as members of multiple systems, and the multiplicity increases for higher stellar masses (Duchêne & Kraus 2013).Therefore, binary star formation is an important process, and several formation scenarios have been proposed so far (Reipurth et al. 2014;Offner et al. 2022).Recent highresolution observations with the Atacama Large Millimeter/Submillimeter Array (ALMA) have revealed the early phases of low-mass binary and multiple star formation (Hioki et al. 2007;Fukagawa et al. 2013;Dutrey et al. 2014;Takakuwa et al. 2014;Alves et al. 2019;Takakuwa et al. 2020).The binary protostars surrounded by circumbinary disks imply a contribution of formation mechanism based on the disk fragmentation scenario (Matsumoto & Hanawa 2003;Kratter et al. 2010).
It is known that magnetic fields play a crucial role in star formation.Magnetic fields give rise to vari-Corresponding author: Tomoaki Matsumoto matsu@hosei.ac.jp ous phenomena in the star formation process, such as launching outflows, exerting magnetic braking on disks, (e.g., Machida et al. 2011;Tsukamoto et al. 2022) and generating turbulence via magneto-rotational instability (MRI) (e.g., Balbus & Hawley 1991;Bai & Stone 2013).However, researches on the role of magnetic fields have primarily focused on single star formation rather than binary star formation.These magnetic processes induce angular momentum transport within the system, and in the context of binary star formation, the magnetic fields are expected to influence fundamental binary parameters such as binary separation.
Binary accretion from a circumbinary disk has been intensively investigated in many studies, not only for binary black holes but also for binary star formation (see recent review Lai & Muñoz 2022).These studies usually assume α viscosity model (Shakura & Sunyaev 1973), where angular momentum transport occurs in the circumbinary disk.The assumed viscosity affects the orbital evolution of binary stars (e.g., Dittmann & Ryan 2022).The origin of the α viscosity is thought to be MRI turbulence.Once a magnetic field is assumed in a disk model, it causes not only MRI (Noble et al. 2021) but also disk winds (Suzuki & Inutsuka 2014) and outflows (Machida et al. 2011), which lead to the redistribution of angular momentum.
Several researchers have considered magnetic fields in binary accretion models; Noble et al. (2012); Shi et al. (2012); Shi & Krolik (2015); Lopez Armengol et al. (2021); Noble et al. (2021) incorporated magnetic fields in their simulations, but the initial magnetic field configuration was confined within the circumbinary disk, where the closed magnetic field lines followed the isodensity surfaces of the disk.Bowen et al. (2018); Avara et al. (2023) also took into account magnetic fields, using the initial conditions derived from snapshots of the simulations by Noble et al. (2012).Although the magnetic fields change during the evolution of the circumbinary disks, the initial configuration of the magnetic fields may affect the disk's evolution.In the case of binary star formation with typical field strength, the natal molecular cloud core collapses along the magnetic fields, and a disk perpendicular to the open magnetic field naturally forms (e.g., Matsumoto et al. 2017;Tsukamoto et al. 2018).Moreover, the perpendicular component (vertical component) of the magnetic field is crucial for launching outflows (Machida et al. 2011;Gerrard et al. 2019), which play an important role in angular momentum transport.
Typical binary accretion models employ an infinitely extended circumbinary disk in a steady state as an initial condition (e.g., Moody et al. 2019).However, in the context of binary star formation, recent ALMA observations reveal that circumbinary disks actually have a finite size (e.g., Takakuwa et al. 2020).Furthermore, in the early stages of star formation, a young star is embedded in a cloud core and accretes gas from an infalling envelope (e.g., Hayashi et al. 1993).Therefore, a circumbinary disk model with a finite size, which accretes gas from an infalling envelope, is more suitable for the context of binary star formation, as adopted by Bate & Bonnell (1997).The assumed infalling envelope has the potential to influence the evolution of the binary system (Bate 2000).
In this paper, we examine the accretion of binary stars from a finite-sized circumbinary disk, considering a magnetic field that perpendicularly threads the disk.Additionally, we incorporate an infalling envelope, which significantly contributes to the mass and angular momentum supply to the system, as well as to angular momentum transport through magnetic braking.Despite its importance, this factor has often been overlooked in many previous studies.Consequently, our models not only reproduce the MRI in the circumbinary disk but also simulate the outflows from binary stars and accretion onto the circumbinary disk from an infalling envelope.By employing this model, we conduct a quanti-tative investigation into angular momentum transport in a binary system.The present model reproduces two types of disks: a circumbinary disk and circumstellar disks.Our focus is primarily on the circumbinary disk, as the circumstellar disks are likely influenced by magnetic diffusion processes, such as Ohmic dissipation, ambipolar diffusion, and the Hall effect.We assume the ideal-MHD, and these processes are not included in our model for simplicity.The outflows are also investigated as a case of the ideal-MHD limit.
This paper is organized as follows.In Sections 2 and 3, the model and methods are shown.The results of the simulations are presented in Section 4, and they are discussed in Section 5. Finally, the conclusions of this paper are given in Section 6.The details of the analyses are shown in Appendix.

MODELS
The numerical models presented in this paper are based on Matsumoto et al. (2019b), which is extended to include magnetic fields.The primary and secondary stars have constant masses of M 1 and M 2 , respectively.It is assumed that the stars rotate around the origin in fixed circular orbits at an angular velocity of Ω ⋆ = (GM tot /a 3 b ) 1/2 , where G is the gravitational constant, M tot (= M 1 + M 2 ) is the total mass of the stars, and a b is binary separation.The rotation period of the binary stars is T ⋆ = 2π(a 3 b /GM tot ) 1/2 .A cylindrical computational domain, as shown in Figure 1, is used in the simulations.The radius of the cylinder is set at L (= 12a b ), and the height of the cylinder is [−L, L].The gas is injected at the boundary surfaces, i.e., the side, top, and bottom of the cylinder.The gas injection mimics a rigidly rotating, spherical infalling envelope around binary protostars.The gas has a density distribution of ρ(r) = ρ 0 (r/L) −1.5 at the boundary surfaces, where r denotes the radius in spherical coordinates.The power index of −1.5 is for the infalling envelope of the mass accretion phase.The injected gas has a constant angular velocity of Ω inf and a radial velocity of v r = (2GM tot /r − (Ω inf R) 2 ) 1/2 at the boundary surfaces, assuming freefall from a distance of infinity.Here, R denotes the cylindrical radius.The injected gas has a range of specific angular momenta 0 ≤ j ≤ j inf , where j inf = Ω inf L 2 is the maximum specific angular momentum of the infalling gas.
The barotropic equation of state is assumed, in which the gas pressure is given by p = c 2 s ρ[1 + (ρ/ρ cr ) 2/5 ].This equation of state is an approximation obtained from radiation transfer calculations (Masunaga et al. 1998) and is commonly used in numerical simulations for protostellar collapse (e.g., Matsumoto et al. 2017).The gas is approximated as isothermal if the density is less than ρ cr and polytropic with an adiabatic index γ = 7/5 if the density is larger than ρ cr .In this paper, the critical density of ρ cr is set at 10 4 ρ 0 .Consequently, the gas is isothermal in an infalling envelope and a circumbinary disk, while it is polytropic in circumstellar disks.According to Masunaga et al. (1998), the critical density is given by ρ cr ∼ 10 −13 g cm −3 (with the corresponding number density n cr = 10 10−11 cm −3 ).The envelope gas is therefore assumed to have a density of ρ 0 ∼ 10 −17 g cm −3 (with the corresponding number density n 0 = 10 6−7 cm −3 ), which is typical for a scale of 1000 au around protostars (Onishi et al. 2002;Tokuda et al. 2016Tokuda et al. , 2020)).The barotropic equation of state considers the energy balance between compression heating and radiation cooling during protostellar collapse.In the isothermal region, this equation assumes that the timescale for radiation cooling is significantly shorter than that for heating processes, such as shock heating and compression heating due to gravity.While the barotropic equation of state accurately replicates the gas temperature in the central region of a collapsing gas cloud, it underestimates the temperature by a factor of 2 − 3 on a scale of approximately 10 au when compared with results from radiation hydrodynamics (Tomida et al. 2010).
The self-gravity of the gas is ignored, implying that the models here can be applied to a binary system where the total mass of the gas is much less than those of the stars.Models that take into account the gas-to-star interaction will be reported in a future paper.
Before starting an MHD simulation, a hydrodynamical calculation is performed for 10 rotation periods of the binary stars without magnetic fields to construct the initial condition.This hydrodynamical calculation provides a quasi-steady state of the circumbinary and circumstellar disks, as shown in Matsumoto et al. (2019b), and the magnetic field is imposed on this state.The simulation clock is set at t = 0 when the MHD calculation begins.The initial magnetic field is uniform and parallel to the z-direction, and its strength is denoted by B z,0 .We assume ideal-MHD in the evolution, and magnetic diffusion processes such as Ohmic dissipation and ambipolar diffusion are not considered.We also assume an inviscid gas, excluding any artificial viscosity and αviscosity, both in constructing the initial condition and during the evolution after the magnetic field imposition.
We adopt the units of a b = 1, GM tot = 1, and ρ 0 = 1.The model parameters are the mass ratio of the binary stars q = M 2 /M 1 , the isothermal gas sound speed c s , the maximum specific angular momentum of the infalling gas j inf , and the initial magnetic field strength B z,0 .In this paper, we examine models with q = 0.2, c s = 0.1, j inf = 1.2, and varying B z,0 to investigate the effects of the magnetic field (Table 1).
The sound speed of c s = 0.1 corresponds to 0.1(GM tot /a b ) 1/2 = 0.30 km s −1 when assuming M tot = 1M ⊙ and a b = 100 au.The disk thickness is related to the sound speed, given by H/R = c s (R/GM tot ) 1/2 .In the case of c s = 0.1, the disk thickness is H/R = 0.1(R/a b ) 1/2 .Note that the simulated circumbinary disk has a thickness H ∼ a b (Figure 3), which is much thicker than that estimation, because of disturbance by the spiral arms.The circumstellar disk around the primary star exhibits a density of ρ ∼ 10 5 ρ 0 at its thickest radius (approximately 0.3a b ).Given the barotropic equation of state assumed here, the corresponding sound speed is 0.21(GM tot /a b ) 1/2 .This leads to an estimated disk thickness of H/R = 0.21[(1 + q)R/a b ] 1/2 ∼ 0.13 when q = 0.2, which is consistent with, but smaller than, the thickness of the simulated disk, approximately 1/3.Note that this paper mainly focus on the circumbinary disk rather than the circumstellar disks.
The initial magnetic field of B z,0 = 0.1 corresponds to 94.2 µGauss, under the same assumptions and ρ 0 = 10 −17 g cm 3 , which corresponds to a number density of n 0 = 2.61 × 10 6 cm −3 for a mean molecular weight of 2.3.The Alfvén speed in the infalling envelope is therefore v A = B z,0 /(4πρ 0 ) 1/2 = 0.028(GM tot /a b ) 1/2 , corresponding to v A = 0.084 km s −1 .The most unstable mode of MRI is given by k z v A /Ω ≃ 1 for a wave number k z (Balbus & Hawley 1991).In the case of B 0,z = 0.1, the corresponding wavelength is λ max = 0.1π 1/2 (ρ 0 /ρ) 1/2 (R 3 /a b ) 1/2 when assuming the Keplerian rotation.At the initial stage, the density in the circumbinary disk is ρ ∼ 10 2 ρ 0 − 10 3 ρ 0 , corresponding to λ max ∼ 0.01a b − 0.05a b at R = 2a b .The wavelength is less than the disk thickness (H ∼ a b ).Therefore, the circumbinary disk is unstable against the MRI.A typical cell width covering the circumbinary disk is ∆x = 0.023a b − 0.0058a b for ℓ = 2 − 4 (see section 3), indicating that the most unstable MRI mode is marginally resolved at the initial stage.In Appendix A, we investigate the numerical resolution for resolving the MRI through the evolution.
To investigate the effects of the infalling envelope, we also computed models without the infalling envelope, where gas injection from the boundary surfaces is halted at t = 0.These models are denoted as "N" in the "Infall" column of Table 1.Additionally, we examined models of a single star (q = 0) to draw comparisons with the binary.

METHODS
The numerical simulations were performed using the adaptive mesh refinement code, SFUMATO (Matsumoto 2007), which employs fixed mesh refinement (FMR) with a static hierarchical grid configuration shown in Figure 1.The computation domain is covered by 16 3 blocks for a grid level of ℓ = 0 (the base grid), and each block has 16 3 cubic cells.The base grid has therefore a resolution of 256 3 cells, and the cell width is ∆x max = 0.0938.The maximum grid level is set at ℓ max = 4, and the minimum cell width is ∆x min = 0.00586.
The simulations use a MHD equation solver with a third-order accuracy in space using the MUSCL method and second-order accuracy in time with the predictor-  (Matsumoto et al. 2019a) is adopted to calculate the numerical flux.This scheme is based on the HLLD Riemann solver (Miyoshi & Kusano 2005) and modified to solve the MHD equation with Boris correction (Gombosi et al. 2002).The modified scheme suppresses the Alfvén speed below the reduced speed of light, which should be several times larger than the maximum velocity in the computation domain, |v| max , for a stable calculation (Matsumoto et al. 2019a).
The reduced speed of light c red is varied according to the gas velocity as c red = max(4|v| max , 50) to avoid a very short timestep when the density becomes very low, allowing for long-term evolution of the simulations.
The sink particles are used to represent the binary stars in the simulations.These sink particles accrete gas within a sphere of radius r sink and density higher than ρ sink (Matsumoto et al. 2015).We set r sink = 4∆x min = 0.0234a b for the sink radius, which is much smaller than the binary star separation, so the sink particles have little hydrodynamical impact on the circumstellar disk.Despite the r sink assumed here being only four times larger than the cell size, influences from finite cell size, such as quadrupole structures, were not observed around the sink regions in the models.The density threshold for accretion onto the sink particles is set to ρ sink = 10 6 ρ 0 , which is a typical density at the center of the circumstellar disks.Although the sink particles accrete mass and angular momentum, their masses and orbits remain fixed during the evolution.
Even though the sink particles accrete gas within the sink regions, gas with a density below the threshold ρ sink remains.If this gas rotates and is coupled with the magnetic fields, outflows would be artificially launched from the sink regions.To avoid this issue, we impose Ohmic dissipation only in the sink regions with a high resistivity, η = v K r sink , where v K = (GM sink /r sink ) 1/2 is the Keplerian velocity for the sink particle with a mass of M sink .With this measure, the artificial launching of the outflows was not observed inside the sink regions.
During the accretion process of a sink particle, we extract the gas with a density exceeding ρ sink within the sink region, while the magnetic field in that area is retained.This process results in an accumulation of magnetic flux in the sink region when the magnetic field lines have a open configuration, as it is in our present models.A more realistic treatment of magnetic fields related to a sink particle will be studied in a future paper.
The circumstellar and circumbinary disks in the simulations launch outflows that eventually reach the top and bottom surfaces of the cylindrical boundary.At this point, the boundary condition is switched from the gas-injection to the out-going boundary condition for outflows by applying the zero-gradient boundary condition.While the out-going boundary condition allows the outflows to flow out of the computational domain, some waves are reflected by the boundary and return to the computational domain.These reflections have little impact on the evolution of the disks, as the outflows reach the boundaries only in the very last stages (see the lower right panel of 2 for the fiducial model).

Overall evolution of the fiducial model
Evolution of the fiducial model with B z,0 = 0.1 is shown in Figures 2 and 3.At the initial stage (t = 0), each star is surrounded by a circumstellar disk, and the binary stars are surrounded by a circumbinary disk.In the circumbinary disk, two spiral arms are excited because of the gravitational torque of the binary stars.The initial condition is constructed by a pure hydrodynamical calculation.Unlike typical binary accretion models (e.g., Moody et al. 2019), this model has no clear cavity or gap at the initial stage.This is because the gas accreting in the vertical direction has small angular momenta to fill the cavity.At t = 0, a uniform magnetic field is imposed, and the MHD calculation starts.
The two spiral arms continue to exist in the circumbinary disk during the MHD calculation.The circumbinary disk gradually extends due to the angular momentum redistribution by the magnetic field.The circumstellar disk around the primary star remains while the disk around the secondary star becomes obscured due to magnetic braking.This phenomenon can be attributed to differences in magnetic field fluxes within the circumstellar disks.The secondary star experiences a higher accretion rate compared to the primary star (see Figure 21).As a result, the secondary circumstellar disk accumulates more magnetic flux than the primary circumstellar disk, leading to more significant angular momentum transport due to stronger magnetic braking in the secondary circumstellar disk.
The outflows begin to be launched from the two circumstellar disks at t ∼ 5T ⋆ (Figure 3).They propagate in the vertical direction and reach the top and bottom boundary surfaces of the computational domain at t = 10T ⋆ (Figure 2).At the launch points at t = 10T ⋆ , the outflow from the primary star has a higher velocity (v z = ±(5 − 8)) than that from the secondary star (v z = ±(1 − 3)) (Figure 4).As showin in the upper left panel of Figure 5, the outflows shown by the red and blue lobes are also twisted due to the orbital motion of the binary stars.As shown in the lower right panel of Figure 5, the magnetic field lines are tightly twisted along the outflow directions, suggesting that magnetic pressure plays a significant role in accelerating the outflows (Tomisaka 2002;Machida et al. 2008;Tomida et al. 2010).The higher angular velocity at the inner parts of the circumstellar disks results in outflows being preferentially launched from these regions, while in the sink region, where the gas is decoupled due to high resistivity, outflows are not initiated (see Section 3).Although centrifugal winds, driven by less twisted magnetic field lines, represent another outflow mechanism (Blandford & Payne 1982;Pudritz & Norman 1986), they are not relevant to this model.
In addition to the fast two outflows stated above, a slow outflow is launched from the circumbinary disk.The slow outflow has velocity of v z = ±(0.2− 0.5) as shown in Figure 4.The slow outflow has a cocoon-like shape with highly twisted magnetic fields as shown in the upper right panel of Figure 5.The rotation of the circumbinary disk is responsible for the twisted magnetic field, which implies that the magnetic pressure accelerates the outflow.
Figure 6 shows the evolution of the plasma β, the ratio of thermal pressure and magnetic pressure denoted by β p .For the fiducial model, the infalling gas has β p ≃ 8πc 2 s ρ 0 /B 2 z,0 = 8π due to c s = 0.1 and B z,0 = 0.1.At the initial stage (t = 0), the circumbinary disk has a typical value of β p ∼ 10 3 −10 4 because the density in the disk reaches (10 2 − 10 3 )ρ 0 .As time passes, β p decreases because the magnetic field strength increases.The increase in the magnetic field is caused by mainly three processes.The first process is accretion due to the infalling envelope, which brings not only gas but also magnetic flux into the computational domain.This accretion process considerably increases the magnetic flux.The second process is amplification of the magnetic field by rotation of disks.The rotation of the circumbinary disk and circumstellar disks amplify the toroidal magnetic field, which are associated with the slow and fast outflows, respectively, as shown in Figure 5.In the fast and slow outflow regions, the toroidal magnetic field dominates over the poloidal fields.Observations by Ching et al. (2016); Lee et al. (2018) suggest similar magnetic field structures associated with outflows, where the toroidal magnetic fields are observed wrapping around the outflows.The third process is the amplification of the magnetic field by turbulence in the circumbinary disks.The MRI excites turbulent flows in the circumbinary disk, as shown in the fluctuated density distribution (Figure 3). Figure 7 shows magnetic field lines entangled within the circumbinary disk, indicating the development of MRI.Furthermore, the magnetic field lines have an overall spiral structure due to the redistribution of angular momentum.The spiral magnetic fields contribute to the spoke-like density structures observed in the face-on views of Figures 2 and 3.These spoke-like density structures are likely caused by channel flows, which are commonly observed in MRI simulations (e.g., Sano & Inutsuka 2001).Layered accretion also occurs on the surface of the circumbinary disk, dragging the magnetic field lines inward in the cylindrical radial direction, as shown in the upper-right panel of Figure 5. Due to this process of angular momentum redistribution, the gas near the mid-plane slowly moves outward in the cylindrical radial direction, leading to the expansion of the outer part of the disk.Similar processes are seen in MHD simulations of an accretion disk around a single star with vertical magnetic fields (Suzuki & Inutsuka 2014).

Turbulence in the circumbinary disks
The turbulence is generated in the circumbinary disk.In order to investigate turbulent level, we measure the αparameters for the Reynolds and the Maxwell stresses of the turbulent components, α R,turb and α M,turb , according to Appendix B.
Figure 8 shows α R,turb and α M,turb for the fiducial model.As shown in lower panels of Figure 8, both the α-parameters exhibit almost the same level, having large values of 0.5 − 0.7 along the spiral arms.This indicates that the spiral arms enhance a turbulent level in the circumbinary disk.Even in the inter-arm regions, the α-parameters exhibit a moderate level of ∼ 0.3, which is larger than that shown in the corresponding single star model (see section 4.4, and panels (h) of Figures 15  and 16).
The α-parameters are measured in the rotating frame in which the binary stars are rest.Even in the rotating frame, the spiral arms are not perfectly steady and fluctuate with time because of inviscid flow.The spiral arm of the secondary star fluctuates more.The high αparameters along the spiral arms is partially attributed to such fluctuation of them.However, they can also disturb gas flows and enhance a turbulent level in the circumbinary disk.
A high value of α R,turb is observed at the outer edge of the circumbinary disk (upper left panel of Figure 8).This is not due to turbulence, but instead can be attributed to the non-circular shape of the disk.This shape causes temporal changes in density at the outer edge, leading to the high value of α R,turb .
Angular momentum transport in the circumbinary disk is not only due to MRI turbulence but also due to mean flows and mean magnetic fields.Figure 9 shows The mean flow of ρ v R v φ has a large amplitude that follows the shape of the spiral arms, as shown in Figure 9 (left).Along the spiral arms, the peak value is α R,mean ∼ 2, while it has large negative values α R,mean ∼ −10 in inter-arm regions.This suggests that the spiral arms efficiently transport angular momentum outward in the spiral arms and inward in the inter-arm regions.The gravitational torque from the binary stars generates the spiral arms and is responsible for the associated angular momentum transport.This transport is also seen in hydrodynamical simulations by Matsumoto et al. (2019b) and MHD simulations by Shi et al. (2012), where gas exhibits expansion motion along the spiral arms and infall motion in the inter-arm regions.These motions have been observed by ALMA for protobinary systems L1551 NE (Takakuwa et al. 2014(Takakuwa et al. , 2017) ) and L1551 IRS 5 (Takakuwa et al. 2020).
The angular momentum transport by the mean magnetic field is shown in Figure 9 (right).The ring-shaped region with a radius of ∼ 3 exhibits α M,mean ∼ 0.5.In this region, the angular momentum is transferred outward by the mean magnetic field of B R B φ , whose magnetic field configuration can be seen in Figure 7.The magnetic field is stretched in the radial direction in the circumbinary disk.
Figure 10 shows the radial distribution of the αparameters for comparing the various components quantitatively.In the region of 1 ≲ R ≲ 4, the α-parameters of the turbulent components (ρv ′ R v ′ φ and B ′ R B ′ φ ) and the mean magnetic field component (B R B φ ) exhibit ∼ 0.5.The large value of the velocity turbulent component (ρv ′ R v ′ φ ) at R ∼ 6 is due to the outer edge of the circumbinary disk, as noted in the upper-left panel of Figure 8.
The mean flow component (ρ v R v φ ) exhibits a large negative value in the circumbinary disk, as expected from the left panel of Figure 9.It overwhelms the other α components and is attributed to gas accretion through  the circumbinary disk to the central circumstellar disks or binary stars.

Gas structures: dependence on binary parameters
Figures 11 and 12 show the surface density distributions in face-on view, comparing the structures of the circumbinary gas among all the models.We compare the models at almost the same stage of t ∼ 10T ⋆ , although the weak field models and non-magnetized models are calculated for a longer time.
When comparing Figure 11a-e, we see that a model with stronger magnetic fields has an extended circumbi-nary disk, which is attributed to angular momentum redistribution by MRI.Because of MRI, more turbulent density structure is seen in a model with a stronger magnetic field, as shown in Figure 12a-e.Furthermore, the cavity of the circumbinary disk and spiral arms are obscure for the strong magnetic field models.
A model with a stronger magnetic field also shows stronger outflows.When we observe outflows at the same stage of t ∼ 10T ⋆ , a model with a stronger magnetic field launches outflows farther, as shown in panels (a)-(e) of Figures 13 and 14.Note that the very weak magnetic field model (Figure 13d) exhibits short  and 11f; Figures 11e and 11g), we observe that the models without an infalling envelope have extended circumbinary disks because there is no ram pressure from the infalling gas.
The models with infalling envelopes exhibit more disturbed density distributions in the circumbinary disks compared to those without, as seen when comparing Figures 12b and 12f, which have the same initial magnetic field strength of B z,0 = 0.1.This difference is mainly due to two factors.The first is difference in magnetic flux between the models with and without an infalling envelope.The infalling envelope not only brings mass but also magnetic flux through the boundary surfaces into the computational domain.As a result, the model with an infalling envelope has 2.5 times more magnetic flux than the model without at the stages shown in the figures, and therefore exhibits a more developed MRI.The increase in the magnetic flux is disscussed in section 5.4.
The second factor is disturbance by the infalling gas onto the circumbinary disk.A fraction of the accretion energy is converted into heat at the accretion shock.However, because a barotropic equation of state is assumed here, this heat is radiated away, keeping the gas temperature constant.A certain part of the kinetic energy of the infalling gas is converted to the kinetic energy of the turbulent motion (Klessen & Hennebelle 2010).This effect is also reported in previous simulations by Matsumoto et al. (2019b).The accretion-  driven turbulence is also observed when comparing nonmagnetized models in Figures 12e and 12g.We also note that the non-magnetized model without an infalling envelope (Figure 12g) exhibits a cavity with the highest contrast because it has the lowest turbulent level among the binary models examined here, and the absence of accreting gas that would otherwise fill the cavity.
Panels (h)-(k) of Figures 11 and 12 show the models of single stars.A comparison between the binary star model and the single star model with the same initial magnetic field (e.g., Figures 12b and 12h) shows that the binary model has a more turbulent density distribution than the single star model, suggesting that the orbital motion of the binary stars disturbs the circumbinary disk.

Alpha parameters: dependence on binary parameters
Figures 15 and 16 compare the α-parameters of turbulent components of Reynolds stress and Maxwell stress, respectively, among the models examined here.
Comparison of panels (a)-(e) in Figures 15 and 16 reveals the influence of magnetic fields on the level of turbulence in the circumbinary disk; models with stronger initial magnetic fields have stronger turbulence except for the vicinity of binary stars and the disk edges.
The effect of the infalling envelope in the models with B z,0 = 0.1 can be observed by comparing panels (b) and (f) in Figures 15 and 16.Similarly, panels (e) and (g) show the effect of the envelope in the models with B z,0 = 0.In both cases, the models with infalling envelopes exhibit higher α-parameters within the cir-   cumbinary disks, which is consistent with the differences in density distributions discussed in section 4.3.
The panels (b) and (h), (e) and (i), (f) and (j), and (g) and (k) demonstrate the impact of binarity on the turbulence levels in disks for various models with/without magnetization and infalling envelopes.These figures indicate that the binary models tend to have a higher level of turbulence than the single star models in their disks.For example, the single star model in panel (h) exhibits a very small turbulent level of α R,turb ∼ 0.01 and α M,turb ∼ 0.05 in the central region (R ≲ 1), and a moderate level of α R,turb ∼ α M,turb ∼ 0.3 at R ∼ 3 (refer to Figure 17), but it is lower than that of the corresponding binary model, which is shown in Figure 10.

Angular momentum transport
In the previous sections, we explored the αparameters, which are indicators of angular momentum transport in the radial or horizontal direction.In this section, we investigate angular momentum fluxes to examine transport in both the horizontal and vertical directions.
Figure 18 displays the angular momentum fluxes in both the R and z-directions for the fiducial model.Negative flux values for F R and F z represent inflow of angular momentum, while positive values represent outflow.A detailed description of the methods for calculating these fluxes is shown in Appendix C.
We found that, in the circumbinary disk region (1 ≲ R ≲ 5), the T Rφ component (blue solid line) is dominant, corresponding to the gas flow associated with spiral arms.The turbulent components of T Rφ and M Rφ (dashed lines) exhibit significant positive values, indicating an outgoing transfer of angular momentum due to MRI turbulence.The turbulent component of M Rφ (orange dashed line) accounts for half of the total positive angular momentum flux of M Rφ (orange solid line) in this region.The other half is attributable to the coherent component of the magnetic field, resulting from the winding of the magnetic field (see Figure 7).The turbulent component of T Rφ peaks at the outer edge of the circumbinary disk (R ∼ 6) due to the the non-circular outline of the disk (see Figure 2).
In the infalling envelope (R ≳ 5), the T Rφ component (blue solid line) dominates the angular momentum flux in the radial direction.Similarly, in the vertical direction, the T zφ component (blue solid line) is the dominant factor across the range of the model.Notably, these values are negative, indicating that the infalling envelope transports both mass and angular momentum towards the central region in both the radial and vertical directions.
The angular momentum flux due to the outflows (green line) is shown in the lower panel of Figure 18.The outflows are responsible for a large portion of the outward flux of angular momentum in the z direction for z ≲ 1, but it does not exceed the angular momentum inflow due to the infalling envelope.Additionally, we note that the outward angular momentum transport in the upper region of the circumbinary disk (z ≳ 1) is carried out by magnetic braking (M zφ component; orange line) rather than outflows (c.f., Marchand et al. 2020;Lee et al. 2021).
The angular momentum flux represents the angular momentum flowing into or out of the region of interest.On the other hand, to consider spin-up or down of the region of interest, it is necessary to consider specific angular momentum, which is approximately angular momentum per mass.We estimate the timescale of specific angular momentum change (τ sam ), following the method provided in Appendix D. In the following, the rate of change in specific angular momentum (τ −1 sam ) is shown, of which definition is given by equation (D34).A positive τ −1 sam represents increase in the specific angular momentum.
Figure 19 shows τ −1 sam for the fiducial model.In the circumbinary disk region of 1 ≲ R ≲ 5, there are a negative contribution due to magnetic field (green line) and a positive contribution due to gas flow (blue line) in τ −1 sam , indicating that the magnetic field, including MRI turbulence, reduces specific angular momentum, and the gas flow associated with the spiral arms increases the specific angular momentum.As a result of these competing contributions, the specific angular momentum decreases in outer region the circumbinary disk and increases in the inner region the disk (dotted line).
In the outer region of R, z ≳ 7, a negative τ −1 sam in total (dotted line) is exhibited, which is responsible for the angular momentum transport in the z-direction, caused mostly by gas flow due to the outflow and magnetic braking (orange and red lines).
Figure 20 shows τ −1 sam for representative models.The binary model with magnetic field and infalling envelope (the fiducial model; blue line) shows high negative values of τ −1 sam at R, z ∼ 4. The models with magnetic fields (blue and orange lines) tend to show higher negative values of τ −1 sam than those without (green and red lines).This is because the magnetic effect, including MRI, plays an important role in the angular momentum transport.The single star models (lower panel of Figure 20) exhibit smaller negative values of τ −1 sam than the binary models, indicating that the orbital motion of the binary stars enhances MRI turbulence, as mentioned above.
. Same as Figure 10 but for the single star model (q = 0) with B0,z = 0.1 and an infalling envelope.
The high negative value of τ −1 sam observed in the magnetic binary models indicates a reduction in specific angular momentum.This suggests that the binary separation could decrease due to the transport of angular momentum facilitated by the magnetic field.A more accurate estimation of the changes in binary separation will be the focus of future investigations.

Mass accretion
The mass accretion onto binary stars has been a topic of debate (Bate & Bonnell 1997;Ochi et al. 2005;Hanawa et al. 2010;Young et al. 2015;Young & Clarke 2015;Satsuka et al. 2017;Matsumoto et al. 2019b), and in general, binary stars accrete gas in a way that increases the mass ratio, leading to the evolution of the system towards a twin binary system.To describe this change in the mass ratio, the Γ parameter is defined as where a positive value of Γ indicates an increase in the mass ratio (Bate & Bonnell 1997;Young et al. 2015;Matsumoto et al. 2019b).Figure 21 depicts the mass accretion rates onto the primary and secondary stars and the Γ parameter as functions of time for the fiducial model.Even in the model with magnetic field, the mass accretion rate onto the secondary star is higher than that onto the primary star in consistent with the previous hydrodynamical simulations (Bate & Bonnell 1997;Young et al. 2015;Matsumoto et al. 2019b).The Γ parameter exhibits time variability, but it remains positive throughout the evolu-  tion, suggesting that the binary system evolves towards twin binaries.The time-average value of Γ is 3.691.

R-contribution z-contribution Total
In Figure 21, the total mass accretion rate for both the primary and secondary stars is also plotted.The total mass accretion rate remains just below unity for most of the time, implying that the total mass of gas in the computational domain is gradually increasing.We observe that the mass of the circumbinary disk also experiences gradual growth over time.This suggests that the disk's evolution is not in a perfect steady state; however, the growth rate is notably slow.Specifically, ṀCBD /M CBD ∼ (10 −2 − 10 −3 )Ω ⋆ at t ∼ 10T ⋆ for the fiducial model, where M CBD denotes the mass of the circumbinary disk.Such a slow growth rate for the circumbinary disk was also noted in the previous hydrodynamical model (Matsumoto et al. 2019b).
Figure 22 shows the Γ parameters for four representative models, including models with and without a magnetic field and an infalling envelope.All the models shown have a positive Γ.The models with a magnetic field exhibit significant variability in Γ due to the time variability in the mass accretion rates.In contrast, the models without a magnetic field show low time variability.It is worth noting that the presence of the infalling envelope affects the evolution of the mass ratio since it provides gas that falls onto the primary or secondary stars/circumstellar disks.However, Γ remains positive in both cases with and without infalling envelopes.
The rate of change in specific angular momentum (τ −1 sam ) inside a cylinder as a function of radius and height of the cylinder.The upper and lower panels are for binary and single star models, respectively.A negative value means a decrease in the specific angular momentum.The time averages are taken in the periods depicted in Figure 15.

Implication for formation of close twin binaries
The current models offer a possible formation scenario for close binary systems within the framework of the disk fragmentation scenario.
Based on the analysis of Gaia data, El-Badry et al. ( 2019) discovered the existence of a sharp excess of equal-mass binaries (twin excess) in wide ranges of binary star masses and separations (see also Raghavan et al. 2010;Moe & Di Stefano 2017).The existence of the twin excess indicates the contribution of mass accretion from circumbinary disks, which has been reproduced by hydrodynamical simulations (Bate & Bonnell 1997;Young et al. 2015;Young & Clarke 2015;Matsumoto et al. 2019b).
Both MHD and hydrodynamical models show the same tendency of accretion rates, where binaries tend to evolve to twin binaries.In contrast to the hydrodynamical models, the MHD models exhibit considerable time-variability in the accretion rates, as well as in the parameter Γ, which is proportional to q.Despite the  minor changes of accretion rates, the MHD models are also consistent with the twin excess.Hwang et al. (2022) reported that twin binaries with wide separations are likely to have high eccentric orbits, proposing a scenario in which such wide twin binaries could form by scattering from dynamical interactions in their birth environment after accreting gas from a circumbinary disk on a ∼ 100 au scale.While this scenario explains the formation of wide twin binaries, the formation scenario of close twin binaries has been unknown.As shown in section 4.5, the magnetic field transfers angular momentum from the system, and it can make the separation of binaries small and provide a possible formation scenario for close twin binaries.Hara et al. (2021) reported the detection of a wiggling structure in an outflow as observed by ALMA.They suggest that the structure is a result of the orbital motion of the primary star, which acts as the driving source.Our simulation, as shown in Figure 5, successfully reproduces the presence of twisted outflows, with the outflow from the primary star being stronger than that from the secondary star.This stronger outflow from the primary star is likely to propagate farther and corresponds to the wiggling structure observed in the ALMA data.Twisted outflows are a natural consequence of binary stars.

Implications for observations
The Class I object BHB07-11 is a binary system showing spiral streamers at the central region and also shows an extended flat envelope according to the ALMA observations (Alves et al. 2017(Alves et al. , 2019)).The outflows are launched at the boundary between the inner dense circumbinary disk and the extended flat envelope.This structure resembles what the magnetized models reproduce.Figure 5 (upper right panel) exhibits outflows from the circumbinary disk, and the circumbinary disk is extended from the launching radius because of angular momentum redistribution due to the MRI.The present models explain not only the density structure but also the velocity structure.Alves et al. (2017) reported that the extended envelope exhibits a rotational velocity slightly higher than the Keplerian velocity.This excess in rotational velocity is approximately of the order of the sound speed, as observed in the positionvelocity diagram of H 12 CO line emission.Similarly, our fiducial model also shows faster rotation than Keplerian by approximately ∼ c s , which is not depicted in a figure, but aligns well with the observed envelope.Consequently, the extended circumbinary disk in our model could correspond to the envelope observed by ALMA Asymmetry has been observed in multiple circumbinary disks, including L1551 NE (Takakuwa et al. 2017), L1551 IRS5 (Takakuwa et al. 2020), and HD 142527 (Fukagawa et al. 2013).Matsumoto et al. (2019b) conducted hydrodynamical simulations that replicated the m = 1 asymmetry in circumbinary disks, which arises due to the gravitational torque of binary stars on the circumbinary disk (see also Shi et al. 2012;Ragusa et al. 2017).This asymmetry also manifests in our current models without a magnetic field (Figures 11e and 12e).In contrast to the hydrodynamical models, models with a magnetic field exhibit reduced asymmetry in circumbinary disks due to the MRI, because it dilutes the asymmetry.A similar effect has also been observed in simulations of protoplanetary disks containing planets (Zhu & Stone 2014), where the presence of asymmetry depends on the type of magnetic field model used; a non-ideal MHD model with ambipolar diffusion reproduces asymmetry in the disk, whereas an ideal-MHD model produces strong turbulence and an axisymmetric feature.Noble et al. (2021) also reported that a higher magnetization weakens substructures, such as lumps.Therefore, the magnetic field appears to regulate the degree of asymmetry in circumbinary disks.We also note that the MRI dilutes substructures of circumbinary disks, e.g., spiral arms and cavities, and the effect of the magnetic field likely controls the density contrast of the substructures.

Planet formation
Planets around binary stars are classified into two groups: those orbiting each star and those orbiting the pair of binary stars.The latter is known as circumbinary planets.To date, only 14 circumbinary planets have been discovered in 12 binary systems by Kepler and TESS (e.g., Doyle et al. 2011;Welsh et al. 2012).The small number of circumbinary planets discovered to date suggests a challenging environment for the formation of planets in circumbinary disks because orbital motion of binary stars generates turbulence in the disk (e.g., Pierens et al. 2020).
Our simulations demonstrate that, in the models without magnetic fields, distinct spiral arms emerge, which produce pressure bumps and vortices that potentially enhance the capture of dust particles (Barge & Sommeria 1995;Klahr & Henning 1997).In the magnetic models, MRI induces turbulence, which hinders dust settlement in a circumbinary disk.Moreover, turbulence disrupts the spiral arm structure, weakening the dust capture mechanism.To confirm the effect of magnetic fields on dust settlement and coagulation, it is necessary to perform MHD simulations that incorporate the advection and growth of dust particles (as has been done in studies of single-star cases; Tsukamoto et al. 2021).

Limitations of the present models
In the present models, the binary system accretes not only gas but also magnetic flux when considering the infalling envelope.This leads to the so-called magnetic flux problem (e.g., Shu et al. 2006;Zhao et al. 2011), wherein a protostar would end up with a much higher magnetic flux than what is typically observed in young stars.A potential solution to this problem involves the inclusion of magnetic diffusion processes, such as Ohmic dissipation, ambipolar diffusion, and the Hall effect (e.g., Wurster et al. 2018).The incorporation of these processes weakens the magnetic effects.Therefore, the real process is likely an intermediate situation between our previous hydrodynamical model (Matsumoto et al. 2019b) and the current ideal-MHD model.
Magnetic diffusion is more effective in the circumbinary and circumstellar disks than in the infalling envelope, due to the comparison of timescales of infall and magnetic diffusion (Nakano et al. 2002).In particular, the circumstellar disk is likely to be significantly influenced by magnetic diffusion due to its high density.This implies that the magnetic diffusion reduces the magnetic flux threading these disks (Tsukamoto et al. 2015), and influences the MRI turbulence and the launching of outflows (e.g., Bai & Stone 2013;Masson et al. 2016).These magnetic diffusion processes will be included in future work.
Even though our current work primarily focuses on the circumbinary disk rather than the circumstellar disks, we acknowledge several issues concerning the circumstellar disks.These disks are considerably thinner than the circumbinary disk.The disk thickness is around 0.1 for the primary circumstellar disk and about 0.05 for the secondary circumstellar disk at the most thick parts.The cell width in these regions is ∆x = 0.0058 (at the FMR grid level ℓ = 4), indicating that the thickest regions of the circumstellar disks are resolved with more than 20 cells.However, in the vicinity of the sink particles, for instance, at a point r sink away from the sink particle surface, the disk thickness reduces to approximately H ∼ 0.01.This thickness is resolved by only a few cells, and it is consistent with the estimated value of H = 0.004 − 0.007.This estimate comes from the relation H/R = c p (R/GM 1 ) 1/2 , where c p = 0.4 − 0.6 represents the sound speed as given by the barotropic equation of state assumed here.Therefore, accurately resolving the disk thickness of the circumstellar disks remains a challenging aspect of global 3D simulations.
One might speculate that the numerical resolution near the sink particles in the circumstellar disks is insufficient to fully resolve the MRI.However, in the case of the secondary circumstellar disk, as discussed in Appendix A, it remains stable against the MRI during most of the evolutionary stages.For the primary circumstellar disk, we confirm that the very thin area mentioned above is also stable against the MRI because of the thinness and strong magnetic field, and that the outflow is accelerated mainly in a region outside the thin area.While the numerical resolution may affect the turbulence in the inner regions of the circumstellar disks, it is unlikely to significantly affect the mean magnetic field structures that drive the outflows from the disks.Consequently, the outflows and angular momentum transport associated with them are likely to be reproduced properly.
The sink particles are utilized in the present models, indicating that structures finer than the sink radius are not resolved in the simulations.To mitigate the influ-ence of magnetic effects caused by the unresolved structures inside the sink radius, Ohmic dissipation is implemented to decouple the gas from the magnetic field.In actual protostars, there should be unresolved circumstellar disks, which could be around r sink ∼ 2 au in size, assuming a binary separation of a b ∼ 100 au.These small-scale disks could drive high-speed jets, as noted by Machida et al. (2008), but such phenomena are not reproduced in the present models.To incorporate these effects, employing prescriptions for a sink particle is a potential method (e.g., Cunningham et al. 2011;Grudić et al. 2021).

SUMMARY
We conducted three-dimensional MHD simulations of accreting binary systems to investigate the impact of magnetic fields on circumbinary materials.We varied the strength of the magnetic fields and compared the results between models with and without magnetic fields, as well as between binary and single-star models.Additionally, we examined the effects of infalling envelopes.Our findings are summarized as follows: 1.The binary models produce twisted twin outflows with high velocity from the circumstellar disks around the binary stars and a wide outflow with low velocity from the circumbinary disk.The twisted structure is attributed to the orbital motion of the binary stars.The angular momentum redistribution by MRI expands the circumbinary disk.These characteristic structures reproduced in the magnetic field models have recently been observed by ALMA.
2. A circumbinary disk in a binary system is more turbulent (α R,turb ∼ α M,turb ∼ 0.5 for the fiducial model) compared to a circumstellar disk in a single-star system.Turbulence is driven by MRI and is intensified by disturbances from spiral arms, which are caused by the orbital motion of binary stars.We also confirmed that accretion from the infalling envelope contributes to the turbulence in the circumbinary disk.
3. The angular momentum is transferred both in radial and vertical directions.In the radial direction, MRI turbulence and gas flow associated with spiral arms are the main drivers of angular momentum transfer in the circumbinary disk.Specifically, MRI turbulence reduces the specific angular momentum of the circumbinary disk.In the vertical direction, both outflows and magnetic braking transfer angular momentum, also reducing the specific angular momentum.4.Even in the MHD cases, the primary star accretes more gas than the secondary star, as observed in previous hydrodynamic models.This suggests that binary systems tend to evolve into twin binaries.In the MHD models, accretion rates exhibit variability.
We would like to thank James M. Stone for fruitful discussions and comments.Numerical computations were carried out on XC50 (ATERUI II) at Center for Computational Astrophysics, National Astronomical Observatory of Japan.This research was supported in part by the Hosei Society of Humanity and Environment, and by JSPS KAKENHI Grant Numbers JP18H05437, JP17K05394, JP23K03464.

A. NUMERICAL RESOLUTION FOR MRI
The development of the MRI in a numerical simulation hinges on the numerical resolution (e.g., Sano et al. 2004).To probe the numerical resolution requisite for accurately reproducing the MRI, we evaluated its characteristic wavelengths and the MRI quality factors (Noble et al. 2010;Hawley et al. 2011;Shiokawa et al. 2012) for the circumbinary disk (CBD), the primary circumstellar disk (PCSD), and the secondary circumstellar disk (SCSD).
The MRI wavelengths are characterized by the following equations: ) pertaining to the z and φ directions, respectively (Noble et al. 2010;Hawley et al. 2011;Shiokawa et al. 2012).
The term λ MRI,z approximately corresponds to the most unstable wavelength in the initial conditions where the magnetic fields align in the z-direction.The MRI quality factors represent the ratios of the characteristic MRI wavelengths, λ MRI , to the numerical cell sizes.They are defined as Q z = λ MRI,z /∆x and Q φ = λ MRI,φ /∆x (Noble et al. 2010;Hawley et al. 2011;Shiokawa et al. 2012).Note that the cell width ∆x depends on location because of the fixed mesh refinement.Figures A1 and A2 depict the MRI wavelengths (λ MRI,z and λ MRI,φ ) and the MRI quality factors (Q z and Q φ ) throughout the evolution for the CBD, PCSD, and SCSD of the fiducial model.As regions of the disks, we consider the regions with ρ > ρ disk (= 10ρ 0 ), and z ∈ [−0.05, 0.05], for SCSD, (A5) where (x 1 , y 1 ) and (x 2 , y 2 ) represent the positions of the primary and secondary stars, respectively, in the x − y plane.The volume averages of λ MRI and Q in these regions are shown in Figure A1 and in the left panel of Figure A2, while the azimuthally and vertically averaged cylindrical radial distribution of Q is presented in the right panel of Figure A2.
For the CBD, at the initial stage, the CBD has ⟨λ MRI,z ⟩ = 0.06, which increases to values between ∼ 1 and 2. The vertical extent of the CBD lies in the range z ∈ [−1, 1] (see Figure 3).Consequently, the CBD undergoes influence from the MRI throughout its evolution.The MRI quality factor starts with ⟨Q z ⟩ = 3.3 at the onset, but it quickly exceeds 10.Since the required Q z value to resolve the characteristic wavelength is Q z ≳ 6 (Sano et al. 2004), the MRI in the CBD is marginally resolved in the early stages of t ≲ T ⋆ .Following this period, both ⟨Q z ⟩ and ⟨Q φ ⟩ evolve within the 10-100 range (left panel of Figure A2), and they exhibit spatially constant distributions (right panel of Figure A2), effectively resolving a typical MRI mode across the entirety of the CBD..
The PCSD begins with a value of ⟨λ MRI,z ⟩ = 0.002, which subsequently increases to values between ∼ 0.05− 0.2.This is comparable to the disk thickness of z ∈ [−0.1, 0.1], implying that the PCSD is also influenced by the MRI.Initially, the MRI quality factor is ⟨Q z ⟩ = 0.3.However, by the early phase at t = 0.8T ⋆ , it exceeds the threshold value of approximately 6. Subsequently, ⟨Q z ⟩ evolves within the range of 9 − 40.The radial distribution of ⟨Q z ⟩ indicates that it exceeds the threshold value of 6, but drops to around 4 at R/a b ∼ 0.1 − 0.2.Consequently, a typical MRI mode is resolved across almost the entire PCSD, although it is marginally resolved at R/a b ∼ 0.1 − 0.2.Volume averages of MRI wavelengths for the vertical and azimuthal directions as a function of time for the circumbinary disk (CBD), the primary circumstellar disk (PCSD), and the secondary circumstellar disk (SCSD) in the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).
For the SCSD, the initial value of ⟨λ MRI,z ⟩ is 0.06, consistent with the disk thickness of the SCSD ∼ 0.05.This value exhibits a rapid increase in the very early stages, rising beyond ∼ 1, which suggests that the SCSD becomes stable against the MRI.Due to these elongated wavelengths, the MRI quality factors are notably high.

B. ESTIMATE OF ALPHA PARAMETERS
The α-parameters for the Reynolds and Maxwell stresses are computed using the fluctuating components of the velocity and magnetic field, represented by v ′ and B ′ , respectively.The calculation is performed as follows.
All physical variables were transformed from the laboratory frame to the rotating frame with an angular velocity of the binary stars, Ω ⋆ , using the equation: where U = (ρ, v, B) denotes a state vector and R(θ) denotes the rotation operator with angle θ.This transformation allows the binary stars to be observed as stationary in the rotating frame.In this frame, time averages are performed for all physical variables using the equation: The MRI quality factors for the circumbinary disk (CBD), the primary circumstellar disk (PCSD), and the secondary circumstellar disk (SCSD) in the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).The left panel shows the volume average of the MRI quality factor for each disk as a function of time.The right panel shows the cylindrical radial distributions of the MRI quality factors at t = 7 T⋆ (corresponding to the stage shown in the third row of Figure 2).For the PCSD and SCSD, the centers of the radial distributions coincide with the positions of the respective sink particles.The legend in the right panel is common for both the left and right panels.
time-averaged values from the rotating-frame variables: where we use the Favre average v, defined as because we consider compressible flows in our simulations.
The time-averaged stress components related to angular momentum transport, v ′ R v ′ φ (the Reynolds stress) and B ′ R B ′ φ (the Maxwell stress), are computed using the fluctuating parts defined by Equations (B8)-(B9), as follows: where integration along the z-direction is performed for a region with density larger than ρ disk : We set ρ disk = 10ρ 0 , which traces the surfaces of the circumbinary disk.Equations (B11) and (B12) are α parameters contributed by turbulence.We also define the α parameters of mean flow and mean magnetic field as follows, where j = Rv φ is a specific angular momentum, and T zφ = ρv z v φ , (C20) where F R (R) is the angular momentum flux through the side of an enclosed cylinder with a radius of R and a range of [z min , z max ], and F z (z) is the angular momentum flux through the base of an enclosed cylinder with a radius of R max and a height of z.
We also define a gravitational torque component as a function of R from the right-hand side of Equation (C23).This is given by S g (R) =  The outflow contribution (green line in the lower panel of Figure 18) is estimated by considering the T zφ term in Equation (C25), but taking into account the cells with v z > 0 in the z > 0 region and v z < 0 in the z < 0 region when evaluating the integration.In the scale of the circumbinary disk, contribution of S g (R) is much less than F r (R), and we only show F r (R) in Figure 18 (upper panel) for simplicity.Other contributions were calculated in a similar method.

Figure 1 .
Figure 1.The schematic diagram illustrates the computational domain.The left panel depicts a 3D view of the computational domain and the imposition of the boundary condition.The four right panels display a static hierarchical grid configuration that overlays the density distribution in cross sections.The gray grid represents the FMR blocks, with each block consisting of 16 3 cells.Black circles indicate the positions of sink particles, which model binary stars.The radii of the circles are two times larger than the actual sink radius for ease of visual confirmation.

Figure 2 .
Figure2.Evolution of the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope) at t = 0 (the initial condition), 5T⋆, 7T⋆, and 10T⋆ from left to right.The colors show the logarithmic density distributions in the z = 0 plane (upper panels) and the y = 0 plane (lower panels).The black circles represent the sink particles, with the right and left particles corresponding to the primary and secondary stars, respectively.The radii of the circles are equal to the sink radius, r sink .The entire computational domain is shown.An animation is available in the HTML version.

Figure 3 .
Figure 3. Same as Figure 2, but the magnifications of the central region of [−3, 3] 2 .An animation is available in the HTML version.

Figure 4 .
Figure 4. Evolution of the outflows for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope) at t = 0 (the initial condition), 5T⋆, 7T⋆, and 10T⋆ from left to right.The colors indicate the distributions of vz in the z = 0 plane.The lower panels are the magnifications of the central region of [−1.5, 1.5] 2 .The black circles show the sink particles; the right and left sink particles correspond to the primary and secondary stars, respectively.The radii of the circles are the sink radius, r sink .

Figure 5 .
Figure 5. 3D views of the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope) at t = 10T⋆.The volume rendering displays the gas density distribution on the logarithmic scale.The bipolar structures in red and blue indicate the velocity distributions of the vertical component (vz), revealing that the outflows from the circumstellar disks of the primary and secondary stars propagate in both positive and negative z-directions.The upper-right panel displays the same view as the upper-left panel, but with magnetic field lines depicted as gray tubes.The upper panels show the entire computational domain of [−12, 12] 3 , while the lower panels are magnifications around the binary stars, showing the region of [−3, 3] 3 .the α-parameters contributed by the mean flow ρ v R v φ and the mean magnetic field B R B φ , which are denoted by α R,mean and α M,mean , respectively (see Appendix B).The mean flow of ρ v R v φ has a large amplitude that follows the shape of the spiral arms, as shown in.Along the spiral arms, the peak value is α R,mean ∼ 2, while it has large negative values α R,mean ∼ −10 in inter-arm regions.This suggests that the spiral arms efficiently transport angular momentum outward in the spiral arms and inward in the inter-arm regions.The gravitational torque from the binary stars generates the spiral arms and is responsible for the associated angular momentum transport.This transport is also seen in hydrodynamical simulations byMatsumoto et al. (2019b) and MHD simulations byShi et al. (2012), where gas exhibits expansion motion along the spiral arms and infall motion in the inter-arm regions.These motions have been observed by ALMA for protobinary systems L1551 NE(Takakuwa et al. 2014(Takakuwa et al. , 2017) ) and L1551 IRS 5(Takakuwa et al. 2020).

Figure 6 .
Figure 6.Evolution of the plasma β for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope) at t = 0 (the initial condition), T⋆, 5T⋆, and 10T⋆ from left to right.The colors show the plasma β distributions on the logarithmic scale in the z = 0 plane (upper panels) and the y = 0 plane (lower panels).The black circles show the sink particles; the right and left sink particles correspond to primary and secondary stars, respectively.The radii of the circles are the sink radius, r sink .The central region of [−3, 3] 2 is shown.

Figure 7 .
Figure 7. Face-on view (left panel) and edge-on view (right panel) of the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope) at t = 10T⋆.The volume rendering shows the gas density distribution on the logarithmic scale.The gray tubes show the magnetic field lines.The left panel shows a region of x, y, z ∈ [−6, 6] 2 × [−2, 2], and the right panel shows a region of [−6, 6] 3 .

Figure 8 .
Figure 8.The α-parameters of the turbulent components of the Reynolds stress α R,turb (left panels) and Maxwell stress α M,turb (right panels) for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).They correspond to the αparameters for ρv ′ R v ′ φ and B ′ R B ′ φ components, respectively.The lower panels are a magnification of the upper panels.The time average is taken in the period of t ∈ [7T⋆, 10T⋆].The density of the disk surface is set at ρ disk = 10ρ0.The black circles show the sink particles; the right and left sink particles correspond to the primary and secondary stars, respectively.outflows at t = 10T ⋆ , because the outflow growth rate is slower compared to the fiducial case.Consequently, they extend to the upper and lower boundary surfaces (z = ±12) at a later stage of t ≃ 20T ⋆ .Panels (f)-(g) of Figures 11 and 12 show the binary models without infalling envelopes.Comparing models with the same magnetic field strength (e.g., Figures 11band 11f; Figures11e and 11g), we observe that the models without an infalling envelope have extended circumbinary disks because there is no ram pressure from the infalling gas.The models with infalling envelopes exhibit more disturbed density distributions in the circumbinary disks compared to those without, as seen when comparing Figures 12b and 12f, which have the same initial magnetic field strength of B z,0 = 0.1.This difference is mainly due to two factors.The first is difference in magnetic flux between the models with and without an infalling

Figure 9 .
Figure 9.The α-parameters of the mean flow/field components of the Reynolds stress αR,mean (ρ vR vφ component) (left panels) and the Maxwell stress αM,mean (BRBφ component) (right panels) for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).The time average is taken in the period of t ∈ [7T⋆, 10T⋆].The color scales are different on the left and right panels.

Figure 10 .
Figure 10.The distribution of α-parameters as a function of radius for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).The radial distribution is obtained by averaging the α-parameters in the azimuth direction.The time average is taken in the period of t ∈ [7T⋆, 10T⋆].The upper panel shows the positive values of α, while the lower panel shows the negative values.

Figure 11 .Figure 12 .
Figure 11.Surface density distributions for all the models at t ∼ 10T⋆.The whole computational domains are shown.The model parameters are indicated in each panel, with the terms "w/ infall" and "w/o infall" referring to the models with and without an infalling envelope, respectively.

Figure 13 .
Figure 13.Density distributions in the y = 0 plane (the meridional plane) for all the models at t ∼ 10T⋆.The whole computational domains are shown.

Figure 14 .
Figure 14.Density distributions in the y = 0 plane (the meridional plane) for all the models at t ∼ 10T⋆.Each panel shows the magnification of Figure 13.

Figure 15 .Figure 16 .
Figure 15.The α-parameters of the turbulent component the Reynolds stress α R,turb (ρv ′ R v ′ φ component) for all the models.

Figure 18 .
Figure 18.Angular momentum fluxes for the the fiducial model (q = 0.2, Bz,0 = 0.1 with the infalling envelope).The upper panel shows the angular momentum flux in the R-direction FR(R) as a function of R. The flux is measured through the side of a cylinder with a radius of R and a height of z ∈ [−2, 2].The lower panel shows the net flux of the angular momentum in the z-direction Fz(z) − Fz(−z) as a function of z.The net flux is measured through the top and bottom surfaces of a cylinder with a radius of 12 and height of ±z.The time average is taken in the period of t ∈ [7T⋆, 10T⋆].A positive flux means that angular momentum is extracted from the region of interest (outgoing flux).Each line represents a contribution of each component; e.g., the solid blue line shows a total contribution of the Reynolds stress TRφ, while the dashed blue line shows the contribution of turbulent component of the Reynolds stress.The methods for evaluating the contributions are described in Appendix C.

Figure 19 .
Figure19.The rate of change in specific angular momentum (τ −1 sam ) inside a cylinder as a function of radius and height of the cylinder for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).The each line show the contribution of each component to τ −1 sam (see Appendix D).The time average is taken in the period of t ∈ [7T⋆, 10T⋆].

Figure 21 .
Figure21.Mass accretion rate of the primary and secondary stars as a function of time for the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).The mass accretion rates of stars are normalized by the mass accretion rate of the envelope at t = 0 (gas injection rate at the boundary surfaces).The rate of change in the mass ratio Γ is also shown.

Figure 22 .
Figure 22.The rate of change in the mass ratio Γ as a function of time for the representative binary models.
Figure A1.Volume averages of MRI wavelengths for the vertical and azimuthal directions as a function of time for the circumbinary disk (CBD), the primary circumstellar disk (PCSD), and the secondary circumstellar disk (SCSD) in the fiducial model (q = 0.2, Bz,0 = 0.1, with the infalling envelope).
B7)    where the period of the time average is set at t ∈ [7T ⋆ , 10T ⋆ ] for the fiducial model.With the time-averaged variables, we can extract the fluctuating parts due to turbulence by subtracting the we also define α-parameters for total stressα R,total = α R,mean + α R,turb ,(B16)α M,total = α M,mean + α M,turb .(B17) C. ANGULAR MOMENTUM TRANSPORT THROUGH AN ENCLOSED CYLINDERWe consider a cylindrical region with a radius of R max and a range of z ∈ [z min , z max ].The angular momentum transport in the R and z directions are expressed by zφ + M zφ )] = Rρg φ , (C18)

F
C21)M zφ = −B z B φ /4π. is applied to Equation (C18), and an integral form is obtained, given by in the left-hand side of Equation (C23) denotes the rate of change of the angular momentum of the region of interest.We define the fluxes of angular momenta in the R and z directions as functions of R and z, respectively, from the second and third terms in the left-hand side of Equation (C23).These are given by (T zφ + M zφ ) ,(C25) where S g (R) is the gravitational torque acting on the region of interest of the cylinder with the radius R and the height, [z min , z max ], although it is not an angular momentum flux.With Equations (C23)-(C26), change in the angular momentum is expressed by∂J ∂t + F R (R max ) + F z (z max ) − F z (z min ) = S g (R max ),(C27) where J = ρjdV .

FdzR 2 TF
Figure 18 displays several contributions to the fluxes separately.The line of T Rφ represents the contribution of the T Rφ term in equation (C24), i.e., Rφ .(for T Rφ component) (C28) The contributions of M Rφ , T zφ , and M zφ are calculated in the same manner.The turbulent contributions, which are shown by the dashed lines in Figure 18, is calculated, taking into account only the fluctuating components of the stresses, e.g., 2 ρv ′ R v ′ φ .(for turbulent contribution of T Rφ )(C29) D. CHANGE IN MEAN SPECIFIC ANGULAR MOMENTUMChange in the mean specific angular momentum in the region of interest is give by angular momentum and mass inside a cylinder with the radius R max and the height [z min , z max ], and they are given by, j = Rv φ , and the volume integration dV is performed over the region of interest.The time derivative of mass ∂M/∂t is given by Gauss' law as follows, without contribution of S g , and Equations (D31)-(D33), Equation (D30) can be evaluated.When the timescale of the specific angular momentum is defined as τ sam , the rate of change in the specific angular momentum is given by 1 (D34) to the models in and 20, both the radius and height of the cylinder were set to R, and τ −1 sam was calculated as a function of R. In Figure19, τ−1 sam are shown for several contribution separately.The line of "T Rφ + R-flow" means the rate of change in specific angular momentum due to T Rφ and radial flow.It is calculated based on Equation (

Table 1 .
Model parameters