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

Yueh-Ning Lee, Barshan Ray, Pierre Marchand, Patrick Hennebelle

► To cite this version:

Yueh-Ning Lee, Barshan Ray, Pierre Marchand, Patrick Hennebelle. Protoplanetary Disk Size under Nonideal Magnetohydrodynamics: A General Formalism with Inclined Magnetic Field. The Astrophysical Journal Letters, 2024, 961, 10.3847/2041-8213/ad192a. insu-04822485

HAL Id: insu-04822485

https://insu.hal.science/insu-04822485v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Yueh-Ning Lee (李悅寧)¹,²,³ ID, Barshan Ray¹,⁴ ID, Pierre Marchand⁵ ID, and Patrick Hennebelle⁶,⁷,⁸ ID

¹Department of Earth Sciences, National Taiwan Normal University, Taipei 116059, Taiwan; ynlee@ntnu.edu.tw

²Center of Astronomy and Gravitation, National Taiwan Normal University, Taipei 116059, Taiwan

³Physics Division, National Center for Theoretical Sciences, Taipei 106319, Taiwan

⁴Department of Physical Sciences, Indian Institute of Science Education and Research Kolkata, Mohanpur 741246, India

⁵Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UT3-PS, CNRS, CNES, 9 av. du Colonel Roche, F-31028 Toulouse Cedex 4, France

⁶Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS, F-91191 Gif-sur-Yvette, France

⁷IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

⁸LERMA (UMR CNRS 8112), Ecole Normale Supérieure, F-75231 Paris Cedex, France

Received 2023 December 3; revised 2023 December 23; accepted 2023 December 28; published 2024 January 22

## Abstract

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.

Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Magnetohydrodynamics (1964)

## 1. 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 (Zhao et al. 2020; Lee et al. 2021). 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.

## 2. Formalism

### 2.1. Properties of the Envelope and the Disk

The prescriptions follow exactly the same as those in Lee et al. (2021). The density in the collapsing envelope is

proportional to the singular isothermal sphere (Shu 1977)

$$ \rho_e(r) = \frac{\delta_{\rho} C_s^2}{2\pi Gr^2}, \qquad (1) $$

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 right inside the envelope-disk boundary is amplified by the accretion shock, such that

$$ \rho(r) = \frac{\delta_{\rho} C_s^2}{2\pi G r^2} M^2 = \frac{\delta_{\rho} C_s^2}{2\pi G r^2} \left( \frac{u_r(r)}{C_s} \right)^2, \qquad (2) $$

with M being the Mach number. The Keplerian velocity is given by

$$ v_{\text{kep}} \approx \sqrt{\frac{GM}{r}} = \sqrt{\frac{G(M_{*} + M_{\text{d}})}{r}}, \qquad (3) $$

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

$$ h(r) \approx C_s \sqrt{\frac{r^3}{GM}}. \qquad (4) $$

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

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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

$$B_z = B_0 \cos i, \qquad (5)$$

$$B_r = B_0 \frac{2}{\pi} \sin i, \qquad (6)$$

$$B_{\phi} = B_0 \chi, \qquad (7)$$

where $\chi$ is a dimensionless coefficient for the strength of the azimuthal component. The factor $2/\pi$ 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_02rh$ over area $\pi rh$. As a result of symmetry, the circular integral of $B_\phi$ 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.

## 2.2. Nonideal MHD Equations

For treating the disk properties, we present below the nonideal magnetohydrodynamic (MHD) equations. The equation of momentum conservation writes

$$\frac{\partial(\rho u)}{\partial t} = -\nabla \cdot \left[ \rho u \mathbf{u} + \left( P + \frac{B^2}{2} \right) \mathbf{I} - \mathbf{B}\mathbf{B} \right] - \rho \nabla \Phi, \quad (8)$$

where $u$ is the velocity vector, $P$ is the pressure, $\mathbf{B}$ is the magnetic vector, and $\Phi$ is the gravitational potential. The induction equation writes

