Protoplanetary Disk Size under Nonideal Magnetohydrodynamics: A General Formalism with Inclined Magnetic Field

Many mechanisms have been proposed to alleviate the magnetic catastrophe, which prevents the Keplerian disk from forming inside a collapsing magnetized core. Such propositions include inclined field and nonideal magnetohydrodynamics effects, and have been supported with numerical experiments. Models have been formulated for typical disk sizes when a field threads the rotating disk, parallel to the rotation axis, while observations at the core scales do not seem to show evident correlation between the directions of angular momentum and the magnetic field. In the present study, we propose a new model that considers both vertical and horizontal fields and discuss their effects on the protoplanetary disk size.


INTRODUCTION
Due to advances in interferometric observations, the increasing number of embedded early phase disks now allows some interpretation of the physical process during the class 0/I young stellar object evolution phases.One major conclusion is that disks seem to be small with size not exceeding 50 AU at this stage.However, there have yet been strong conclusive results on whether the orientation of the magnetic field is correlated with the disk rotation axis (see reviews e.g.Maury et al. 2022;Tsukamoto et al. 2023, and references therein).The disk size at this stage is regulated by its magnetic interaction with the accreting envelope, and some scaling relations has been suggested assuming aligned field and rotation (Hennebelle et al. 2016).The ambipolar diffusion has been concluded to be the major effect regulating the disk size at longer timescales (Lee et al. 2021;Zhao Corresponding author: Yueh-Ning Lee ynlee@ntnu.edu.twet al. 2020).Numerical simulations suggested no correlation at the core scale (Kuznetsova et al. 2020).It is therefore important to understand how the inclination angle between the angular momentum and global magnetic field affects the formation of the disk.In this current work, we propose a model that describes the disk size as function of mass, magnetic field strength, and the misalignment angle, as a result of ambipolar diffusion.

Properties of the envelope and the disk
The prescriptions follow exactly the same as that in Lee et al. (2021).The density in the collapsing envelope is proportional to the singular isothermal sphere (SIS Shu 1977) where r is the radial distance, C s = 200 m/s is the thermal sound speed, G is the gravitational constant, and δ ρ is a numerical factor of the order of unity.The density Lee et al.
right inside the envelope-disk boundary is amplified by the accretion shock, such that with M being the Mach number.The Keplerian velocity is given by where M is the total mass of the system, with M * and M d representing the mass of the star and that of the disk.The rotational velocity u ϕ = δ ϕ v kep and infall velocity u r = δ r v kep both scale with the Keplerian rotation velocity.Vertical hydrostatic equilibrium results in the disk scale height In a disk where the Keplerian rotation velocity is significantly supersonic, it follows naturally that h ≪ r.The numerical factors of order unity δ ρ ≳ 1, and δ ϕ , δ r ≲ 1.
It can be demonstrated that the results depend only sublinearly on these factors (Lee et al. 2021) and thus they will be ignored in the following discussions for conciseness.In the present work, we discuss a general scenario where the magnetic field can be oriented in any direction.Assigning external field B 0 and field inclination angle i with respect to the disk rotation axis, we can express the characteristic value of the three components as where χ is a dimensionless coefficient for the strength of the azimuthal component.The factor 2/π in the radial component results from equally distributing the horizontal flux through the lateral cross-section of the disk into a cylindrical area, i.e., flux B 0 2rh over area πrh.As a results of symmetry, the circular integral of B ϕ along the disk edge is zero.The strength of the azimuthal field component is only a result of interplay between the magnetic flux and the disk rotation.

Nonideal MHD equations
For treating the disk properties, we present below the non-ideal magnetohydrodynamic (MHD) equations.The equation of momentum conservation writes where u is the velocity vector, P is the pressure, B is the magnetic vector, and Φ is the gravitational potential.
The induction equation writes where η A , η H , and η O are the resistivities of the ambipolar diffusion, Hall effect, and Ohmic dissipation, respectively.The Hall effect is more relevant for transient disk behaviors Zhao et al. (2020) and the Ohmic dissipation operates at higher densities.The present work focuses on the effects of the ambipolar diffusion.

