Magnetohydrodynamic Model of Late Accretion onto a Protoplanetary Disk: Cloudlet Encounter Event

Recent observations suggest late accretion, which is generally nonaxisymmetric, onto protoplanetary disks. We investigated nonaxisymmetric late accretion considering the effects of magnetic fields. Our model assumes a cloudlet encounter event at a few hundred au scale, where a magnetized gas clump (cloudlet) encounters a protoplanetary disk. We studied how the cloudlet size and the magnetic field strength affect the rotational velocity profile in the disk after the cloudlet encounter. The results show that a magnetic field can either decelerate or accelerate the rotational motion of the cloudlet material, primarily depending on the relative size of the cloudlet to the disk thickness. When the cloudlet size is comparable to or smaller than the disk thickness, magnetic fields only decelerate the rotation of the colliding cloudlet material. However, if the cloudlet size is larger than the disk thickness, the colliding cloudlet material can be super-Keplerian as a result of magnetic acceleration. We found that the vertical velocity shear of the cloudlet produces a magnetic tension force that increases the rotational velocity. The acceleration mechanism operates when the initial plasma $\beta$ is $ \beta \lesssim 2\times 10^1 $. Our study shows that magnetic fields modify the properties of spirals formed by tidal effects. These findings may be important for interpreting observations of late accretion.

Protoplanetary disks are natural byproducts of the formation of stars through the gravitational collapse of rotating gas clouds. Direct three-dimensional (3D) magnetohydrodynamic (MHD) simulations of gravitational collapse of molecular cloud cores revealed several fundamental processes of formation and growth of disks (e.g. Machida et al. 2007;Hennebelle & Ciardi 2009;Tomida et al. 2013) and driving jets and outflows (e.g. Banerjee & Pudritz 2006;Machida et al. 2008). Previous studies have found that magnetic fields play essential roles in removing angular momentum from accreting gas at different evolutionary stages and radii (e.g. Mouschovias & Spitzer 1976;Tomisaka 2002;Hennebelle & Fromang 2008). The deceleration of rotating gas by magnetic tension is called magnetic braking.
Theoretical studies of star and disk formation have been progressing along with advances in numerical simulations. Spherically symmetric models were used to study the gravitational collapse of molecular clouds (e.g. Larson 1969;Masunaga & Inutsuka 2000). The formation of protoplanetary disks has been investigated using a gravitationally unstable, isolated rotating molecular cloud with a spherically symmetric density structure as the initial condition (e.g., see reviews by Inutsuka 2012;Li et al. 2014). With this type of setup, accretion onto the disks occurs nearly in an axisymmetric manner and becomes monotonically weaker with time. Numerical simulations that relax the assumption of axisymmetry have been performed to study the effects of the misalignment between the rotational axis and the background magnetic fields in an axisymmetric density structure (e.g. Matsumoto & Tomisaka 2004;Joos et al. 2012;Tsukamoto et al. 2015). Recently, more complicated star and disk formation processes have been investigated that consider the effects of turbulence on the core scale or larger (e.g. Joos et al. 2013;Lee & Hennebelle 2016;Kuffmeier et al. 2017;Lam et al. 2019).
Heterogeneous late infall will commonly occur because star-forming regions are intrinsically inhomogeneous. Starforming molecular clouds must be turbulent (see the review by Hennebelle & Falgarone 2012). Kuffmeier et al. (2017) performed 3D MHD simulations of star formation starting from the giant molecular cloud (GMC)-scale and showed that the accretion rate onto the disk significantly varies with time, depending on the condition of the surrounding environment. They also found accretion of gas clumps or "cloudlets" onto a disk (Kuffmeier et al. 2018). Star-disk systems formed in turbulent regions will naturally have a finite relative velocity to the ambient gas, which is another important factor to determine the accretion rate. The accretion rate becomes higher when star-disk systems enter higher density regions or with an increase in the infall rate. The importance has been extensively studied based on Bondi-Hoyle-Lyttleton accretion (Padoan et al. 2005;Throop & Bally 2008;Moeckel & Throop 2009;Scicluna et al. 2014;Wijnen et al. 2016). From the observational census toward NGC 3603, Beccari et al. (2010) argued the necessity of late infall in old (>10 Myr) pre-main sequence stars with disks.
Some observations support the idea of the heterogeneous accretion onto disks. Tail structures connecting to disks with a size of 100-1000 au are found in several systems such as AB Aur (Nakajima & Golimowski 1995;Grady et al. 1999), HD 100546 (Ardila et al. 2007), Z CMa (Nakajima & Golimowski 1995;Liu et al. 2016), SU Aur (Akiyama et al. 2019;Ginski et al. 2021), and RU Lup (Huang et al. 2020). These asymmetric structures may be the result of late infall from the remnants of envelopes or giant molecular clouds. The sulfur monoxide (SO) emissions provide evidence that asymmetrically infalling gas forms accretion shocks around disks (Sakai et al. 2016;Garufi et al. 2021). Those observations suggest that late infall onto disks will be common and that disks will be subject to asymmetric accretion. Moreover, connection to the origin of the misalignment of inner and outer disk regions has been discussed observationally (Ginski et al. 2021), which could be consistent with hydrodynamic models (Thies et al. 2011;Kuffmeier et al. 2021).
Motivated by observations and simulations, 3D hydrodynamic simulations of late infall onto protoplanetary disks have been performed to study the detailed process of late encounter events. Dullemond et al. (2019) investigated the late encounter between a star and a cloudlet on a scale of several thousand au and showed the encounter can lead to the formation of arc or tail-like structures. Kuffmeier et al. (2020) also performed a set of hydrodynamic simulations of the cloudlet encounter on a similar scale and studied the formation of second-generation disks. These 1000 au scale simulations highlight how the stellar gravity captures GMC scale gas. However, detailed interaction between infalling gas and pre-existing protoplanetary disks on a 10-100 au scale remains unclear. Magnetic fields play important roles at this scale because the cloudlets strongly bend and amplify magnetic fields during the encounter and increase the importance of magnetic tension. Because the encounter process at this scale determines the mass and angular momentum supply to the disks, detailed investigations based on MHD models are necessary.
By performing a set of 3D MHD simulations, we investigate the magnetic effects during the cloudlet encounter event on several 100 au scale. The rest of this paper is structured as follows: Section 2 describes our model setup. The disk and cloudlet structures are explained. Section 3 presents numerical results. The role of magnetic tension is clarified by varying the initial size of a cloudlet. We explain that magnetic tension can not only decelerate but also accelerate the rotation motion of infalling gas depending on the size of the encountering cloudlet. In addition, a brief comparison to a hydrodynamic model is presented. Section 4 will briefly discuss the magnetic field strength of the cloudlet required for magnetic acceleration and the impact of the non-ideal MHD effect. Section 5 summarizes our results.