$$\begin{aligned} \frac{\partial \mathbf{B}}{\partial t} &= \nabla \times (\mathbf{u} \times \mathbf{B}) - \nabla \times \left\{ \eta_A \frac{\mathbf{B}}{\|\mathbf{B}\|} \times (\nabla \times \mathbf{B}) \times \frac{\mathbf{B}}{\|\mathbf{B}\|} \right. \\ & \qquad \left. + \eta_H (\nabla \times \mathbf{B}) \times \frac{\mathbf{B}}{\|\mathbf{B}\|} + \eta_O \nabla \times \mathbf{B} \right\} \\ &= \nabla \times [(\mathbf{u} + \mathbf{u}_A + \mathbf{u}_H) \times \mathbf{B} - \eta_O \nabla \times \mathbf{B}] \\ &= \nabla \times (\mathbf{v} \times \mathbf{B} - \eta_O \nabla \times \mathbf{B}), \end{aligned} \qquad (9)$$

where $\eta_A$, $\eta_H$, and $\eta_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.

## 2.3. Governing Equations of the Disk

Lee et al. (2021) discussed the strong-field and weak-field 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 inertia-dominated 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

$$0 \approx -\frac{u_r}{r} \frac{\partial r u_\phi}{\partial r} + \frac{\partial (B_z B_\phi)}{\rho \partial z} + \frac{\partial (r B_r B_\phi)}{\rho r \partial r}. \qquad (10)$$

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),

$$\begin{aligned} 0 \approx & \frac{\partial u_{\phi} B_z}{\partial z} + \frac{\partial u_{\phi} B_r}{\partial r} \\ & + \frac{\partial}{\partial z} \left[ \frac{\eta_A}{B^2} B_{\phi} \left( B_r \frac{\partial B_r}{\partial z} + B_{\phi} \frac{\partial B_{\phi}}{\partial z} \right) \right] \\ & + \frac{\partial}{\partial z} \left( \frac{\eta_A B_z^2}{B^2} \frac{\partial B_{\phi}}{\partial z} \right). \end{aligned} \qquad (11)$$

For conciseness, only the nonideal 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 nonideal MHD diffusion in the radial direction is subdominant in the disk geometry (Lee et al. 2021) and is neglected in the following discussions.

# 3. Results

## 3.1. Disk Threaded by Horizontal Field: Inertia-dominated Regime

Lee et al. (2021) and Hennebelle et al. (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 \ll B_\phi$, $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

$$\frac{\rho u_{\phi} u_{r}}{r} = \frac{B_{\phi} B_{r}}{r} \qquad (12)$$

from Equation (10) and

$$\frac{u_{\phi} B_{r}}{r} = \frac{\eta_{A} B_{\phi}}{h^{2}} \qquad (13)$$
from Equation (11), we obtain

$$ 
\begin{aligned}
r_h &= \left[ (2\pi)^{-2} G^3 C_s^{-4} \eta_A^2 M^5 \left(B_0 \frac{2}{\pi}\right)^{-4} \right]^{\frac{1}{11}} \\
&= 53.6 \text{ au} \left[ \frac{\eta_A}{10^{19} \text{cm}^2 \text{s}^{-1}} \right]^{\frac{5}{11}} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{5}{11}} \left[ \frac{B_0}{0.1 \text{G}} \right]^{-\frac{4}{11}}.
\end{aligned} 
\quad (14) $$

When simplifying from Equations (11) to (13), it is presumed that $B^2 = B_{\phi}^2 + B_r^2$. 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).

### 3.2. Horizontal versus Vertical Field

Recall that the stationary disk size threaded by a vertical field is (Lee et al. 2021)