Governing equations of the disk
Lee et al. ( 2021) discussed the strong-field and weakfield cases, where the diffusion of field lines is predominantly due to the magnetic tension and magnetic pressure gradient, respectively.This does not necessarily mean that the magnetic field is strong or not in absolute values, but actually indicates its strength with respect to the fluid inertia and the extent of field line winding due to induction.In reality, a growing disk quickly enters the weak field regime, where the azimuthal field generated by differential rotation is the dominant component, although the field strength could actually increase in time.
Here we rename them magnetic-dominated and inertiadominated regimes, which are more physically intuitive.The field lines tend to straighten themselves due to the magnetic tension in the magnetic-dominated regime.On the other hand, field lines tend to passively follow the rotating flow and have more winded configuration in the inertia-dominated regime.
The angular momentum is maintained by the balance between the accretion from the envelope and the braking due to the magnetic tension, such that the azimuthal component of Equation ( 8) simplifies to The azimuthal magnetic field results from the balance between the induction from the vertical and radial components and the magnetic diffusion that redistribute the field lines.From Equation (9), For conciseness, only the non-ideal terms with second derivative in z−direction are kept since they are the dominating terms in the disk geometry.The two diffusion terms correspond to the vertical diffusion due to pressure gradient and the azimuthal diffusion due to magnetic tension that tends to straighten the field lines.
The non-ideal MHD diffusion in the radial direction is subdominant in the disk geometry (Lee et al. 2021) and is neglected in the following discussions.2016) only discussed disk threaded by an external vertical field.Here we start by discussing the case where the external field is purely horizontal.This implies that B z ≪ B ϕ , B r since the vertical velocity is negligible and no vertical component is generated through induction.
In the regime where field lines are significantly winded, the azimuthal magnetic field is lost through diffusion due to vertical pressure gradient.This corresponds to the weak field case in Lee et al. (2021).Solving from Equation ( 10) and from Equation ( 11), we obtain Since we are already estimating properties at the disk edge without describing the exact profile, a factor 1/2 that would result from the integration along z is neglected for simplicity.The field line straightening due to magnetic tension in the r − ϕ plane has a much larger timescale and is therefore not relevant for the equilibrium.In the scenario where this tension term dominates, a different solution can be found (see Appendix A).

Horizontal vs. vertical field
Recall that the stationary disk size threaded by a vertical field is (Lee et al. 2021) Combining Equations ( 14) and ( 15), one derives the ratio . This ratio is only weakly dependent on the disk parameters.In the stellar mass regime, the disk size is always larger when threaded by a horizontal field.

Disk threaded by inclined field
Without loss of generality, we replace the derivatives with characteristic values in Equations ( 10) and ( 11) to obtain the simplified system, which should still stay valid within correction by some numerical factor: ) Combining the two equations to eliminate B ϕ , the equilibrium radius will satisfy This equation is general to any prescription of density and velocity profiles.Inserting the disk properties from Sect.2.1 leads to It is straightforward to show that the radius reduces to r v or r h when i = 0 or π/2.The results of disk radius as function of i are shown in Figure 1 for several values of r h /r v .There exists a minimal disk size, r min , that is smaller than r v (see derivation in Appendix B).This is particularly evident when r h is close to r v .However, such minor effect might be difficult to detect through observations.The effect of horizontal field component only becomes significant at large inclinations.This might be able to explain that Lee et al.only some, while not many, large disks exist when the field inclination is indeed uniformly distributed.
Combining Equations ( 17) and ( 18) allows to express χ in terms of the disk parameters.The strength of the induced azimuthal field is .
There is a critical value 33 AU, below which B z > B ϕ for i = 0. We recover what was derived in Table 2 of Lee et al. (2021).This is the separation between the magnetic-dominated and inertia-dominated regimes.In the perpendicular configuration i = π/2, lower mass or more strongly magnetized disks can be dominated by the radial field component.These conditions are less difficultly achieved in common disk environments.This is not unexpected given that the magnetic flux from the external field naturally gives rise to both radial and azimuthal components.
The values of χ are plotted against the inclination i for several values of r h /r v in Figure 2. The black lines represent the strength of the vertical (solid) and radial (dashed) field components.The azimuthal field 0.0 0.1 0.2 0.3 0.4 0.5 inclination( ) strength, presented with colored curves, is not only function of r h /r v , but also has to be rescaled according to the actual value of r v .Theses curves move upwards as r v increases, and this is why large disks are always dominated by the toroidal field.
We compared some simulation results from the literature and the difference with our model prediction is always within 30%.For the aligned configuration, Tu et al. ( 2023) presented three cases with varied ambipolar diffusion strength, and we took values from their profiles of B z and ρ at M = 0.25 M ⊙ with their analytical formula for η A .Xu & Kunz (2021) showed profiles of B z and ρ for four snapshots of varying accreted mass.Their ambipolar diffusion coefficient is not displayed thus we used the canonical value η A = 10 19 cm 2 s −1 as reference.Masson et al. (2016) considered two initial values of magnetization.Only total B is shown and we took the saturation value to approximate B z and also used the canonical value for η A .The disk in the weakly magnetized case shows spirals much larger than the model size, while the more regular central part still fits quite well with our model prediction.On the other hand, their two cases misaligned by 40 • , where the derived r/r v ≈ 1.1, also agree within 20% of the model prediction.Hennebelle et al. (2020) presented a series of runs with varied initial magnetization and mass.The cases with 30 • misalignment are well described by the model with aligned field.While the value of B r is not presented, using the canonical value B = 0.1 G also reproduces very well their two cases with 90 • alignment at different levels of initial magnetization.