Basic Equations
We performed 3D MHD simulations to examine the nonaxisymmetric accretion process of a magnetized gas clump (cloudlet) onto a protoplanetary disk (as shown in Section 2.2). It is assumed that the cloudlet and disk consist of cold molecular gas. The cold cloudlet and disk are surrounded by a warm neutral atomic gas. We solve the following ideal MHD equations, where ρ, p, v, B and Φ denote the mass density, the thermal pressure, the velocity vector, the magnetic field vector, and the gravitational potential by the protostar, respectively. The equation of state is where T , k B , and m H denote the temperature, the Boltzmann constant, and the mass of a hydrogen atom, respectively. We assume that the gas is nearly isothermal, considering that the thermal relaxation timescale is significantly shorter than the dynamical timescale (e.g. Aota et al. 2015). For this reason, we adopt an ideal gas that has a specific heat ratio, γ = 1.05. We ignore explicit heating and cooling in this study. Non-ideal MHD effects such as ambipolar diffusion are also ignored. We will discuss the validity of this assumption in section 4.2.
The source of gravity in our model is only the (spatially unresolved) central protostar with a mass of M = 0.5M . We ignore the self-gravity of the gas disk and the cloudlet. The gravitational potential of the protostar is softened artificially within a radius of a s from the protostar. The functional form is expressed as in the cylindrical coordinates, (r, ϕ, z), with the origin at the central protostar. A smaller softening length a s requires a larger computational resources. Due to the limitation of our resources, we take a s = 100 au. The softening of the gravitational potential introduces some artificial structures including a ring-like structure around r = a s in the density map. However, we argue in Appendix A that the cloud-disk interaction is not significantly affected by the softening.