$$ 
\begin{aligned}
r_v &= [(2\pi)^{-2} G^1 C_s^0 \eta_A^2 M^3 B_0^{-4}]^{\frac{1}{9}} \\
&= 19.2 \text{ au} \left[ \frac{\eta_D}{10^{19} \text{cm}^2 \text{s}^{-1}} \right]^{\frac{2}{9}} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{1}{3}} \left[ \frac{B_0}{0.1 \text{G}} \right]^{-\frac{4}{9}}.
\end{aligned} 
\quad (15) $$

Combining Equations (14) and (15), one derives the ratio

$$ 
\begin{aligned}
\frac{r_h}{r_v} &= [(2\pi)^4 G^{16} C_s^{-36} \eta_A^{-4} M^{12} B_0^8]^{\frac{1}{99}} \left(\frac{2}{\pi}\right)^{-\frac{4}{11}} \\
&= 2.80 \left[ \frac{\eta_A}{10^{19} \text{cm}^2 \text{s}^{-1}} \right]^{-\frac{4}{99}} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{4}{33}} \left[ \frac{B_0}{0.1 \text{G}} \right]^{\frac{8}{99}}.
\end{aligned} 
\quad (16) $$

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.

### 3.3. 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:

$$ \frac{\rho u_{\phi} u_r}{r} \approx \frac{B_{\phi} B_z}{h} + \frac{B_{\phi} B_r}{r}. \quad (17) $$

$$ 
\begin{aligned}
\frac{u_{\phi} B_z}{h} + \frac{u_{\phi} B_r}{r} &\approx \frac{\eta_A (B_{\phi}^2 + B_r^2) B_{\phi}}{B^2 h^2} \\
&\quad + \frac{\eta_A B_z^2 B_{\phi}}{B^2 h^2} = \frac{\eta_A B_{\phi}}{h^2},
\end{aligned} 
\quad (18) $$

Combining the two equations to eliminate $B_{\phi}$, the equilibrium radius will satisfy

$$ \frac{B_z \sqrt{r}}{\sqrt{\eta_A \rho u_r}} + \frac{B_r h}{\sqrt{\eta_A \rho u_r r}} = 1. \quad (19) $$

This equation is general to any prescription of density and velocity profiles. Inserting the disk properties from Section 2.1

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The disk radius, divided by $r_v$, plotted against the inclination angle $i$. Several values of $r_h/r_v$ are shown, while 2.8 is characteristic of the common disk parameters.</figcaption></figure>

leads to

$$ \left(\frac{r}{r_v}\right)^{\frac{9}{4}} \cos i + \left(\frac{r}{r_h}\right)^{\frac{11}{4}} \sin i = 1. \quad (20) $$

It is straightforward to show that the radius reduces to $r_v$ or $r_h$ when $i=0$ or $\pi/2$.

The results of the disk radius as a 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 the derivation in Appendix B). This is particularly evident when $r_h$ is close to $r_v$. However, such a minor effect might be difficult to detect through observations. The effect of the horizontal field component only becomes significant at large inclinations. This might be able to explain that only some, while not many, large disks exist when the field inclination is indeed uniformly distributed.

Combining Equations (17) and (18) allows to express $\chi$ in terms of the disk parameters. The strength of the induced azimuthal field is

$$ 
\begin{aligned}
\chi &= [r^{-5} (2\pi)^{-2} G^1 C_s^4 \eta_A^{-2} M^3 B_0^{-4}]^{\frac{1}{4}} \\
&= \left(\frac{r}{r_v}\right)^{-\frac{5}{4}} \left(\frac{r_v C_s}{\eta_A}\right)^{-\frac{5}{4}} \left[ \frac{r_v}{33.3au} \right].
\end{aligned} 
\quad (21) $$