Disk size population
Our model describes the disk size as function of the inclination angle between the angular momentum and the external magnetic field.Observations of disk size distribution of a population can thus be used to test whether this angle i is randomly distributed or follow some certain distribution P (i).The numerical integration requires some special care since there exists a minimum radius, r min , at inclination i min .By changing the variables, it is possible to perform the integration in i: where i 1 (r) < i min and i 2 (r) > i min are the two values corresponding to the same radius r. 1 In Figure 3, we plot the number fraction of disks with size larger than r, which writes while assuming uniform distribution of inclinations, that is, P (i) ∝ sin(i), or a more aligned distribution, P (i) ∝ 1.As already pointed out in Figure 1, the disk size only approaches r h when the field is significantly inclined.Therefore, the size distribution is dominated by small ones even with a uniform inclination distribution, and worsens with a more aligned distribution.This effect is even more pronounced for higher mass and more magnetized disks (represented by large r h /r v ).This is proba-  bly the reason why large disks are rare from observations (Maury et al. 2022).
For the size distribution of the complete proto-stellar population, we integrate over the IMF and display the results in Figure 4. We relate the disk parameters, including the field strength and η A , to the total mass in order to express everything with one single parameter (see Appendix C).With such parametrization, higher-mass stars have wider distribution of disk sizes.We integrated from 0.05 to 300 M ⊙ , with a lognormal distribution around 0.2 M ⊙ and powerlaw with slope dN/d log M ∝ −1.35 above 1 M ⊙ (Chabrier 2005).To examine the effect of the IMF variation, we also used a slope of −0.75 (Lee & Hennebelle 2018).Three alignment distributions are presented: uniform (P (i) ∝ sin(i), solid curves), moderately aligned (P (i) ∝ 1, dashed curves), and perfectly aligned (dotted curves).The population is seriously dominated by the low-mass stars and thus no significant variation results from different IMF slopes.On the other hand, the disk size distribution seems to be a good proxy to test alignment between the rotation and external magnetic field.Lebreuilly et al. (2023) presented disk size distributions for series radiative non-ideal MHD simulations.Their results generally agree our model prediction for moderately aligned distribution.However, one case with lower initial global magnetization yields larger disks.This might come from the fact that in such conditions the disk magnetization does not follow our model simplification of field saturation regulated by non-ideal MHD (Appendix C).Bate (2018) analyzed the disk size distribution in a clustered star formation stimulation, and   (Chabrier 2005), and orange curves have the IMF slope changed to −0.75.Misalignment angle following uniform distribution, P (i) ∝ sin(i), is plotted with solid lines, and a more aligned distribution, P (i) ∝ 1, is plotted with dashed lines.
they found surprisingly good agreement with observed class 0 disks without introducing magnetic fields.They identified dynamical encounters as an explanation for the compact disk sizes.The size distribution they discovered (their Fig. 28) falls in between our moderately aligned and perfectly aligned cases.These numerical works suggest that there could be multiple factor regulating the embedded disk size.Knowing complete disk parameters including mass, field strength, and alignment angle, rather than reducing to one single parameter, would definitely give more confidence in the prediction of disk sizes.

Role of the Hall effect
In the early phases of disk formation, the Hall effect causes a bimodality in disk sizes, depending on the relative orientations of the angular momentum and the magnetic field (Tsukamoto et al. 2015).However, as already demonstrated analytically in Lee et al. (2021) and numerically in Zhao et al. (2020), this effect is transient and disappears after a few thousands years.Irrespective of the initial angular momentum, a disk with antiparallel angular momentum to the external field always forms at the end.This is due to the conversion of poloidal field in toroidal field, that weakens or strengthens the magnetic braking in the infalling envelope.This is described with an additional term in the induction equation ( 11).When threaded by a horizontal field, the Hall drift of charged particles brings the field lines off the disk plane, which is described with the induction term This generated vertical field is parallel to the rotation in the envelope and anti-parallel inside the disk.The disk and the magnetic field lines will finally adjust to a configuration such that the rotation axis becomes antiparallel to the local magnetic field.In the end, the quasistationary disk size mostly results from the effect of ionneutral (ambipolar) friction.