Warm neutral medium
The warm neutral medium is assumed to be in hydrostatic equilibrium under the gravitational potential given by Equation (8). The pressure and density distributions outside the disk and the cloudlet are given as  where p 0 and ρ 0 denote the pressure and density at a large distance from the star, respectively. We take ρ 0 = 2.0 × 10 −19 g cm −3 . The temperature and mean molecular weight of the warm neutral medium are assumed to be T 0 = 410 K and µ 0 = 1.27, respectively. Substituting these values into Equation (7), we obtain p 0 = 5.3×10 −9 erg cm −3 .

Protoplanetary disk
The disk is set in the region of r ≤ R d and |z| ≤ z s (r), where R d and z s (r) denote the disk outer radius and the disk upper boundary, respectively. We also set the length unit r 0 as 100 au. The upper boundary is defined as z s (r) = 2 5 r 2 0 + r 2 tanh 3 2 The disk outer radius is set to R d = (5/2)r 0 = 250 au. The disk is mildly flared and has the maximum height, z s,max = 66 au at r = 163 au. The disk has a uniform temperature of T d = 49 K, if the mean molecular weight is µ d = 2.3. To ensure the pressure balance with the warm neutral medium, we impose the following boundary condition: Combining the hydrostatic equations and Equation (12), we obtain Substituting Equation (13) into the radial component of the hydrostatic equation, we obtain The rotational velocity is slightly lower than the Keplerian velocity because the disk is supported against the gravity in part by the gas pressure. The disk mass is M d = 1.6 × 10 −3 M .