There is a critical value 33 au, below which $B_z > B_{\phi}$ 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=\pi/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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relative strength of the magnetic field components. Colored lines show $(r/r_v)^{-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/\pi$ (r-component). The azimuthal field is dominant in most of the cases.</figcaption></figure>

not unexpected given that the magnetic flux from the external field naturally gives rise to both radial and azimuthal components.

The values of $\chi$ 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 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$. These 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 $\rho$ at $M = 0.25 M_\odot$ with their analytical formula for $\eta_A$. Xu & Kunz (2021) showed profiles of $B_z$ and $\rho$ for four snapshots of varying accreted mass. Their ambipolar diffusion coefficient is not displayed thus we used the canonical value $\eta_A = 10^{19} \text{cm}^2 \text{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 $\eta_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^\circ$, where the derived $r/r_v \approx 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^\circ$ 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

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Disk size distribution for varied values of $r_h/r_v$. Misalignment angle following uniform distribution, $P(i) \propto \sin(i)$, is plotted with solid lines, and a more aligned distribution, $P(i) \propto 1$, is plotted with dashed lines.</figcaption></figure>

reproduces very well their two cases with $90^\circ$ alignment at different levels of initial magnetization.

# 4. Discussion

## 4.1. 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$:

$$ 
\begin{aligned}
N(<r) &= \int_{r_{\min}}^{r} \frac{dN}{dr} dr = \int_{r_{\min}}^{r} \frac{dN}{di} \frac{di}{dr} dr \\
&= \int_{r_{\min}}^{r} P(i) \frac{di}{dr} dr = \int_{i_1(\min(r, r_v))}^{i_2(r)} P(i) di,
\end{aligned} 
\quad (22) $$

where $i_1(r) < i_{min}$ and $i_2(r) > i_{min}$ are the two values corresponding to the same radius $r$.<fnref n="9" /> In Figure 3, we plot the number fraction of disks with size larger than $r$, which writes

$$ \text{CDF}(>r) = 1 - \frac{N(<r)}{N(<r_h)}, \quad (23) $$

while assuming uniform distribution of inclinations, that is, $P(i) \propto \sin(i)$, or a more aligned distribution, $P(i) \propto 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$).

<footnote n="9">⁹ We derived such expressions to stay mathematically coherent in the discussions, while $r_{min}$ is usually very close to $r_v$ in the physically relevant regime (see Figure 1). This small difference is easily overwhelmed by uncertainties due to our simplification in some geometrical constants, and should not have strong physical implications.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4. 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 power law of slope -1.35 (Chabrier 2005), and orange curves have the IMF slope changed to -0.75. Misalignment angle following uniform distribution, $P(i) \propto \sin(i)$, is plotted with solid lines, and a more aligned distribution, $P(i) \propto 1$, is plotted with dashed lines.</figcaption></figure>

This is probably the reason why large disks are rare from observations (Maury et al. 2022).

For the size distribution of the complete protostellar population, we integrate over the initial mass function (IMF) and display the results in Figure 4. We relate the disk parameters, including the field strength and $\eta_A$, to the total mass in order to express everything with one single parameter (see Appendix C). With such parameterization, 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 a power law with a slope of $dN/d \log M \propto -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) \propto \sin(i)$, solid curves), moderately aligned ($P(i) \propto 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 nonideal MHD simulations. Their results generally agree with 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 nonideal MHD (Appendix C). Bate (2018) analyzed the disk size distribution in a clustered star formation stimulation, and 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 Figure 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.

## 4.2. 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

$$ \frac{\partial B_{\phi}}{\partial t_{\text{Hall}}} = - \frac{\partial}{\partial z} \left[ \frac{\eta_{\text{H}} B_z}{B} \frac{\partial B_r}{\partial z} \right] \quad (24) $$

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

$$ \frac{\partial B_z}{\partial t_{\text{Hall}}} = - \frac{\partial}{\partial r} \left[ \frac{\eta_{\text{H}} B_r}{B} \frac{\partial r B_{\phi}}{r \partial r} \right]. \quad (25) $$

This generated vertical field is parallel to the rotation in the envelope and antiparallel 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 quasi-stationary disk size mostly results from the effect of ion-neutral (ambipolar) friction.

## 5. Conclusions

We proposed a general formalism for size of a disk threaded by inclined magnetic field, as a result of self-regulation of angular momentum and field strength by the nonideal 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 a 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. Y.N.L. 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. P.M. received financial support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (ERC Starting Grant “Chemtrip,” grant agreement No. 949278). B.R. 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. P.H. acknowledges funding from the European Research Council synergy grant ECOGAL (Grant: 855130).

## Appendix A

### Magnetic-dominated Regime with Horizontal Diffusion (Strong-field Case)

In contrast to Section 3.1, we consider the possible scenario at the very beginning of disk formation, where the toroidal field is not well developed. Under such conditions, the azimuthal component is removed by the field line straightening due to magnetic tension, instead of vertical diffusion. In the absence of vertical field, the only magnetic tension comes from the product between the radial and azimuthal field components. The induction-generated azimuthal field is balanced by horizontal diffusion

$$ \frac{\partial B_{\phi}}{\partial t} \approx \frac{\partial u_{\phi} B_r}{\partial r} + \frac{\partial}{\partial r} \left[ \eta_{\Lambda} \frac{\partial (r B_{\phi})}{\partial r} \right]. \quad (A1) $$

By solving simultaneously Equation (12) and

$$ \frac{u_{\phi} B_r}{r} = \frac{\eta_{\Lambda} B_{\phi}}{r^2}, \quad (A2) $$

we obtain the radius

$$ \begin{aligned} r_h &= \left[ (2\pi)^{-2} G^1 C_s^0 \eta_A^2 M^3 \left(B_0 \frac{2}{\pi}\right)^{-4} \right]^{\frac{1}{9}} \\ &= 23.4 \text{ au} \left[ \frac{\eta_A}{10^{19} \text{cm}^2 \text{s}^{-1}} \right]^2 \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{1}{3}} \left[ \frac{B_0}{0.1 G} \right]^{-\frac{4}{9}}. \end{aligned} \quad (A3) $$