CONCLUSIONS
We proposed a general formalism for size of a disk threaded by inclined magnetic field, as a result of selfregulation of angular momentum and field strength by the non-ideal MHD effects, notably the ion-neutral friction.While we presented results with some prescriptions for density, velocity, and field strength, this formalism is directly applicable to any other profiles of these physical parameters.Our model allows to infer a disk size distribution from a distribution of the misalignment angle between rotation and magnetic field.We also reduced the disk size expression to function of the mass only.It is then possible to integrate over a stellar mass distribution function to obtain the size distribution of a complete protostellar population.Such model can be tested with size observations of a disk population.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for the thorough reading and constructive comments that helped to better shape the manuscript.The third segment is almost flat.By combining Equations ( C11) and (C13), we can obtain the saturated field strength value The values are consistent with MHD simulations (Masson et al. 2016;Commerçon et al. 2022), and the results of Hennebelle et al. (2016, Figure 1) is essentially recovered with this analytical formalism without need to perform numerical integration.While we leave η A , B 0 , and M and independent parameters in the main text, one can also take into account Equations (C10) and (C15), and apply Equation (C10) with n > n crit when solving Equation ( 19).This leads to . (C17) The density corresponding to this solution is much below the presumed magnetic field saturation density, and is thus not a self-consistent solution.Therefore, instead one should consider Equations (C10) and (C13) when solving Equation (19).By using B ≈ B 0 χ in Equation (C10) with n < n crit , this leads to (C19) Whether the solution is valid depends on the corresponding density.The actual radius is the larger among the two values.This radius is subject to the upper limit imposed by a pure hydrodynamic treatment, which conserves all accreted angular momentum (Lee et al. 2021): where β is the ratio between the rotational and gravitational energy in the core.This hydrodynamic radius could be over-estimated, since some braking should already happen within the envelope, particularly in the horizontal field case.Strictly speaking, there is an inconsistency because the density and the field strength experience an amplification across the shock at the disk edge.However, as shown by Hennebelle et al. (2020, Fig. 3), the amplification of B z is weaker than the amplification of the density.Using Equation (C13) for the sublinear dependence of B z on ρ turns out to be a reasonable estimate, even across the shock, due to its diffusive nature.
Disk threaded by horizontal field: inertia-dominated regime Lee et al. (2021) and Hennebelle et al. (

Figure 1 .
Figure 1.The disk radius, divided by rv plotted against the inclination angle i.Several values of r h /rv are shown, while 2.8 is characteristic of the common disk parameters.

Figure 2 .
Figure 2. Relative strength of the magnetic field components.Colored lines show (r/rv) −5/4 , which are representative of the azimuthal field strength.Solid black line shows cos(i) (z-component) and dashed black line shows sin(i)2/π (r-component).The azimuthal field is dominant in most of the cases.

Figure 4 .
Figure4.Disk size distribution for a complete population following IMF distribution.Blue curves correspond to a population that follows a lognormal around 0.2 M⊙ and a powerlaw of slope −1.35(Chabrier 2005), and orange curves have the IMF slope changed to −0.75.Misalignment angle following uniform distribution, P (i) ∝ sin(i), is plotted with solid lines, and a more aligned distribution, P (i) ∝ 1, is plotted with dashed lines.
YNL acknowledges funding from the National Science and Technology Council, Taiwan (NSTC 112-2636-M-003-001) and the grant for Yushan Young Scholar from the Ministry of Education, Taiwan.PM received financial support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC Starting Grant "Chemtrip", grant agreement No. 949278) BR acknowledges financial support from the National Science and Technology Council, Taiwan (NSTC) under the International Internship Pilot Program (IIPP), and the Kishore Vaigyanik Protsahan Yojana (KVPY) program, India.PH acknowledges funding from the European Research Council synergy grant ECOGAL (Grant: 855130).

Figure 5 .
Figure 5. Magnetic field strength-density relation for collapsing envelope around various central stellar mass, from Equations (C13) and (C14).The gray line shows the relation from Equation (C11) which seperates different regimes of ηA behaviors.The kink is due to the discontinuous expression in Equation (C12), which can be easily removed by applying some smoothing.
4.0 Figure 3. Disk size distribution for varied values of r h /rv.
Misalignment angle following uniform distribution, P (i) ∝ sin(i), is plotted with solid lines, and a more aligned distribution, P (i) ∝ 1, is plotted with dashed lines.