Cloudlet
We consider an ellipsoidal cloudlet accreting onto the protoplanetary disk though a warm neutral medium. The initial position of its center is located at (x, y, z) = (5r 0 , 0, 0). The surface of the cloudlet is defined as where a x , a y and a z are half-lengths of the principal axes in the x, y and z directions, respectively. The cloudlet has a uniform temperature of T c = 37 K if the mean molecular weight is set to µ c = 2.3. The cloudlet is in the pressure balance with the warm neutral medium. The pressure and density distributions inside the cloudlet are described as The initial cloudlet mass, M c is listed in Table 1. The velocity distribution inside the cloudlet is determined as follows. The initial velocity of the cloudlet is set by assuming that the orbit is parabolic with a periastron distance, d = 100 au ( Figure 1). The cloudlet has a uniform specific angular momentum l z = √ 2GM d = 4.5 × 10 20 cm 2 s −1 such that the pericenter is d. This requirement determines the velocity distribution of the cloudlet ( The initial magnetic field is assumed to be where B 0 denotes the initial magnetic field strength. Only the cloudlet is threaded by a straight magnetic field (parallel to the z axis), and the disk is unmagnetized ( Figure 1). If we initially had imposed a magnetic field to the disk, the disk structure would have been largely affected by the magnetic field via e.g. magneto-rotational instability (e.g. Velikhov 1959;Balbus & Hawley 1991). As we wish to focus on the dynamical interaction between the infalling cloudlet and the disk, we ignore the disk magnetic field to simplify the situation. The plasma β of the cloudlet is ∼ 10 so that the magnetic pressure is minor in and around the cloudlet. Therefore, the cloudlet expansion due to the magnetic pressure is insignificant. We conducted simulations for eleven cloudlet models of various sizes as summarized in Table 1. The model name consists of the half size of the cloudlet, a z , and the initial field strength. The values of a x and a y are fixed to 150 au. As shown later, the relative size between the cloudlet and the disk thickness, a z /z s,max is a crucial parameter for the resulting rotational velocity profile.
We describe two major limitations: the timescale and hysteresis of our models. The dynamical interaction between the cloudlet and the disk occurs on a much shorter timescale (∼ 10 3 yr) than the viscous timescale (∼ 10 6 yr). Considering this, we ignore the details associated with disk accretion, such as an effective viscosity and a disk magnetic field. Our models are therefore incapable of studying the long-term evolution. In addition, our models focus on the single cloudlet encounter event, though cloudlet capture events are likely to occur repeatedly during the viscous timescale (see, e.g., Kuffmeier et al. (2018)). Our model cannot deal with the situations where the hysteresis of previous events is significant. Note-az/zs,max denotes the relative size between the cloudlet and the disk thickness.

Numerical Methods
We solved Equations (1) through (6) using CANS+ (Matsumoto et al. 2019). The basic equations are solved in the Cartesian coordinates, (x, y, z). We adopt the Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver of Miyoshi & Kusano (2005) and the hyperbolic divergence cleaning method (Dedner et al. 2002). We employ MP5 (Suresh & Huynh 1997) and a third-order Runge-Kutta method to achieve the fifth-order accuracy in space and third-order accuracy in time.
The computation domain is a rectangular box that covers |x| < L x , |y| < L y and |z| < L z , where L x , L y , and L z take different values depending on the models. The values of L x , L y and L z are listed in Table 1. We fix the spatial resolution for all the models, and ∆x = ∆y = ∆z = 7 au. For simplicity, we apply the fixed boundary conditions to all the physical variables (ρ, p, v, B) in all the directions. The fixed boundary condition for the magnetic field may not be natural. We discuss the effects of the fixed boundary condition in Appendix B.

RESULT
We first present the results of two models with and without magnetic fields, S60 B116 and S60 NB, respectively, to highlight the importance of magnetic fields. We then show how the impact of magnetic fields depends on the cloudlet size and the cloudlet field strength.

Models with and without magnetic fields
We compare the results of Models S60 B116 and S60 NB, which only differ in the presence of magnetic fields. In both models, the relative size of the cloudlet to the disk thickness is close to unity (0.91). The initial plasma β at the center of the cloudlet of S60 B116 is 14, where the plasma β is defined as the ratio of the gas pressure to the magnetic pressure. The field strength is B 0 = 116 µG. Figure 2 displays the evolution of the density (left) and rotational velocity (right) distributions on the midplane for the hydrodynamic model, S60 NB. The rotational velocity is normalized by the local Keplerian orbital velocity, v Kep = GM/r. At t = 1126 yr, an outer part of the disk is broken by the cloudlet. After the encounter, a part of the colliding cloudlet material is ejected away from the disk, forming a spiral structure as shown in the panels of t = 3378 yr and 5630 yr. A large fraction of the gas in the spiral falls back onto the disk. A diffuse remnant of the ejected material remains outside the disk. The rotational velocity is not higher than the local Keplerian orbital velocity except for a narrow region of the cloudlet impact in an early stage (t = 1126 yr). Figure 3 is the same as Figure 2; however, for model S60 B116, the spiral arm forms but does not grow in size in model S60 B116. The growth is likely to be suppressed by the magnetic tension force acting on the cloudlet. Gas x [au] ejection is not appreciable, and accordingly, the fallback of the ejecta is insignificant. Instead, there appears an arc-like region where the rotational velocity is half of the local Keplerian orbital velocity at t = 5630 yr. The reduction is a direct consequence of magnetic braking. The arc is wound up to be narrower at t = 9008 yr. When the field strength is weak (B 0 = 58 µG (S60 B58) and B 0 = 82 µG (S60 B82)), the results are almost the same with S60 B116 although the rotational velocity reduction is small. We visualized the three-dimensional structure of magnetic fields in Figure 4. Line colors denote the field strength. We also show the mass density distribution by volume rendering. Magnetic fields are dragged and stretched by the infalling cloudlet. As a result, magnetic tension force decelerates the rotational motion of the colliding cloudlet material.

Dependence on the cloudlet size
We investigate the dependence of the resulting rotational velocity profile on the cloudlet size for the models with magnetic fields.
Here, we present the results of Model S150 B116 as an example to highlight the importance of the relative size of the cloudlet to the disk thickness. In this model, the vertical cloudlet size is approximately twice as large as the disk thickness. Therefore, the top and bottom parts of the cloudlet do not collide with the disk material at the time of impact and slide on the disk surfaces. We will show that magnetic fields accelerate the rotational motion in such cases.   Figure 5 displays the midplane density (left) and rotational velocity (right) distributions at different times for Model S150 B116. A spiral structure is formed even in the presence of magnetic fields, suggesting magnetic acceleration. The rotational velocity maps indicate that the colliding cloudlet material is accelerated to a super-Keplerian velocity. A fraction of the cloudlet becomes gravitationally unbound as a result of magnetic acceleration. Figure 6 describes the 3D structure of magnetic fields that accelerate the colliding cloudlet material. Panels (a) and (b) show the birds-eye-view images of the density structure at two different times. As in the case of Model S60 B116, the magnetic fields of the cloudlet are highly stretched and amplified. However, the magnetic field geometry is different between the two models; in Model S150 B116, the magnetic tension force is operating to accelerate the disk rotational motion while in Model S60 B116 it decelerates the region. Panel (c) is an enlargement of Panel (a) but from a different viewing angle, which is shown by the black arrow in panel (a). The cloudlet and disk are colored in blue (35 ∼ 45 K) and grey (45 ∼ 55 K), respectively. Vector arrows show the normalized velocity of cloudlet on the slice at x = 3.5 au. After the cloudlet encounters the disk, the cloudlet material around the midplane is shock-compressed and decelerates in the radial direction. However, the top and bottom parts of the cloudlet do not collide with the disk body. Instead, they slide on the disk surfaces along the parabolic orbit (see panel (c) of Figure 6). Namely, the top and bottom parts retain a significant radial component of the velocity. As a result, magnetic fields are transferred toward the center more quickly around the disk surfaces than around the disk midplane. The flows sliding on the disk surfaces move toward the center to spin up, twisting up magnetic fields. The vertical velocity shear results in a magnetic field structure that accelerates the colliding cloudlet material through the magnetic tension force. This example demonstrates the importance of the relative size of the cloudlet to the disk thickness. Figure 7 shows a schematic diagram of this process. For comparison, we also investigated the non-magnetized model (S150 NB). We confirmed that the model has almost the same result as S60 NB; a part of colliding cloudlet material is ejected and forms a sub-Keplerian spiral structure. We will summarize the magnetic and non-magnetic cases in Section 5.
To demonstrate that magnetic tension force indeed accelerates the colliding cloudlet material on a timescale ∼ 10 3 yr, which is shorter than the Keplerian orbital period, we evaluated τ mag , the timescale of acceleration by magnetic tension force required for the gas to reach the local escape velocity. τ mag is given by the work rate of magnetic tension force (v r e r + v ϕ e ϕ ) · ((B · ∇) B)/4π and expressed as, τ mag can take either the positive and negative signs, and the positive and negative signs indicate the acceleration and deceleration timescales, respectively. x [au] Figure 5. Same as Figure 2, but for Model S150 B116. Figure 8 displays the midplane distributions of τ Kep /τ mag at different times, where τ Kep = 2π(r 3 /GM ) 1/2 is the Keplerian orbital period. The bottom panels show that τ Kep /τ mag is larger than two in the spiral and super-Keplerian regions where τ Kep = 7.4×10 3 yr and 1.6×10 4 yr at r = 300 au and 500 au, respectively. It means that the acceleration timescale is τ mag ≤ 3.7 × 10 3 and 8.0 × 10 3 yr at r = 300 au and 500 au, respectively. The acceleration timescale can be shorter if we take account of the preexisting kinetic energy before magnetic acceleration. Therefore, we confirmed the rapid acceleration by magnetic tension force.

Dependence of the mass of the gravitationally unbound gas on the cloudlet size
We have shown that magnetic fields can either decelerate or accelerate the rotational motion of the colliding cloudlet material, depending on the relative size of the cloudlet to the disk thickness. When the cloudlet size is larger than the disk thickness, the vertical velocity shear in the colliding cloudlet can result in the acceleration of rotational motion. As a result, a fraction of accelerated gas becomes gravitationally unbound.
To evaluate the size dependence of the mass of the gravitationally unbound gas M unbound , we investigated the time evolution of the mass of the gravitationally unbound gas in all models. The results are shown in Figure 9, where the mass is normalized by the initial cloudlet mass M c . The left panel shows the dependence of cloudlet size. The right panel shows the dependence of magnetic field strength. The mass of the gravitationally unbound gas is evaluated in the regions with the conditions of ρ ≥ 10 −18 g cm −3 and |z| ≤ h, where the definition of h depends on the relative size of the cloudlet: Approximately 1-10% of the initial cloudlet mass is found to be gravitationally unbound in Models S120 B116, S150 B116, and S180 B116. In these models, the cloudlet sizes are larger by a factor of ∼2 to 3 when compared to the disk thickness, and the field strengths are larger than 1 × 10 2 µG. In the other models, the mass of the gravitationally unbound gas is negligibly small. The above results provide constraints on the relative size of the cloudlet and the field strength of magnetic acceleration.

DISCUSSION
We performed 3D MHD simulations with different cloudlet sizes and magnetic field strengths. The simulations demonstrate that the rotational velocity of the colliding gas can either be sub-Keplerian or super-Keplerian, depending mainly on the relative size of the cloudlet to the disk thickness. When the cloudlet size is comparable to or smaller than the disk thickness, the rotation motion of the colliding cloudlet material is only decelerated by magnetic braking. However, if the cloudlet size is larger than the disk thickness, the colliding cloudlet material can rotate at a super-Keplerian velocity as a result of magnetic acceleration. We showed that the vertical velocity shear of the cloudlet produces a magnetic tension force that increases the rotational velocity (Figures 6 and 7). Our model shows that the spiral or arm can rotate at a super-Keplerian velocity as a result of magnetic acceleration ( Figure 5). Such a super-Keplerian (non-Keplerian) spiral structure is found around RU Lup, a class II object (e.g. Alcalá et al. 2017;Andrews et al. 2018;Bailer-Jones et al. 2018;Gaia Collaboration et al. 2018;Yen et al. 2018). According to Huang et al. (2020), gravitationally unbound clumps are distributed along the spiral arms. Their mass is estimated to be ∼ 0.1 − 150 M ⊕ . It has been considered that the clumpy spirals are results of gravitational instability x [au] Figure 8. The ratio of Keplerian orbital period τKep to acceleration timescale by magnetic tension force τmag at the midplane for Model S150 B116. τKep/τmag is larger than two in the spiral and super-Keplerian regions where τKep = 7.4 × 10 3 yr and 1.6 × 10 4 yr at r = 300 au and 500 au, respectively. It means that the acceleration timescale is τmag ≤ 3.7 × 10 3 and 8.0 × 10 3 yr at r = 300 au and 500 au, respectively. The acceleration can be explained by magentic tension force.
(e.g. Durisen et al. 2007;Vorobyov 2016). However, our results suggest that such structures can also be formed through the cloudlet capture event and magnetic acceleration. In Section 3.3, we showed that 1-10% of the initial cloudlet mass goes to the gravitationally unbound gas if magnetic acceleration efficiently operates (see also Figure 9). If the initial cloudlet mass is ∼ 1 × 10 −4 M , then the gravitationally unbound gas with the mass of ∼ 0.3 − 3 M ⊕ can be produced.

Munbound/Mc
S60 NB S60 B58 S60 B82 S60 B116 S150 NB S150 B58 S150 B82 S150 B116 Figure 9. Time evolution of the ratio of the mass of the gravitationally unbound gas to the initial cloudlet mass for all models.°5 In the following, we discuss some conditions under which magnetic acceleration can be expected.

Magnetic field strength required for magnetic acceleration
We estimate the lower limit of the initial field strength B 0 of the cloudlet required for magnetic acceleration. The field strength of the cloudlet is mainly amplified after the encounter via shock compression. The field strength will be further increased by the vertical velocity shear (Figures 6 and 7), however, we assume that the dominant amplification process is shock compression. For the amplified magnetic fields to produce gravitationally unbound structures, the magnetic energy density should be larger than the gravitational energy density at the accelerated region (spiral in Figure 5). Considering this, we require the relation where B post is the field strength after shock compression, and the value is B 0 times the compression ratio. The compression ratio can be obtained from the Rankine-Hugoniot relation for perpendicular shocks (e.g., Equation (5.35) of Priest (2014)). For r ∼ 400 au and ρ ∼ 5 × 10 −18 g cm −3 , the lower limit is estimated as B 0 ∼ 1 × 10 2 µG when the cloudlet has v pre ∼ 1.7 km s −1 , T pre ∼ 38 K, ρ pre ∼ 1 × 10 −17 g cm −3 just before the collision. The upper limit of the initial plasma β at the center of the cloudlet is ∼ 2 × 10 1 . The right panel of Figure 9 supports the validity of the above estimation. Among the models S150 B58, S150 B82 and S150 B116, only Model S150 B116 has a field strength larger than the lower limit. The figure shows that Model S150 B116 produces a significantly larger mass of unbound gas than the other two models. However, the magnetic field becomes strong enough before the collision of the cloudlet with the disk, magnetic braking should suppress the dynamical infall of the cloudlet. Therefore, the acceleration occurs when the field strength is stronger than the above lower limit but not strong enough to suppress the infall of cloudlet before the collision.

Non-ideal MHD effect
The protoplanetary disk is generally in a low ionization state (e.g. Gammie 1996), which makes the non-ideal MHD effects important. As magnetic diffusion results in a weak magnetic acceleration, we estimate the impact. We only consider ambipolar diffusion. In the plasma composed of ions, electrons, and neutrals, the ambipolar diffusion coefficient η AD (e.g. Wardle 2007;Masson et al. 2016) is given by, where γ AD , σ in v i , ρ n , ρ i , m n , m i are the drag coefficient, the ion-neutral collision rate, the neutral density, the ion density, the mass of neutral particle and the mass of ion particle. The ion-neutral collision rate is σ in v i ∼ 2 × 10 −9 cm 3 s −1 (Osterbrock 1961) for m n ∼ 2m H and m i ∼ 10m H . We also introduce the ionization fraction χ = n i /n n where n n and n i are the neutral number density and the ion number density. The magnetic diffusion timescale τ AD based on η AD is expected to be τ AD = L 2 η AD ∼ 5 × 10 3 yr L 10 2 au 2 B 10 2 µG −2 n n 10 6 cm −3 2 χ 10 −8 where L denotes a typical length. We take the disk thickness ∼ 100 au as a typical length L because the magnetic field lines are curved on that length scale. If the ionization fraction χ is smaller than 10 −8 , τ AD is smaller than acceleration timescale ∼ 10 3 yr. Therefore, magnetic acceleration may be suppressed. However, we note that there is large uncertainty in the magnitude of the ambipolar diffusion coefficient. The value of the diffusion coefficient is sensitive to the details of dust grains, which are not clearly understood.

SUMMARY
We investigated nonaxisymmetric late accretion onto the protoplanetary disk considering magnetic fields. We modeled the accretion as a cloudlet encounter event at a few hundred au scale. We summarize the results in the phase diagram ( Figure 10) under different conditions of cloudlet size (a z ) and magnetic fields (B 0 ) at t = 9008 yr. When the cloudlet is not magnetised (Gray region of Figure 10), a part of the cloudlet is ejected away from the disk after the cloudlet-disk collision and forms the sub-Keplerian spiral. When the cloudlet size is comparable or smaller than disk thickness (a z < 2z s,max ) and magnetised (Blue region of Figure 10), magnetic braking is effective and a small sub-Keplerian spiral is formed. When the cloudlet is larger than the disk thickness (a z > 2z s,max ) and magnetised (Yellow region of Figure 10), colliding cloudlet material is accelerated to a super-Keplerian velocity by magnetic fields.
A part of accelerated material is gravitationally unbound. The mass of the gravitationally unbound gas is regulated by the lower limit of field strength (See subsection 4.1). In earlier studies by e.g., Dullemond et al. (2019) and Kuffmeier et al. (2020), spiral structures are formed by the tidal effect. In this study, we found that a sufficiently high magnetic field causes magnetic acceleration, which can rotate the spiral at super-Keplerian velocities.
This work was supported in part by the JSPS KAKENHI grant Nos. JP18K13579, JP19K03906, JP21H04487, and JP22K14074. Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

A. DEPENDENCE OF SOFTENING LENGTH OF THE GRAVITATIONAL POTENTIAL
We describe the effect of softening length (a s ) of the gravitational potential on the disk structure and its evolution. We perform test calculations for two models with a s = 75 and 100 au in the absence of the cloudlet. The calculated period is 11260 yrs, which corresponds to approximately two rotational periods at r = 250 au. Figure 11 shows the time evolution of the surface density distribution, where one can find that both models nearly maintain the initial density distributions during the calculated periods. Therefore, the softening length does not significantly affect the disk evolution at a scale larger than the softening length.
The softening length affects the gas profile around the center. As a smaller softening length results in a deeper gravitational potential, the inner surface density is higher for the model with the smaller a s . The surface density profile outside the softening length is unchanged between the two models. We note in Figure 11 that a peak appears in the surface density around the transition region (r = a s ). This structure is an artifact produced by the artificial change in the gravitational potential and can be discerned as a ring-like structure (see panel (b) of Figure 1). Nevertheless, as the cloud-disk interaction occurs in the outer part of the disk, the inner disk structure does not affect the dynamics.
We also describe the effects on the cloud-disk interaction. We perform two calculations for models of S60 B116 (magnetic braking) and S150 B116 (magnetic acceleration) with a s = 75 au, and compare these results with the models with a s = 100 au. The magnetic field structure, which is essential for both magnetic braking and acceleration, is very similar between the cases with the two different softening lengths ( Figure 12). The mass of the gravitationally unbound gas is found to be also similar ( Figure 13). Therefore, the softening length does not significantly affect the cloud-disk interaction.

B. IMPACTS OF BOUNDARY CONDITIONS ON THE CLOUDLET MAGNETIC FIELD
The models in the main text adopt fixed boundary conditions. We describe that the top and bottom boundary conditions have little impact on the main results of this study by comparing the results with different boundary conditions. We note that the dynamics before the cloudlet collision is unaffected by the boundary conditions because the Alfvén transit timescale outside the cloudlet (approximately a few 1,000 yrs) is longer than the infall timescale. Indeed, Figures 4 and 6 show that the magnetic field remains straight near the top and bottom boundaries approximately at 1,000 yrs when the cloudlet encounters the disk.
We compare the results of S60 B116 and S150 B116 with two different boundary conditions for the top and bottom boundaries; the fixed boundary and outgoing boundary conditions. When the outgoing boundaries are adopted, the flows entering the simulation domain from the outside are forbidden (v z is set to zero and the zero-gradient boundary conditions are applied to the other variables). The boundaries allow the outgoing flows and apply the zero-gradient boundary conditions to all the variables of the flows. The 3D structures for the cases of the outgoing boundary conditions are shown in Figure 14. Comparing these results with the results for the fixed boundary cases (Figures 4 and 6), one will find that the magnetic field structures, which are essential for both magnetic braking and acceleration, are very similar to each other. We also show the mass of the gravitationally unbound gas in Figure  15. For S150 B116 Outgo (Outgoing boundary conditions), the result is almost the same as S150 B116 Fix (Fixed boundary conditions). For S60 B116 Outgo, the mass of the gravitationally unbound gas has a small peak around
10000 yr against S60 B116 Fix. However, it is less than 1% of the initial cloudlet mass, much smaller than magnetic acceleration models. Therefore, we consider that the impacts of fixed boundary conditions on the results is insignificant. Munbound/Mc S150 B116 as = 75 au S150 B116 as = 100 au Figure 13. Same as Figure 9, but for S60 B116 and S150 B116 with as = 75, 100 au.

S60_B116 S150_B116
Magnetic acceleration Magnetic braking Figure 14. 3D structure of S60 B116 (left) and S150 B116 (right) with outgoing boundary conditions at t = 3378 yr.  Figure 15. Same as Figure 9, but for S60 B116 and S150 B116 with fixed boundary conditions (S60 B116 Fix and S150 B116 Fix) and outgoing boundary conditions (S60 B116 Outgo and S150 B116 Outgo). The vertical ranges of the two panels are different.