This is almost identical to what was found by Hennebelle et al. (2016), while the only difference comes from the rescaling due to Equation (6). The right-hand side of Equation (A2) comes from the radial diffusion due to gradient of $B_{\phi}$ and azimuthal diffusion due to the $B_{\phi}B_r$ tensor in Equation (9), which are not shown in the simplified expression in Equation (11). In disk geometry, the diffusion timescale, $r^2/\eta_{\Lambda}$, is much larger than that in the inertia-dominated regime, $h^2/\eta_{\Lambda}$. The magnetic tension diffusion is only relevant at the very early stage of prestellar core collapse, where a disk is not yet well developed. This solution is therefore not very relevant for a physical disk.

## Appendix B

### Minimal Disk Size

Taking the derivative of Equation (20) with respect to $i$,

$$ \left[ \frac{9}{4} \left( \frac{r}{r_v} \right)^{\frac{9}{4}} \cos i + \frac{11}{4} \left( \frac{r}{r_h} \right)^{\frac{11}{4}} \sin i \right] \frac{dr}{di} - \left( \frac{r}{r_v} \right)^{\frac{9}{4}} \sin i + \left( \frac{r}{r_h} \right)^{\frac{11}{4}} \cos i = 0, \quad (B1) $$

allows to find the minimum radius at which $dr/di = 0$. The critical radius satisfies the condition

$$ \left(\frac{r_{\min}}{r_v}\right)^{\frac{9}{4}} + \left(\frac{r_{\min}}{r_h}\right)^{\frac{11}{4}} = 1, \quad (B2) $$

and is only slightly inferior to, $r_v$. The corresponding inclination

$$ i_{\min} = \arctan \left[ \left( \frac{r_{\min}}{r_h} \right)^{\frac{11}{4}} \left( \frac{r_{\min}}{r_v} \right)^{-\frac{9}{4}} \right]. \quad (B3) $$

## Appendix C

### Magnetic Field Strength and the Stellar Mass

Due to the ambipolar diffusion, the magnetic field strength saturates at the center of a collapsing core. In order to estimate this value at saturation for cores of various mass, we solve the radial equilibrium between advection and diffusion while assuming the radial component to be negligible:

$$ u_r B_z = \frac{\eta_{\mathrm{AD}}}{B^2} \left( B_z^2 \frac{\partial B_z}{\partial r} + \frac{B_{\phi} B_z}{r} \frac{\partial r B_{\phi}}{\partial r} \right). \quad (C1) $$

For a shallow $B_{\phi}$ profile (shallower than $r^{-1}$), the second term on the right-hand side operates in synergy with the advection term. Therefore, the left-hand side can only be balanced by the gradient of $B_z$ component, whatever the dominant component is. When $B_{\phi} \gg B_z$, $B_{\phi}$ should obey a profile (probably close to $\propto r^{-1}$) such that the inward diffusion due to its gradient is dominated by the outward diffusion due to the gradient of $B_z$ and

$$ u_r B_z = \eta_{\mathrm{AD}} \frac{B_z^2}{B^2} \frac{dB_z}{dr} \quad (C2) $$

always stays true. Replacing the expressions with Equations (1) and (3) for a core that follows the description in Section 2.1, one obtains

$$ \frac{dB_z}{d\rho} = \frac{u_r B^2}{\eta_{\mathrm{AD}} B_z} \frac{dr}{d\rho} = \frac{B^2 \sqrt{GM}}{B_z 2\rho \eta_{\mathrm{AD}}} \left( \frac{C_s^2}{2\pi G \rho} \right)^{\frac{1}{4}}. \quad (C3) $$

The ambipolar diffusion coefficient is dominated by positively charged ions in the regime of our interest and thus

$$ \begin{aligned} \frac{4\pi}{c^2} \eta_A &\approx \sigma_{\perp}^{-1} = \frac{1 + (\omega_M \tau_{Mn})^2}{\sigma_M} \\ &\approx \left\{ \begin{array}{ll} \frac{B^2}{n_M n c^2 m_{H_2} a_{MHe} \langle \sigma_{\mathrm{coll}} w \rangle_{\mathrm{M}}} & \text{for } n < n_{\mathrm{crit}} \\ \frac{n M_{H_2} a_{MHe} \langle \sigma_{\mathrm{coll}} w \rangle_{\mathrm{M}}}{n_M e^2} & \text{for } n > n_{\mathrm{crit}} \end{array} \right. , \end{aligned} \quad (C4) $$

where

$$ n_{\mathrm{crit}} = \frac{e}{c m_{\mathrm{H}_2} a_{\mathrm{MHe}} \langle \sigma_{\mathrm{coll}} w \rangle_{\mathrm{M}}} B = k_{\mathrm{crit}} B \quad (C5) $$

is the density above which $\omega_M \tau_{Mn} < 1$. $\omega_M$ is the cyclotron frequency of the ion, $\tau_{Mn}$ is its characteristic time of collision with the neutral hydrogen, $c$ is the speed of light, $e$ is the elementary charge, and $m_{H_2}$ is the mass of a hydrogen
<figure><img src="image_8.png" /><figcaption>Figure 5. Magnetic field strength–density relation for collapsing envelope around various central stellar mass, from Equations (C7) and (C8). The gray line shows the relation from Equation (C5), which separates different regimes of η_A behaviors. The kink is due to the discontinuous expression in Equation (C6), which can be easily removed by applying some smoothing.</figcaption></figure>

molecule. We adopt the same constants used in Marchand et al. (2016) for positive ions, $a_{MHe} = 1.14$ (Desch & Mouscho vias 2001) is a constant factor, and $\langle \sigma_{coll}w \rangle_M = 9.1 \times 10^{-10}$ is the rate constant for collision with H$_2$ molecules (Pinto & Galli 2008). For the number density of the ions, we follow the fiducial case in (Marchand et al. 2016; Figure 3) and parameterize as

$$ n_M = n_{M0} \left( \frac{n}{n_{\text{ref}}} \right)^\alpha, \quad \text{with} \\ \alpha = \begin{cases} 0.3, & \text{for } n < n_{\text{ref}} \\ 0, & \text{for } n > n_{\text{ref}} \end{cases} \qquad (\text{C6}) $$

where $n_{M0} = 10^{-2} \text{ cm}^{-3}$ and $n_{\text{ref}} = 10^9 \text{ cm}^{-3}$. At high densities, the ion-electron recombination ($\propto n_{H_2}^2$) becomes more efficient than their creation by cosmic rays ($\propto n_{H_2}$), leading to an ion/electron fractional abundance scaling as $n_{H_2}^{-1}$, and therefore a constant abundance (Marchand et al. 2021). The system can be solved for three segments:

$$ B(n) = \left[ \frac{(GM)^{\frac{1}{2}}}{(\alpha + 3/4)} \left( \frac{C_s^2}{2\pi G \mu m_p} \right)^{\frac{1}{4}} \times 4\pi m_{H_2} a_{MHe} \langle \sigma_{coll} w \rangle_M \frac{n_{M0} n^{\alpha+3/4}}{n_{\text{ref}}^{\alpha}} \right]^{\frac{1}{2}}, \quad (\text{C7}) $$

for the two $\alpha$ values, separated at $n_{\text{ref}}$, up to density $n_{\text{crit}}$, and then

$$ B(n) = \frac{n_{\text{crit}}}{k_{\text{crit}}} \exp \left[ \frac{(GM)^{\frac{1}{2}}}{(5/4 - \alpha) \left( \frac{C_s^2}{2\pi G \mu m_p} \right)^{\frac{1}{4}}} \times \frac{2\pi}{c^2} \frac{e^2 n_{M0}}{m_{H_2} a_{MHe} \langle \sigma_{\text{coll}} w \rangle_M n_{\text{ref}}^{\alpha}} \right] (n_{\text{crit}}^{\alpha-5/4} - n_{\text{ref}}^{\alpha-5/4}). \quad (\text{C8}) $$

The expression for density $n < n_{\text{crit}}$ in Equation (C7) is always true, irrelevant to the strength of $B_φ$, since the dependence of $\eta_A$ cancels out with the total field strength. However, Equation (C7) is only true when negligible rotation is present and the vertical field component dominates throughout the whole domain. In the

presence of significant toroidal field, vertical equilibrium has to be solved simultaneously. If $B_φ$ saturates at some value, a very shallow profile decreasing toward high densities will be recovered. We plot the $B-\rho$ relations for several values of core mass in Figure 5. This $B-\rho$ relation is derived for the collapsing envelope, and is therefore also valid for the scenario with field lines perpendicular to the rotation axis.

The third segment is almost flat. By combining Equations (C5) and (C7), we can obtain the saturated field strength value

$$ B_{\text{sat}} \approx \left[ \frac{(GM)^{\frac{1}{2}}}{(\alpha + 3/4)} \left( \frac{C_s^2}{2\pi G \mu m_p} \right)^{\frac{1}{4}} \frac{4\pi e k_{\text{crit}}^{\alpha - \frac{1}{4}} n_{\text{M0}}}{c n_{\text{ref}}^{\alpha}} \right]^{\frac{1}{\frac{5}{4} - \alpha}} \\ = 5.5 \times 10^{-2} G \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{2}{5}}. \qquad (\text{C9}) $$

The values are consistent with MHD simulations (Masson et al. 2016; Commercon 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 $\eta_A$, $B_0$, and $M$ and independent parameters in the main text, one can also take into account Equations (C4) and (C9), and apply Equation (C4) with $n > n_{\text{crit}}$ when solving Equation (19). This leads to

$$ \left(\frac{r}{r_v}\right)^{\frac{15}{4}} \cos i + \left(\frac{r}{r_h}\right)^{\frac{17}{4}} \sin i = 1, \qquad (\text{C10}) $$

where

$$ r_v = 26.8 \text{ au} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{17}{75}} \quad \text{and} \quad r_h = 50.2 \text{ au} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{27}{85}}. \qquad (\text{C11}) $$

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 (C4) and (C7) when solving Equation (19). By using $B \approx B_0 \chi$ in Equation (C4) with $n < n_{\text{crit}}$, this leads to

$$ \left(\frac{r}{r_v}\right)^{1-\frac{9}{4}\alpha} \cos i + \left(\frac{r}{r_h}\right)^{\frac{3}{2}-\frac{9}{4}\alpha} \sin i = 1, \qquad (\text{C12}) $$

where

$$ r_v = \begin{cases} 4.50 \text{ au} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{1}{13}} & \text{and } r_h = \begin{cases} 334 \text{ au} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{21}{33}} \\ 199 \text{ au} \left[ \frac{M}{0.1 M_{\odot}} \right]^{\frac{1}{2}} \end{cases}. & (\text{C13}) \end{cases} $$

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):

$$ r_{\text{hydro}} = 111 \text{ au} \left[ \frac{\beta}{0.02} \right]^{\frac{1}{2}} \left[ \frac{M}{0.1 M_{\odot}} \right], \qquad (\text{C14}) $$
where $\beta$ is the ratio between the rotational and gravitational energy in the core. This hydrodynamic radius could be overestimated, 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; Figure 3), the amplification of $B_z$ is weaker than the amplification of the density. Using Equation (C7) for the sublinear dependence of $B_z$ on $\rho$ turns out to be a reasonable estimate, even across the shock, due to its diffusive nature.

## ORCID iDs

Yueh-Ning Lee (李悦寧) https://orcid.org/0000-0003-3497-2329

Barshan Ray https://orcid.org/0009-0003-6115-1419

Pierre Marchand https://orcid.org/0000-0002-4577-8292

Patrick Hennebelle https://orcid.org/0000-0002-0472-7202

## References

* Bate, M. R. 2018, MNRAS, 475, 5618

* Chabrier, G. 2005, in Astrophysics and Space Science Library, Vol. 327, The Initial Mass Function 50 Years Later, ed. E. Corbelli, F. Palla, & H. Zinnecker (Dordrecht: Spinger), 41

* Commerçon, B., González, M., Mignon-Risse, R., Hennebelle, P., & Vaytet, N. 2022, A&A, 658, A52

* Desch, S. J., & Mouschovias, T. C. 2001, ApJ, 550, 314

* Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJL, 830, L8

* Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

* Kuznetsova, A., Hartmann, L., & Heitsch, F. 2020, ApJ, 893, 73

* Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2023, 19672, arXiv:2310.19672

* Lee, Y.-N., & Hennebelle, P. 2018, A&A, 611, A88

* Lee, Y.-N., Marchand, P., Liu, Y.-H., & Hennebelle, P. 2021, ApJ, 922, 36

* Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50

* Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

* Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32

* Maury, A., Hennebelle, P., & Girart, J. M. 2022, FrASS, 9, 949223

* Pinto, C., & Galli, D. 2008, A&A, 484, 17

* Shu, F. H. 1977, ApJ, 214, 488

* Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, ApJL, 810, L26

* Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, in Astronomical Society of the Pacific Conf. Ser. 534, ed. Y. Inutsuka, A. Maury, B. Commercon et al. (San Francisco, CA: ASP), 317

* Tu, Y., Li, Z.-Y., Lam, K. H., Tomida, K., & Hsu, C.-Y. 2023, MNRAS, 527, 10131

* Xu, W., & Kunz, M. W. 2021, MNRAS, 508, 2142

* Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, SSRv, 216, 43