© The Authors 2025

1 Introduction

During its formation, a low mass star accretes through a flattened centrifugally supported circumstellar disk. Such a disk inevitably forms due to the angular momentum contained in the parent molecular cloud or through local subsonic turbulence (Yorke et al. 1993; Li et al. 2014; Verliat et al. 2020). For the accretion of gas to proceed through the disk, the accreting material must efficiently lose its angular momentum. Gravitational instability (GI) facilitates angular momentum redistribution during protostellar disk formation from a magnetized prestellar core (Tsukamoto et al. 2015; Masson et al. 2016; Xu & Kunz 2021; Mauxion et al. 2024) and continues to play a dominant role during the early embedded phases of the disk’s evolution (Vorobyov & Basu 2006, 2009; Kratter & Lodato 2016; Rice 2016). This occurs when the disk is sufficiently massive as compared to the host star, and the self-gravity of the disk leads to the formation of spiral arms, density fluctuations, and fragmentation, resulting in efficient global mass transport.

Later in the Class II stage, a protoplanetary disk (PPD) is canonically thought to evolve viscously, with magnetorotational instability (MRI) providing the turbulence necessary for this viscosity (Shakura & Sunyaev 1973; Pringle 1981; Balbus & Hawley 1991; Turner et al. 2014). When a weak magnetic field embedded in the disk interacts with the partially ionized material, the Keplerian shear in the disk gives rise to MRI. The action of MRI requires sufficiently ionized gas, and cosmic rays are thought to be the primary source of ionization in PPDs (Umebayashi & Nakano 1981, 1988; Thi et al. 2019). The cosmic rays penetrate the disk surface from the outside, passing through the entire disk thickness at large radii. The innermost regions have a high enough surface density that the cosmic rays cannot reach the midplane, resulting in insufficient ionization to couple to the magnetic field. This suppresses the MRI and causes the PPD to accrete through a layered structure, wherein the accretion occurs only through the ionized surface layers, and a magnetically dead zone is formed at the midplane (Gammie 1996). The dead zone forms a bottleneck for mass and angular momentum transport and gives rise to such disk substructures as rings and vortices, which form pressure maxima and sites for dust growth (Dzyurkevich et al. 2010; Pinilla et al. 2012; Kadam et al. 2022; Regály et al. 2023). Spatially resolved images of PPDs indeed show an abundance of disk substructure in the majority of large, bright disks (Andrews et al. 2009; Flock et al. 2015; Andrews et al. 2018; Zhang et al. 2018; Law et al. 2021; Hsieh et al. 2024).

However, several lines of reasoning point toward a dramatically different picture of PPD evolution, at least in the Class II stage. The Shakura-Sunyaev α-parameter quantifies the efficiency of angular momentum transport in the standard theory of viscous accretion (Shakura & Sunyaev 1973). In the absence of a dead zone, accretion within the disk is most efficient, and a typical value of α is about 0.01 (Hartmann et al. 1998; Hughes et al. 2011; Andrews et al. 2010). However, recent observations of protoplanetary disks indicate that the inferred value of the α parameter is much lower, which is insufficient for the disk to accrete within its typical observed lifespan of a few million years (Williams & Cieza 2011). For example, the dust rings observed within PPDs by the Atacama Large Millimeter/ Submillimeter Array (ALMA) are significantly thin, and studies of dust properties suggest the turbulence to be α ~ 10⁻³–10⁻⁴ (Pinte et al. 2016; Dullemond et al. 2018; Muro-Arena et al. 2018; Franceschi et al. 2023). However, some uncertainty still exists from the assumption of dust properties since the dusty rings can only constrain the ratio of α to the Stokes number. Spatially and spectrally resolved molecular line measurements also indicate that the contribution of turbulence to the line broadening is consistent with weak turbulence (Flaherty et al. 2015, 2017, 2020). Conversely, from a theoretical viewpoint, the disks in recent magnetohydrodynamic (MHD) simulations also point toward marginal excitation of MRI turbulence. When the nonideal MHD effects are included, especially ambipolar diffusion, the simulations show an increasingly laminar flow and much less overall turbulence. Instead, magnetic disk winds are launched near the disk surface, which carry a significant amount of angular momentum (Suzuki & Inutsuka 2009; Bai & Stone 2013; Fromang et al. 2013; Simon et al. 2013; Béthune et al. 2017; Zhu & Stone 2018; Cui & Bai 2020; Lesur 2021b).

The mechanism for angular momentum loss via disk winds was first described by Blandford & Payne (1982), and it is similar to magnetic braking of rapidly rotating stars. When the disk magnetic field makes a sufficiently large angle with the normal (>30°), the outflowing gas, which is constrained to move along the open magnetic field lines, is accelerated centrifugally. The gas flow eventually bends the poloidal magnetic field component into a toroidal field, and the gradient of the toroidal magnetic pressure assists in the collimation of the wind (Pudritz & Ray 2019). We note that these magnetocentrifugal winds are fundamentally different from photoevaporative winds; the latter are launched thermally from the disk surface due to external irradiation and do not exert a torque on the disk (Hollenbach et al. 1994; Alexander et al. 2014; Ercolano & Pascucci 2017). Photoevaporation is negligible in the early disk evolution, although it becomes significant in the late stages, when the stellar accretion rate becomes comparable to the photoevaporation rate. Photoevap- orative winds are considered to be responsible for the eventual dispersal of the gas in PPDs (Alexander et al. 2006, 2014). In reality, these two mass ejection processes operate simultaneously and can be collectively considered as a “magneto-thermal wind" (Bai et al. 2016; Lesur et al. 2023).

High spectral resolution observations of gas lines that trace unbound disk material (forbidden lines such as [OI] λ6300 and [SII] λ6731) show the ubiquitous presence of outflows around PPDs (Pascucci et al. 2023). These outflows are divided into a high-velocity component (HVC) traveling at hundreds of kilometers per second and a low velocity component (LVC) of 10–50 km s⁻¹ (Natta et al. 2014; McGinnis et al. 2018; Nisini et al. 2018; Rab et al. 2022; Fang et al. 2023). The HVC is interpreted as a collimated axial jet, which is launched very close to the surface of the accreting object via magnetic processes. The LVC component originating beyond the dust evaporation radius and spanning a majority of the disk is typically identified as magnetocentrifugal wind, although it may also be interpreted as X-ray driven photoevaporative flow (Rab et al. 2023). Spatially resolved outflows of young objects show small-scale rotating flows as slow molecular winds in addition to fast moving jets (Güdel et al. 2018; Lee et al. 2021). These winds show a nested onion-like morphology at large opening semi-angles (20–50°) and may originate within the inner ten au of the disk (de Valon et al. 2020; Pascucci et al. 2025). The mass flux estimates of the disk winds are exceedingly significant, with one to ten times the jet mass flux and 0.3 to three times the accretion rate onto the central protostar (Tabone et al. 2020). Thus, an emerging paradigm of protoplanetary disk evolution suggests that in addition to redistribution of angular momentum through gravito- viscous torques, extraction of mass and angular momentum via magnetic disk winds plays a crucial role as a driver of accretion.

As the winds fundamentally arise from the interaction of partially ionized gas with the magnetic field, theoretical insights into disk winds primarily come from 3D non-ideal MHD simulations. This includes both local shearing box studies (Bai & Stone 2013; Simon et al. 2013; Suzuki & Inutsuka 2009) and global simulations, which resolve a significant fraction of the disk (Béthune et al. 2017; Zhu & Stone 2018; Cui & Bai 2020; Lesur 2021b). Since long integration times of the order of the disk lifetime (a few million years) are not feasible with 3D simulations, disk winds can be studied with simplified semi-analytical global evolutionary models in 1D (Kunitomo et al. 2020; Alessi & Pudritz 2022; Somigliana et al. 2022; Tabone et al. 2022). However, certain important processes occurring in PPDs cannot be faithfully represented in 1D axisymmetry. For example, in the early stages of disk evolution, GI leads to the formation of large-scale spirals and clumps, which may eventually form wide orbit gas giants (Boss 1997; Vorobyov & Basu 2010a; Kratter & Lodato 2016). In the later stages of evolution, disk-planet interactions are also intrinsically axially asymmetric (Paardekooper et al. 2023). Thus, theoretical studies of long-term accretion disk evolution are often performed in the plane of the disk (r, ϕ) in 2D, with the thin-disk approximation. The vertical hydrostatic equilibrium is imposed throughout the disk so that the vertical temperature gradients and motions are neglected. Such an approximation typically works very well for PPDs since the ratio of gas scale height to the local radius (H_g/ r) is much less than unity. This approach is powerful for studying the long-term behavior of PPDs because it relaxes many of the simplifying assumptions of the 1D Pringle (1981)-like viscous equation of disk evolution while also being computationally feasible.

Almost all of the thin-disk models of protoplanetary disks to date assume a viscous evolution of the disk. The disk selfgravity is often considered, but magnetic winds are essentially neglected. Recently, Kimmig et al. (2020) and Elbakyan et al. (2022) have considered the effects of wind-driven accretion in the thin-disk limit in the context of planetary migration and gap opening, respectively. Nevertheless, these models leave much room for improvement. The wind parameters, such as its escape velocity and the angle of the magnetic field, cannot be held constant over the disk lifetime. Similarly, the vertical angular momentum transport of magnetic disk winds is fundamentally different from the Shakura & Sunyaev (1973) α parameter and cannot be modeled as an additional viscosity, as the latter only redistributes the angular momentum within the disk. In addition, the co-evolution of the dust component is not explicitly considered. Furthermore, accretion variability during star formation complicates the situation, especially the powerful FUor-type outbursts, wherein the accretion luminosity of the protostar increases by hundreds of times (Hartmann & Kenyon 1996; Audard et al. 2014; Fischer et al. 2023). These events are important over the long term because they can significantly alter the disk composition, dynamics, and snow lines, while as much as 10–35% of the main-sequence stellar mass can be accreted via episodic accretion (Offner & McKee 2011; Dunham & Vorobyov 2012; Fischer et al. 2019; Molyarova et al. 2018; Vorobyov et al. 2020a, 2022). The accretion variability during these events should also directly lead to a corresponding variability in the wind outflow rate (Konigl & Pudritz 2000; Ellerbroek et al. 2014). With these issues in mind, a model framework of magnetic disk winds for protoplanetary disks in the thin-disk limit, which incorporates all of the aforementioned phenomena, is desirable.

In this work, we propose a new global model of magnetic disk winds that takes into account the insights gained from local shearing box MHD simulations. The results of the shearing box simulations conducted at several radial distances can be stitched together in order to estimate the wind behavior throughout the PPD (Bai 2013; Bai & Stone 2013). We improve upon the wind mass loss rate and surface stress estimates of Bai et al. by taking into account the corrections that are essential for long-term disk evolution. The resulting prescription of winds is incorporated self-consistently into the global MHD code FEOSAD, which is capable of simulating the formation and evolution of PPDs over long timescales in the thin-disk limit (Vorobyov et al. 2020b). Thus, all three major sources of angular momentum transport are now included, namely, turbulent viscosity, gravitational instability, and magnetic disk winds. We follow the disk evolution from the earliest embedded Class 0 or I stages through to the later Class II stage, demonstrating that the proposed model of magnetic disk winds is a reasonable approximation in the thin- disk limit. We confirm that the disk winds have a significant impact on the disk evolution and global disk properties. Most notably, the disks with winds tend to be smaller in both size and mass. The simulation results are post-processed with the radiation thermo-chemical disk code PRODIMO (Woitke et al. 2009, 2016), and the obtained synthetic observations are compared with ALMA survey data of Class II disks. Our analysis demonstrates that such a comparison can be used to guide and constrain the disk wind model parameters.

In the interest of the conciseness and readability of the study, we present in this paper our disk wind model in detail and explain the steps of obtaining synthetic observations and the subsequent comparison with large-scale surveys. In a second companion paper, we will comprehensively investigate the properties of magnetic disk winds, for example, characteristic bulk and radial properties of the winds, their effects on substructure formation, and comparison of different modes of angular momentum transport. The structure of the current paper is as follows. Section 2.1 describes our MHD disk evolution code, and Sect. 2.2 elaborates on the magnetic wind model. In Sect. 2.3, we explain how the synthetic observations are obtained, while Sect. 2.4 describes the initial conditions. In Sect. 3, the results of the numerical simulations as well as radiation thermo-chemical modeling are presented in the context of observational constraints and theoretical consensus. The findings are summarized in Sect. 4.

2 Methods

2.1 Disk evolution with FEOSAD

In this section, we summarize the underlying MHD model that is used for the evolution of the PPD in the thin-disk limit. The magnetic disk winds described in Sect. 2.2 are added to this base model in a modular fashion. The hydrodynamic simulations are conducted using the code Formation and Evolution Of a Star And its circumstellar Disk (FEOSAD, Vorobyov et al. 2018a, 2020b). The numerical implementation of hydrodynamic algorithms is based on the operator-split procedure similar in methodology to the ZEUS code (Stone & Norman 1992). We take into account the self-gravity of the disk, co-evolution of the coupled two-part dust component, calculation of magnetic fields in the flux-freezing approximation, and a dead zone modeled with an adaptive-α formulation, where the effective α is derived from ionization-recombination balance equations. Here we describe the minimum set of equations, while a more detailed description of the model can be obtained in Vorobyov et al. (2020b) and references therein.

The equations of mass continuity, conservation of momentum, and energy transport for the gas component in the thin disk limit can be described as follows: $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}+\nabla \cdot \left({\Sigma }_{\text{g}}v\right)=-{\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}g}},$$(1) $$\begin{array}{l}\frac{\partial }{\partial t}\left({\Sigma }_{\text{g}}v\right)+\nabla \cdot \left({\Sigma }_{\text{g}}v\otimes v\right)=-\nabla 𝒫+{\Sigma }_{\text{g}}\left(𝑔+{𝑔}_{*}\right)+\nabla \cdot \Pi \\ -{\Sigma }_{\text{d},\text{gr}}{f}_{\text{b}}+\frac{B_{ɀ}{B}_p^{+}}{2\pi }-H_{\text{g}}\nabla \left(\frac{B_{ɀ}^2}{4\pi }\right)-T^{\text{Max}},\end{array}$$(2) $$\frac{\partial e}{\partial t}+\nabla \cdot (ev)=-𝒫(\nabla \cdot v)-\mathrm{\Lambda }+\Gamma +\nabla v:\Pi -{\dot{e}}_{\text{w}}.$$(3)

We note that the equations are solved in the plane of the disk in polar coordinates (r, ϕ) so that the vectors and tensors have only the planar components. Here Σ_g is the gas surface density, v is the gas velocity in the disk plane and ${\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}g}}$ is the rate of mass loss (per unit disk area) from magnetic winds. The vertically integrated gas pressure, P, is calculated via the ideal equation of state: P = (γ − 1)e, where e is the internal energy per surface area and the assumed value of γ is 7/5. The term g_* is the gravitational acceleration by the central protostar, while the gravitational acceleration in the disk plane, g, takes into account self-gravity of both the gaseous and the dusty components by solving the Poisson equation (see Binney & Tremaine 2008). Thus, the gravitational torques, which are a major source of mass and angular momentum transport in the early stages of disk evolution, are considered self-consistently. The turbulent viscosity is taken into account via the viscous stress tensor Π (Vorobyov & Basu 2010b), while the Maxwell’s stress tensor, T^Max originates in the disk winds. The only component relevant for the surface stress is $T_{ɀ\varphi }^{\text{Max}}$ , which forms a sink term in the ϕ component of the momentum equation. The term f_b is the drag force per unit mass due to the back-reaction of dust on to the gas, and Σ_d,gr is the surface density of grown dust. Later in this section, the same back-reaction term enters the dynamics of grown dust. The vertical scale height of the gas disk, H_g = c_s/Ω, is calculated assuming vertical hydrostatic balance. The magnetic field in the disk has only a vertical component B_ɀ, while planar components exist at the top and bottom surfaces. The planar component of the magnetic field at the top and bottom surfaces of the disk are denoted by ${\mathit{B}}_p^{+}$ and ${\mathit{B}}_p^{-}$, and a symmetry about the disk midplane is assumed, implying ${\mathit{B}}_p^{-}=-B_p^{+}$. This approach of evolving the magnetic field is similar to evolution of a thin magnetized disk (Ciolek & Mouschovias 1993; Basu & Mouschovias 1994; Vorobyov & Basu 2006).

The energy transport equation (Eq. (3)) pertains to the evolution of the internal energy per surface area, e. The disk midplane temperature (T_mp), which is assumed equal for the dust and gas components, is obtained by solving this energy balance equation. The heating (Γ) and cooling (Λ) rates are based on the analytical solution of the radiation transfer equations in the vertical direction (Dong et al. 2016). The heating function takes into account the irradiation of the disk surface from the stellar as well as background irradiation and radiation cooling by dust grains is also considered (Vorobyov et al. 2018a). The irradiation from the central star is calculated using incidence angle of radiation arriving at the disk surface, which in turn is obtained using the local scale height and a fixed flaring angle (Vorobyov & Basu 2010b; Dong et al. 2016). These calculations neglect the self-shielding of the disk, which may arise in the innermost regions due to formation of concentric rings. However, this effect is considered to be insignificant, as viscous heating dominates the energy balance in the midplane of this inner disk region. The resulting model has a flared structure, wherein the disk vertical scale height increases with radial distance. Both the disk and the infalling envelope receive a fraction of the irradiation energy from the central protostar, both in terms of accretion and photospheric luminosity (Vorobyov & Basu 2009). The optical depths in the calculations are proportional to the total dust surface density; however, they do not take into account the dust growth. The Planck and Rosse- land mean opacities needed for these calculations are taken from Semenov et al. (2003). The term involving viscous stress tensor accounts for the viscous heating, which typically dominates the inner disk and becomes significant for an extended region during outbursts (Vorobyov & Basu 2010b; Kadam et al. 2020). The last term in Eq. (3) corresponds to the internal energy of matter carried away by the magnetic winds.

The dust is modeled with two components and evolves via the processes of coagulation, fragmentation, and drift (Vorobyov et al. 2018a). The small dust component is assumed to be fully coupled with the gas, while the grown dust drifts with respect to the gas and contributes to the back-reaction term in Eq. (2). The small dust has a grain size of a_min < a < a_*, and grown dust corresponds to a size of a_* ≤ a < a_max, where a_min = 5 × 10⁻³ µm and a_* = 1 µm. Here, a_max is a dynamically varying maximum size of the dust grains, which depends on the efficiency of radial dust drift and the rate of dust growth. The dust grains are assumed to have a constant density of ρ_s = 3.0g cm⁻³. The equations of continuity and momentum conservation for small and grown dust components are $$\frac{\partial {\Sigma }_{\text{d},\text{sm}}}{\partial t}+\nabla \cdot \left({\Sigma }_{\text{d},\text{sm}}v\right)=-S\left(a_{\max }\right)-{\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}d}},$$(4) $$\frac{\partial {\Sigma }_{\text{d},\text{gr}}}{\partial t}+\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}u\right)=S\left(a_{\max }\right),$$(5) $$\begin{array}{l}\frac{\partial }{\partial t}\left({\Sigma }_{\text{d},\text{gr}}u\right)+\nabla \cdot \left({\Sigma }_{\text{d},\text{gr}}u\otimes u\right)={\Sigma }_{\text{d},\text{gr}}\left(\mathit{𝑔}+{\mathit{𝑔}}_{*}\right)\\ +{\Sigma }_{\text{d},\text{gr}}{f}_{\text{b}}+S\left(a_{\text{max}}\right)v,\end{array}$$(6)

where Σ_d,sm and Σ_d,gr are the surface densities of small and grown dust, respectively. The term u denotes the planar components of the grown dust velocity and ${\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}d}}$ accounts for the small dust blown away with the disk winds. Since the winds originate in the upper layers of the disk surface, we assume that they do not carry grown dust. S is the rate of conversion from small to grown dust per unit surface area, which is a function of the maximum size of the dust (a_max). The dust is assumed to mix vertically with the gas, which is a reasonable approximation for a young disk evolving under gravitational and viscous torques. S (a_max) is derived from the assumption that the size distributions of both the dust populations can be described by a power law with an exponent of −3.5. The discontinuity in the dust size distribution at a_* is assumed to get smoothed out, which implies the dominant role of dust growth as compared to the dust flow (Molyarova et al. 2021). The evolution of the maximum size of the grown dust can be expressed as an advection equation modified for the presence of a source term: $$\frac{\partial a_{\max }}{\partial t}+(u\cdot \nabla )a_{\max }=\mathcal{D}.$$(7)

Here, D = ρ_dv_rel/ρ_s is the growth rate that accounts for coagulation and ρ_d is the total dust volume density. The relative dust-to-dust velocity, v_rel, is calculated by considering the main sources of relative velocities – the Brownian motion and turbulence (Vorobyov et al. 2018a). The maximum size that the dust can achieve is limited by the fragmentation barrier $$a_{\text{frag\hspace{0.17em}}}=\frac{2{\Sigma }_{\text{g}}{v}_{\text{frag\hspace{0.17em}}}^2}{3\pi {\rho }_{\text{s}}{\alpha }_{\text{turb\hspace{0.17em}}}{c}_{\text{s}}^2},$$(8)

where v_frag is the fragmentation velocity, c_s is the sound speed and α_turb corresponds to the turbulence parameter (Birnstiel et al. 2012), wherein the possible contribution from gravitotur- bulence is not considered (see, e.g., Vorobyov et al. 2023). We set the fragmentation velocity to a canonical value of 3 m s⁻¹, although there is significant ambiguity due to sticking properties of different grain compositions, for example, bare silicate versus icy grains (Blum & Wurm 2008; Wada et al. 2009; Steinpilz et al. 2019). When a_max exceeds a_frag , the growth rate D is set to zero and a_max is set equal to a_frag , which guarantees that the particle size is capped at a_frag . We note that in our approach, the dust growth is limited to keep the size of dust particles within the Epstein regime so that the back-reaction term is expressed as $$f_{\text{b}}=\frac{v-u}{t_{\text{stop}}},$$(9)

where the stopping time is t_stop = a_maxρ_s/c_sρ_g. The update of gas and dust velocities due to the friction force is calculated using a fully implicit scheme (Stoyanovskaya et al. 2018, 2020).

The viscosity in a protoplanetary disk is thought to originate primarily due to the turbulence generated by MRI (Balbus & Hawley 1991; Turner et al. 2014). As the cosmic rays externally penetrate the disk surface, the inner disk (≲10 au) accretes through the MRI-active surface layers, while a magnetically dead zone is formed at the disk midplane (Gammie 1996; Simon et al. 2013). Although there may be quantitative uncertainty, this picture of layered accretion is considered generally valid, even in the presence of disk winds (Armitage & Kley 2019). We modeled the magnetically dead zone using an effective and adaptive Shakura & Sunyaev (1973) α parameter (Bae et al. 2014; Kadam et al. 2019): $${\alpha }_{\text{turb}}=\frac{{\Sigma }_{\text{MRI}}{\alpha }_{\text{MRI}}+{\Sigma }_{\text{dz}}{\alpha }_{\text{dz}}}{{\Sigma }_{\text{MRI}}+{\Sigma }_{\text{dz}}},$$(10)

where Σ_MRI is the gas surface density of the MRI-active layer and Σ_dz is that of the magnetically dead layer so that Σ_g = Σ_MRI + Σ_dz. Here, α_MRI and α_dz correspond to the strength of the turbulence in the MRI-active layer and the dead zone, respectively. The Σ_MRI is obtained from the ionization fraction (x), which in turn is determined by solving the ionization balance equation $$(1-x)\xi ={\alpha }_{\text{r}}{x}^2{n}_{\text{n}}+{\alpha }_{\text{d}}x{n}_{\text{n}}$$(11)

where ξ is the ionization rate, α_r is radiative recombination rate, n_n is the number density of neutrals, and α_d is the total rate of recombination on to the dust grains (Dudorov & Khaibrakhmanov 2014; Dudorov & Sazonov 1987; Balduin et al. 2023). We consider the ionization by cosmic rays and radionuclides; however, effects of stellar far-ultraviolet and X-ray radiation are neglected since their contribution is typically significantly smaller in the inner disk and limited to the upper disk atmosphere, which is not resolved in our thin-disk model (Igea & Glassgold 1999; Bergin et al. 2007; Perez-Becker & Chiang 2011).

In the simulations, α_MRI is set to the canonical value of 0.01, while lower values of 10⁻³ and 10⁻⁴ are also explored, the latter being consistent with recent observations of low viscosity in PPDs (Flaherty et al. 2015; Dullemond et al. 2018; Muro-Arena et al. 2018; Flaherty et al. 2020; Franceschi et al. 2023). 3D MHD simulations of FU Orionis type outbursts suggest that the α values in the inner disk region directly involved in the outburst can be much larger than that in the outer disk (e.g., Zhu et al. 2020). During the MRI-type outbursts, the inner disk gets significantly hot and exceeds ionization temperature for alkali metals of ≈1000 K (Umebayashi & Nakano 1981, 1988). The local thermal ionization then increases significantly and dominates the total ionization fraction. Thus, we set α_MRI to the maximum value of 0.1, when the local thermal ionization fraction exceeds 10⁻¹⁰. The larger value of the α parameter results in a better fit to the duration of FUor-type outbursts that is inferred from observations (Zhu et al. 2009; Vorobyov et al. 2020b). The ionization threshold value is larger than that required for MRI activation in the midplane; however, it is justified since the upper regions of the disk are significantly more ionized from external irradiation (Bai 2011; Desch & Turner 2015). The method of transitioning αMRI as a function of ionization fraction is a better approximation of disk physics, as compared to the previous criterion in Vorobyov et al. (2020b), where this switch was made for the entire innermost five au region. Inside the dead zone, a small value of residual viscosity corresponding to α_dz = 10⁻⁵ is considered, which may arise from purely hydrodynamic turbulence driven by the Maxwell stress in the active layer (Fleming & Stone 2003; Okuzumi & Hirose 2011).

The co-evolution of the magnetic field is implemented in the flux-freezing limit, while the non-ideal MHD effects (Balbus & Terquem 2001; Kunz & Balbus 2004) are neglected due to prohibitively expensive computational costs. However, while calculating the thickness of the MRI active surface layer above the midplane (Σ_crit), we take into account Ohmic diffusivity, which dominate the inner disk region. Equating the wavelength of the most unstable MRI mode to the gas scale height of the disk yields $${\Sigma }_{\text{crit\hspace{0.17em}}}={\left[{\left(\frac{\pi }{2}\right)}^{1/4}\frac{c^2{m}_{\text{e}}{\langle \sigma v\rangle }_{\text{en}}}{e^2}\right]}^{-2}{B}_{ɀ}^2{H}_{\text{g}}^3{x}^2,$$(12)

where e is the charge of an electron, m_e is the mass of an electron and 〈σv〉_en is the slowing-down coefficient (see Vorobyov et al. 2020b, for more details). We note that 〈σv〉_en is now updated to the empirically verified 10⁻⁹ cm3 s⁻¹ (Draine 2011; Das & Basu 2021), instead of the previously used value of 10⁻⁷ cm3 s⁻¹ that was based on a classical estimate (Nakano 1984). The dead zone is present only when Σ_g > 2Σ_crit, in which case Σ_MRI is set to 2Σ_crit. If Σ_g ≤ 2Σ_crit, the disk is considered to be fully MRI active and there is no dead zone. The MRI-active surface layer, as calculated with Eq. (12), depends sensitively on the ionization balance, thermal balance and the magnetic field evolution in the disk. However, it results in a low-viscosity dead zone in the approximately innermost ten au of the disk (Eq. (10)), as anticipated by simulations (Gammie 1996; Dzyurkevich et al. 2010; Vorobyov et al. 2020b). The evolution of B_z is computed by solving the advection equation and the planar components of magnetic field, B_p , is evaluated as the solution of a magnetic analogue of the Poisson integral (Vorobyov et al. 2020b). We note that the assumption of flux-freezing limit leads to significantly overestimated B_z , which has consequences for the disk wind model.

2.2 Magnetic disk wind model

The magnetic disk wind model that we adopt is primarily based on the insights gained from the local MHD shearing box simulations of Bai (2013), which we modify for long-term disk evolution. The results of a series of shearing box simulations conducted at different radii can be stitched together to obtain a global picture of wind properties throughout the disk. Bai (2013) derived the following fitting formulae for the radial dependence of wind mass loss rate (more accurately, rate of gas mass loss per unit area of the disk) and the wind-driven surface stress (i.e., the ɀϕ component of Maxwell’s stress tensor) for a protoplanetary disk (his Eqs. (9) and (10)): $$\frac{{\dot{\Sigma }}_{w,𝑔}}{{\rho }_0{c}_s}\approx 3.08\times {10}^{-5}\frac{1}{N}{\left(\frac{r}{\text{au}}\right)}^{0.70}{\left(\frac{{\beta }_0}{{10}^{5}N}\right)}^{-0.46},$$(13) $$\frac{T_{ɀ\varphi }^{\text{Max}}}{{\rho }_0{c}_s^2}\approx 1.16\times {10}^{-4}\frac{1}{N}{\left(\frac{r}{\text{au}}\right)}^{0.46}{\left(\frac{{\beta }_0}{{10}^{5}N}\right)}^{-0.66}.$$(14)

Here, the quantities on the left hand side are normalized with respect to local midplane density (ρ₀) and sound speed, while r is the stellocentric distance and β₀ is the ratio of gas to magnetic pressure at the disk midplane. The normalization factor N = Σ_g/Σ_MMSN applies to a disk which diverges from minimum mass solar nebula (MMSN) density distribution. With these two equations, one can implement a phenomenological prescription of disk winds; the ${\dot{\Sigma }}_w$ constitutes a sink term in the gas continuity equation, while a sink term in the momentum equation can be obtained from $T_{ɀ\varphi }^{\mathrm{Max}}$ , accounting for the loss of mass and angular momentum through the winds, respectively.

Before Eqs. (13) and (14) could be used in a magnetohydrodynamic model, however, they need several modifications, and here we elaborate on these changes. Firstly, the radial term in these formulae is essentially an MMSN radius, in the sense that it corresponds to the stellar mass of one $M_{\odot }\left(r\equiv r_{M_{\odot }}\right)$. A shearing box is assigned a radial location in the disk based on its average angular velocity, Ω. Thus, to modify Eqs. (13) and (14) for an arbitrary stellar mass, the $\left(r_{{\text{M}}_{\odot }}/\text{au}\right)$ term needs to be modified so as to have the same angular velocity: $${\text{Ω}}^2=\frac{G{M}_{\star }}{r^3}=\frac{G{M}_{\odot }}{r_{{\text{M}}_{\odot }}^3},$$(15)

where M_★ is the mass of the central star. This implies that the radial dependence in Eqs. (13) and (14) should be replaced by $$r_{{\text{M}}_{\odot }}=r{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-1/3}$$(16)

in order to generalize these equations for a system with varying stellar mass.

Calculation of mass and angular momentum loss due to winds as prescribed by Eqs. (13) and (14) requires the disk plasma β₀, that is, the midplane ratio of gas to magnetic pressure. Since FEOSAD evolves the magnetic field, β₀ can be obtained throughout the disk (Sect. 2.1). We note that β₀ is defined differently in Bai (2013) (Heaviside-Lorentz units, ${\beta }_{\text{wind\hspace{0.17em}}}=2P_0/B_{ɀ}^2$ and in the hydrodynamic model of Vorobyov et al. (2020b) (Gaussian units, ${\beta }_{\text{MHD}}=8\pi P_0/B_{ɀ}^2$, where P₀ and B_ɀ are midplane pressure and vertical magnetic field, respectively. In the context of calculating the wind quantities, β₀ = β_wind is used. When the mass loss rate due to the winds over the entire disk is calculated by post-processing FEOSAD simulations (i.e., Eq. (13) applied to surface density of a gravitoviscous simulation), we found that this integrated mass loss rate was consistently about two orders of magnitude larger than the accretion rate onto the star. Such a high rate of mass loss is unrealistic and not sustainable, as the entire disk will be lost within 10⁴ yr.

There are two reasons for a substantially large wind mass loss rate while using Bai (2013)’s equations as they stand. Firstly, we found that the local β₀ in FEOSAD simulations is generally much lower than the values widely accepted in the literature. Measuring the magnetic field strength in protoplanetary disks remains challenging, although it is expected to be strongly subthermal, with a typical T Tauri disk exhibiting β₀ ≲ 10⁴ (Vlemmings et al. 2019; Lesur 2021b). Similarly large values of β₀ ~ 10²− 10⁸ are often considered in shearing box simulations (e.g., Bai 2013; Simon et al. 2018), which are in alignment with global disk simulations that include non-ideal MHD effects (Masson et al. 2016; Xu & Kunz 2021; Mauxion et al. 2024). However, in FEOSAD simulations, β₀ is consistently lower, indicating the presence of a strong magnetic field. This discrepancy occurs due to the ideal MHD limit for evolving the magnetic field, where the accreting gas tends to drag the field with it, increasing B_z and thus decreasing local values of β₀ . Secondly, the shearing boxsimulations may severely overestimate the mass loss rate, as the mass loss rate consistently decreases with an increase in the box height in the direction normal to the disk midplane. This term is highly uncertain and Bai (2013) proposed an additional correction proportional to H_g/r, in order to account for this dependence on the finite size of the shearing box. In our model, we attempt to remedy the large magnetic fields in the MHD simulations via a constant correction factor, $$C_{\beta }=1000,$$(17)

that is multiplied with β₀ while calculating the winds. This factor reflects the limitations of current simulations and may be dropped in the future when the non-ideal MHD equations are solved. Parameter C_β is constrained such that C_β ×β₀ ~ 104 over majority of the disk, which is close to the expected values of β₀ in protoplanetary disks. However, we note that there is relatively large amount of uncertainty in this estimate from both observational and theoretical side. In order to account for the finite size of the shearing box, we include a term proportional to the local disk aspect ratio in the equation for wind mass loss ${\dot{M}}_w$: $$C_{\text{Box}}=C_H\frac{H_{\text{g}}}{r},$$(18)

where C_H is a correction factor with the default value of 0.5. When both the C_β and C_Box terms are considered, the integrated disk wind mass loss rate is about the same order of magnitude as the accretion rate onto the star. The correction C_Box only applies to the mass loss rate and not the surface stress, as the latter is well constrained in the shearing box simulations. This parameter is constrained by the total mass loss rate from the winds such that the latter is of the order of the accretion rate onto the central protostar (Pascucci et al. 2023). Concerning the total mass loss from the system, we additionally assume formation of a high velocity jet close to the star in all models, including those with only grav- itoviscous evolution. The jet is formed concurrently with the star and it lies inside the central sink cell, and 10% of the accreting mass is assumed to be ejected via the jet (Hartigan et al. 1994; Bally 2016).

The phenomenon of episodic accretion plays an important role in the long-term evolution of a protoplanetary disk. In particular, FEOSAD simulations show luminous MRI outbursts, which are analogous to FUor type eruptions observed in young stellar objects (Hartmann & Kenyon 1996; Kadam et al. 2020; Vorobyov et al. 2020b). These eruptions are thought to be accompanied by enhanced ejections, as indicated by observations of variability in outflow signatures, manifesting as clumpy structures in jets and outflows (Reipurth 1989; Vorobyov et al. 2018b; Fischer et al. 2023). Approximately 10% of the accreting mass may be expected to be lost via molecular outflows that are associated with the low velocity, wide angle winds (Audard et al. 2014; Cruz-Sáenz de Miera et al. 2023). For our model, we interpret the enhanced outflows during outbursts as magnetic disk winds that are linked to the increased far-ultraviolet (FUV) luminosity of the central accreting star (Bethell & Bergin 2011). The FUV ionization near the disk surface is essential for making the gas in the disk well-coupled to the magnetic field (Gorti & Hollenbach 2009). The efficiency of wind transport increases with the increase in FUV penetration depth, Σ_FUV, which is proportional to incident radiation (Simon et al. 2013; Bai 2016).

In the shearing box simulations that result in the relations specified in Eqs. (13) and (14), the FUV penetration depth is considered to be a constant (Σ_FUV = 0.03 g cm⁻², in accordance with Perez-Becker & Chiang 2011) and effects of changing Σ_FUV are not parameterized. To the first approximation, when Σ_FUV is increased by an order of magnitude, the wind mass loss increases by about the same factor, while the wind stress also increases, but only by a factor of a few (Bai & Stone 2013). This guides the calibration of the exponents in the following correction terms that are introduced in order to account for varying FUV luminosityin the wind mass loss rate and surface stress equations, respectively: $${\text{C}}_{\text{FUV},{\dot{M}}_w}=\frac{1}{2}{\left(\frac{L_T}{L_{\odot }}\right)}^{0.8},$$(19) $${\text{C}}_{\text{FUV},T_{ɀ\varphi }}=\frac{1}{2}{\left(\frac{L_T}{L_{\odot }}\right)}^{0.5}.$$(20)

Here, L_T is total, that is, the accretion plus stellar luminosity of the system. The exponent for the ${\text{C}}_{\text{FUV},{\dot{M}}_w}$ term is chosen to be marginally lower than unity, since a larger value results in a significant decrease in the wind velocity during an outburst (see Appendix A). We note that this dependence incorporates several assumptions implicitly, for example, the FUV penetration depth is proportional to the central FUV luminosity, which in turn is proportional to the total luminosity of the accreting star, along with the exact exponents of the wind mass loss and stress. Additionally, the dusty winds can self-shield the FUV irradiation; the winds from the inner regions can absorb FUV photons and interfere with the wind launching from the remainder of the disk (Panoglou et al. 2012; Rab et al. 2022). With the exponents in the above equations being less than unity, selfshielding is inherently assumed in the model. We neglect any additional uncertainties, for example, effects of photochemistry in the upper layers on the ionization (Bergin et al. 2003; Woitke et al. 2009). Inclusion of these terms (Eqs. (19) and (20)) results in an increased wind mass loss rate during outbursts that is approximately an order of magnitude less than the accretion rate onto the star. A caveat of this approach is that the radiation takes a finite time to reach the disk, which should cause a delay in producing magnetic winds. However, this delay is a small fraction of the dynamical time at a given radius, and thus, it can be safely neglected.

As the winds originate centrifugally from a protoplanetary disk, it is essential to determine the precise extent of this disk. During the collapse of molecular cloud to form a star, a gas parcel of specific angular momentum j will be accelerated toward the central star in a parabolic orbit. This orbit intersects the plane of the disk at $$r_{\text{eq}}=\frac{j^2}{G{M}_{\text{enc}}},$$(21)

where M_enc is the mass enclosed within the orbit of the gas parcel (Dominik 2015). A gas parcel can be part of the centrifugal disk only if it has sufficient specific angular momentum to resist the inward fall. This criterion implies that material at a given radius r forms a centrifugal disk, if the ratio r/r_eq is less than unity. In a protoplanetary disk, however, the gas experiences additional pressure support, and hence, it moves with a sub-Keplerian velocity. Thus, for calculating the outer boundary of the centrifugal disk, we use a less stringent criterion of r/r_eq < 1.2. The disk winds are active only within this region and are completely inactive outside. In addition, the disk winds are turned off below a gas surface density threshold of 0.1g cm⁻², in order to limit unintended consequences of the winds at large distances. We note that due to the diminishing gas surface density, the wind mass loss rate and the surface stress decrease with radius and become progressively insignificant. We obtain a single “centrifugal radius” (R_cf) by azimuthally averaging the extent of the centrifugal disk at a given time.

The action of the disk winds needs to be taken into account self-consistently in the hydrodynamic equations that are solved in FEOSAD. For this, the equations that need to be modified are – mass continuity equation for gas and small dust (Eqs. (1) and (4)), momentum equation (Eq. (2)) and energy transport equation (Eq. (3)). For solving the MHD equations, FEOSAD uses operator splitting technique, which divides problem into two partial substeps- the source and the transport step (e.g., Stone & Norman 1992). During the numerical calculations, the corrections due to wind mass loss and surface stress $\left({\dot{\Sigma }}_w\text{\hspace{0.17em}and\hspace{0.17em}}{T}_{ɀ\varphi }^{\text{Max}}\right)$ are computed in each cell at the beginning of the source step. Corrections to gas and small dust surface density as well as energy are done at the end of the source step, but before the transport step begins. A fraction, dΣ_w,, is removed from the gas surface density such that $$d{\Sigma }_{\text{w},\text{\hspace{0.17em}g}}={\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}g}}dt,$$(22)

where dt is the current time step. Even though the winds originate well above the midplane, aerodynamically coupled dust grains can get entrained in the outflows and leave the disk (Miyake et al. 2016; Booth & Clarke 2021). Hence, an amount proportional to the small dust-to-gas mass ratio is subtracted from the small dust surface density, which is $$d{\Sigma }_{\text{w},\text{\hspace{0.17em}d}}=\frac{{\Sigma }_{\text{d},\text{sm}}}{{\Sigma }_{\text{g}}}d{\Sigma }_{\text{w},\text{\hspace{0.17em}g}}.$$(23)

The surface stress $T_{ɀ\varphi }^{\text{Max}}$ has the units of momentum per unit area per unit time and it is directly subtracted from the ϕ component of the momentum equation. The internal energy equation is balanced in such a way that the local gas temperature (energy per unit mass) is unchanged when the effect of the wind is applied, which results in the correction term $$d{e}_{\text{w}}=\frac{d{\Sigma }_{\text{w},\text{\hspace{0.17em}g}}}{{\Sigma }_{\text{g}}}e.$$(24)

In order to prevent a sudden, large amount of mass loss during a time step, the mass lost due to the wind is capped at 5% of the gas mass contained within the cell. In such a case, the disk $T_{ɀ\varphi }^{\text{Max}}$ is capped such that the velocity of the escaping wind remains unchanged (see Appendix A for details).

To summarize, our model assumes that the magnetic winds are symmetrical across the disk midplane on both sides. We also assume that the wind properties can be expressed as power law functions of the local quantities and stellar parameters. The final modified fitting formulae describing wind mass loss rate and surface stress are $$\frac{{\dot{\Sigma }}_{w,𝑔}}{{\rho }_0{c}_s}=C_M\frac{C_H{H}_{\text{g}}/r}{N}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-0.23}{\left(\frac{r}{\text{au}}\right)}^{0.7}{\left(\frac{C_{\beta }{\beta }_0}{{10}^{5}N}\right)}^{-0.46}{\left(\frac{L_T}{L_{\odot }}\right)}^{0.8},$$(25) $$\frac{T_{ɀ\varphi }^{\text{Max}}}{{\rho }_0{c}_s^2}=C_{T_{z\varphi }}\frac{1}{N}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-0.15}{\left(\frac{r}{\text{au}}\right)}^{0.46}{\left(\frac{C_{\beta }{\beta }_0}{{10}^{5}N}\right)}^{-0.66}{\left(\frac{L_T}{L_{\odot }}\right)}^{0.5},$$(26)

where $C_{\dot{M}}=1.54\times {10}^{-5}$ and $C_{T_{ɀ\varphi }}=5.8\times {10}^{-5}$ are obtained by consolidating all other constants in the respective equations. With this framework of magnetic disk winds, all of the exponents in Eqs. (25) and (26) are informed by the shearing box simulations. The only “free" parameters are C_H and C_β, which arise from the discrepancy with respect to non-ideal MHD effects and uncertainties in the wind mass loss rate, respectively. We constrain these parameters with the help of observations in Sect. 3.3. We note that with the assumption of spatially uniform mass to flux ratio, $\lambda ={\Sigma }_{\text{g}}/B_{\text{ɀ}},{\dot{\Sigma }}_{w,𝑔}~{\Sigma }_{\text{g}}{c}_{\text{s}}\text{\hspace{0.17em}and\hspace{0.17em}}{T}_{ɀ\varphi }^{\text{Max}}~{\Sigma }_{\text{g}}^{5/4}{c}_{\text{s}}^{2/3}$, which explains how the disk wind properties vary with the local quantities.

2.3 Interface to PRODIMO

The disk structures computed by FEOSAD are passed to the radiation thermo-chemical disk model PRODIMO (Woitke et al. 2009, 2016, 2024) for post-processing. This way, we can predict the optical appearance of these disks at millimeter wavelengths in the continuum and in CO lines to simulate ALMA Band 6 and Band 7 observations, and read off the apparent sizes of the disks as the radii that encircle 90% of the respective continuum and line fluxes. We choose the inclination to be face-on to avoid the necessity to de-project while producing ray-traced continuum images and CO line maps.

The main challenge for this interface is that FEOSAD uses 2D polar coordinates (r, ϕ), considering vertically integrated quantities, whereas PRODIMO uses 2D cylinder coordinates (r, ɀ), assuming the disk to be axisymmetric. We have therefore developed the following scheme. First, the FEOSAD disk structure at a given time is azimuthally averaged. A PPD structure is not strictly axisymmetric; however, the azimuthal asymmetries increasingly become smoothed out as the disk evolves and the GI diminishes. Second, the following quantities are passed to PRODIMO: the gas column density Σ_g , the column density of small dust Σ_d,sm, the column density of grown dust Σ_d,gr, the viscosity parameter α_turb, the scale height H_g, and the maximum grain radius a_max. All these quantities are passed to PRODIMO as vectors on approximately 200 log-equidistant radial grid points used by FEOSAD. In addition, the stellar luminosity L_★, the effective temperature of the star T_eff and the stellar mass M_★ are passed as scalars. Third, PRODIMO sets up the 2D structure of the gas disk in (r, ɀ) on the provided radial grid points, while adding a few additional grid points by extrapolation near the inner rim so that the penetration of stellar light into the disk is traced correctly. We then use the radius-dependent scale height H_??(r) to set up the vertical gas density structure ρ(r, ɀ) at each radius using a Gaussian function.

The dust grains in FEOSAD have size distribution between a fixed minimum size a_min = 0.005 µm and a variable maximum size a_max , while assuming a continuous power-law size distribution of index –3.5. In each vertical column, PRODIMO sets up the unsettled grain size distribution function using 100 size bins, and then the grains are settled according to the prescription of Riols & Lesur (2018) in each bin with the settling parameter α_turb (Woitke et al. 2024). We use the DIANA standard dust opacities Woitke et al. (2016), assuming the grains to be composed of 60% silicate, 15% amorphous carbon, with 25% porosity and without any ice. It is possible that ice formation can locally increase the dust opacities by more than a factor of 100, especially at UV to mid-IR wavelengths (Arabhavi et al. 2022). However, the ice exists only deep in the disk (A_V ≳ 10), where it is protected from UV photons, that is, where the refractory dust is optically thick. Therefore, ice opacities do not have a strong influence on the resulting dust temperature structure and the spectral appearance of the disk concerning CO sizes and millimeter fluxes. Although some uncertainties remain when vertical mixing is included in the models (Woitke et al. 2022), in which case the grains can be mixed up faster than photodesorption can destroy the ice. Based on the settled dust structure, stellar and interstellar irradiation, and internal viscous heating, PRODIMO performs continuum radiative transfer calculations. This results in the disk internal dust temperature structure T_d(r, ɀ) and produces ray-traced continuum images and CO line maps using a polar grid in the image plane with 225 × 72 segments.

We use the large DIANA standard chemical network (Kamp et al. 2017) for the thermo-chemical modeling part, which has 235 chemical species, with reaction rates mostly taken from the UMIST 2012 database McElroy et al. (2013), with added H₂ formation on grains, a simple freeze-out and desorption ice chemistry, X-ray processes including doubly ionized species Aresu et al. (2011), excited molecular hydrogen, and polycyclic aromatic hydrocarbons in five different charging states; altogether there are 4832 reactions. The most relevant chemical processes for this paper are the CO photo-dissociation in the upper disk regions, and the CO freeze-out in the outer midplane regions. For simplicity we solve for the chemical concentrations of all ice and gas phase species in kinetic chemical equilibrium. Recent code improvements concerning escape probability theory and photo-rates are explained in Woitke et al. (2024).

The viscosity parameter α_turb(r) is also used for computing the viscous heating in PRODIMO. Here we use a 2D diffusion solver in the optically thick core of the disk as described in the appendix of Oberg et al. (2022) and use the following formula for the local dust heating (Frank et al. 2002), $${\text{Γ}}_{\text{col}}(r)=\frac{9}{8}{\alpha }_{\text{turb\hspace{0.17em}}}{H}_{\text{g}}^2{\Sigma }_{\text{g}}{\text{Ω}}^3,$$(27) $$\text{Γ}(r,ɀ)={\text{Γ}}_{\text{col}}(r)\frac{{\rho }_{\text{d}}(r,ɀ)}{{\displaystyle {\int }_0^{\infty }{\rho }_{\text{d}}}\left(r,{ɀ}^{\prime }\right)d{ɀ}^{\prime }},$$(28)

where Γ_col(r) [erg s⁻¹ cm⁻²] is the viscous heating rate of a column and Γ(r, ɀ) [erg s⁻¹ cm⁻³] is the local viscous heating rate per volume $\left({\displaystyle {\int }_0^{\infty }\Gamma }\left(r,{ɀ}^{\prime }\right)d{ɀ}^{\prime }={\Gamma }_{\text{col}}(r)\right)$. All PRODIMO models have been performed with git version 4fee3902 from 2023/06/19.

2.4 Initial conditions and model parameters

The MHD simulations in FEOSAD start with the gravitational collapse phase of a starless molecular cloud core. The initial mass distribution as well as angular velocity is consistent with axisymmetric core collapse, where the angular momentum remains constant and magnetic fields are expelled due to ambipolar diffusion. The collapse of the core continues such that a star is formed within the central sink of the computational domain. A surrounding protostellar disk simultaneously forms while the envelope is infalling and the system continues to evolve through the embedded phase. Detailed description of the initial structure and angular velocity profile can be found in Basu (1997) and Vorobyov et al. (2020b). Since the focus of this study is on the effects of the magnetic winds on the disk formation and its evolution, all simulations start with identical initial conditions, which are summarized in Table 1. The cloud core of 0.83 M_⊙ forms the initial gas mass reservoir and results in a young star of mass between 0.5 and 0.7 M_⊙, depending on the efficiency of the particular model. The ratio of kinetic to gravitational energy is set to E_rot/E_grav = 0.23%, consistent with the observations of prestellar cores (Caselli et al. 2002; Li et al. 2012). The disks formed with these conditions are gravitationally unstable, while being stable to fragmentation and clump formation (Vorobyov 2013). These particular initial conditions are chosen because formation of disks in similar conditions is extensively studied with FEOSAD simulations (e.g., Kadam et al. 2020; Vorobyov et al. 2020b). Additionally, we compare our results with the ALMA survey of low mass star forming regions, wherein the stellar mass of the sample centers at M0 type star of mass ≈0.5 M_⊙ (Ansdell et al. 2016; Manara et al. 2023).

The initial bulk dust-to-gas mass ratio, ξ_d2g, is set to the canonical interstellar medium (ISM) value of 1% and all the dust mass is in the fully coupled submicron-sized small particles. The initial cloud core gas temperature is 20 K, which was also the temperature of the background radiation (Kennicutt & Evans 2012). A uniform background magnetic field of 10⁻⁵ G normal to the plane of the disk is also assumed (Liu et al. 2021). The mass-to-flux ratio (λ_c = Σ_g/B_ɀ, in the units of the critical value 1/2πG^1/2) in dense star-forming cores is observed to be λ_c ≃ 2 (Crutcher 2012), with a cloud or core being supercritical for values of λ_c above unity (Mouschovias & Spitzer 1976; Nakano & Nakamura 1978). A commensurate value of λ_c = 2 was used for our earlier investigations; however, some flux should be lost during the collapse phase due to ambipolar diffusion (Shu et al. 1987; Dapp et al. 2012). Since FEOSAD does not include nonideal MHD effects, we can justify a somewhat larger value of λ_c = 5 and thus partially compensate for the model limitations. The resolution of all simulations is 184 × 128 in polar coordinates (r, ϕ), with a logarithmic spacing in the radial direction and uniform spacing in the azimuthal direction at each radius. The increased resolution in the r-direction is chosen to minimize the elongation of grid cells, preventing potential numerical artifacts. The convergence of FEOSAD code at different spatial resolutions has been previously confirmed (e.g., Kadam et al. 2021) and the current resolution sets dr = 2.6 × 10⁻² au at the inner computational boundary at 0.53 au. The location of this inner boundary or sink cell excludes the actual inner disk radius of a PPD occurring at ≈0.1 au, which is determined by the magnetospheric accretion or dust sublimation radius (Dullemond & Monnier 2010; Hartmann et al. 2016). However, the marginally larger boundary relaxes the Courant condition for computational feasibility of the simulations (Courant et al. 1928), while being sufficient for capturing important dynamical phenomena occurring in the inner disk, such as episodic accretion. A carefully constructed inflow-outflow boundary condition is implemented at the inner sink cell in order to prevent the artificial drop in the gas surface density (Vorobyov et al. 2018a; Kadam et al. 2022).

In this study, we present results from six simulations as listed in Table 2. The fiducial simulations, GV-2 and GV-4, evolve via the combined action of gravitational and viscous torques. The wind quantities in these models are calculated; however, no corrections are made to the disk evolution equations. The value of background αMRI in GV-2 is 0.01. As discussed in Sect. 1, protoplanetary disks may exhibit significantly lower value of turbulent α and this possibility is explored in model GV-4, with α_MRI = 10⁻⁴. Model WI-2 is identical to GV-2, wherein the disk winds are taken into account self-consistently. The gravitovis- cous forces are active in all disk models, including those where the magnetic winds are included. In the case of WI-4, winds are considered with a lower value of α_MRI , and thus, it is otherwise identical with GV-4. Model WI-3 considers winds with the intermediate value of α_MRI = 10⁻³. In the case of WI-3a, the wind parameters C_H and C_β are adjusted, in order to demonstrate how a better fit to the observations can be obtained. With this limited parameter space study some of the most important possibilities in the disk wind scenario are explored.

3 Results

3.1 Global picture of disk evolution

In this section we describe the process of disk formation and its evolution on a global scale. FEOSAD simulations with similar initial conditions, model parameters and disk physics have been studied previously in detail in our group. In particular, the fiducial simulation GV-2 is analogous to model 2 in Vorobyov et al. (2020b) and Kadam et al. (2022). The motivation behind reanalyzing the gravitoviscous evolution of GV-2 and GV-4 is to confirm that our current disk model behaves as expected with the inclusion of incremental improvements and establish it as a benchmark in order to study and contrast the effects of including the action of magnetic winds.

As soon as a simulation begins, the initial prestellar cloud core starts to collapse, since it is constructed to have a supercritical mass-to-flux ratio and is gravitationally unstable. The time is measured from the onset of collapse, that is, the start of the simulation. In the initial stages of the collapse, surface density increases obeying a nearly self-similar law (Basu & Mouschovias 1994; Basu 1997), while the rotation profile remains rigid within about 50 au. As the cloud continues to contract, the mass of the gas in the sink cell exceeds 0.02 M_⊙, leading to the assumption of second collapse due to dissociation of hydrogen molecule (Larson 1969). A protostar thus forms at the coordinate center at this time, with a point-like gravitational potential. The stellar luminosity now illuminates the disk and the gas temperature in the innermost regions exceeds 600 K due to both stellar irradiation and compressional heating. The infalling gas in the innermost regions soon starts to decelerate while approaching the star and a centrifugal disk larger than the sink cell forms at about 0.022 Myr. This initial evolution before the formation of the circumstellar disk is almost identical across all simulations.

Figure 1 depicts the gas surface density of the disk models as listed in Table 2, in the inner 800 × 800 au box, at different instances of time. We note that the intervals between these instances are not uniform, but increase as the disk evolves, so as to highlight the viscous spreading of the disk. The dashed line depicts extent of the centrifugal disk at a given time, as specified in Sect. 2.2. Consider the first two fiducial gravitoviscous models, GV-2 and GV-4. The general evolution of the two models looks similar at this scale. The early stages are dominated by GI, as the disk mass is a large fraction of that of the central protostar. This phase is characterized by formation of large-scale spirals, which are seen as azimuthal asymmetries in the gas surface density. With the current choice of initial conditions, formation of self-gravitating clumps is not expected. As the disk continues to evolve, its viscous spread is clearly seen in the increasing centrifugal radius, which marks the theoretical outer extent of the disk. Similar to Kadam et al. (2020), we define the Class 0/I boundary as the time when the mass of the star-disk system exceeds half the mass of the initial cloud core. The Class I/II boundary is considered to be the time when the envelope accretion rate drops below 10⁻⁸ M_⊙yr⁻¹ . In Fig. 1, the first two frames capture the disks in the embedded Class 0 and Class I stages, respectively. In the last three frames, a disk can be considered as a T Tauri object of Class II. These class boundaries occur at roughly the same time in all simulations.

The major difference between models GV-2 and GV-4 at this global scale is the extent of the disk, especially during the Class II phase. GV-2 uses the value of α_MRI as 0.01, while that for GV-4 is 10⁻⁴. We note that α_turb can achieve a much lower value in the dead zone (see Eq. (10)). However, outside of the dead zone (≳10au), the disk maintains the maximum possible value of α_turb = α_MRI. The higher viscosity in the bulk of the disk causes GV-2 to evolve on a faster viscous timescale and show a larger viscous spread. The azimuthal asymmetries and disk substructure in the outer regions also get smoothed out with higher viscosity, producing a nearly perfect axisymmetric disk. The centrifugal radius of the disk is several hundred au, which is comparable with the gas radii of some of the largest T Tauri disks such as DM Tau or IM Lup (Dartois et al. 2003; Cleeves et al. 2016). On the other hand, the disk in GV-4 shows a lesser viscous spread. Lower viscosity in the disk is unable to efficiently get rid of angular momentum and transport the gas inward. Thus, the gas remains in the disk for a longer time and forms a smaller (≈250 au) but relatively massive disk, which shows spiral substructure even at 0.8 Myr. The disk radii and masses are discussed quantitatively in Sect. 3.3, with respect to synthetic observations.

The next four rows show models, WI-2, WI-4 and WI-3, and WI-3a, which include the disk winds self-consistently in the evolution equations. As specified in Table 2, the first three models have an identical prescription of magnetic disk winds, while WI-3a considers adjusted wind parameters. Consider the evolution of model WI-2; the disk shows asymmetry and spirals during the initial stage of evolution. During the Class II stage, the disk is much more compact as compared to analogous gravitoviscous model GV-2 and also smaller than the low viscosity model GV-4. In the last frame, the disk radius shrinks marginally, instead of expanding, resulting in a disk radius of ≈150 au. A purely viscous disk transports the angular momentum radially outward, redistributing it within the disk, which leads to a gradual expansion. Similarly, the gravitational torques work within the disk, but unlike viscosity, gravitational forces are inherently non-local and work over larger distances via non-axisymmetric perturbations and spiral density waves (Kratter & Lodato 2016). With winds, the disk evolution is driven by removal of angular momentum instead of its transport leading to so called “advective disk”, where the disk size tends to contract with time (Alessi & Pudritz 2022). In our models, the exact evolution of the disk size depends on the balance between gravitoviscous forces, which work toward expanding the disk and magnetic winds, which oppose this spread.

The evolution of the low and intermediate viscosity models, WI-4 and WI-3, is peculiar; we obtain centrifugal disks that are very compact. In the case of WI-4, the disk extends to a maximum size of ≈25 au during the Class I phase. However, it soon shrinks to a much smaller size of ≈6 au, which is unusually small for a Class II disk. The disk of model WI-3 achieves a maximum size of 45 au and then gradually shrinks to about 25 au. Although disk sizes of tens of au are well within the observational expectations for a gas disk, we later show in Sect. 3.3 that both models WI-4 and WI-3 produce dust disks that are unreasonably small in comparison with observations. The trend with respect to disk size across different α_MRI values confirms our hypothesis of the opposing action of disk winds and gravito- viscous expansion. With the current prescription of winds, their action is significantly stronger as compared to the viscous spread produced at low α_MRI values of 10⁻⁴−10⁻³.

With WI-3a, we explore the possibility of changing the wind strength by adjusting the wind-related parameters C_H and C_β . For this model, a larger value of C_H results in a more efficient mass loss. The larger value of C_β used implies weaker winds, as it results in a decrease in both the wind mass loss and the surface stress (see Table 2 and Eqs. (25), (26)). As seen in the last row of Fig. 1, the disk in WI-3a achieves a larger size as compared to the other two wind models at low viscosity, while also reaching a maximum size that is smaller than the GV models. This again indicates that the effect of reducing the magnetic wind efficiency via the C_β parameter is as expected; the balance shifts in favor of viscous expansion and a larger disk size is obtained. We note that the size estimates discussed here come from the theoretical centrifugal disk and thus are not directly observable. We quantitatively discuss the disk sizes and compare synthetic observations with demographics from an ALMA survey in Sect. 3.3. Additionally, we note that the gas surface density depicted in Fig. 1 gives us only a limited understanding of the total disk masses. For example, the disk sizes of models GV-4, WI-2 and WI-3a are comparable at any given time. However, their masses diverge significantly after the Class 0 stage. At one Myr, GV-4 produces a massive disk of 0.2M_⊙, which is nearly twice as massive as the disk in WI-2, and over 7 times more massive than that in WI-3a. The evolution of disk masses is discussed in detail in Sect. 3.3 (see Fig. 5 and Table 4).

Fig. 1 Evolution of the gas surface density for the protoplanetary disk models showing the global picture. The white contours show the extent of the centrifugal disk. We note that the time intervals are not uniform and are chosen to depict viscous spread of the disks. (See Table 4 for a summary of the properties of the resulting systems.)

3.2 Vertical and thermo-chemical structure of the disk

In this section, we discuss the results for the physical and thermochemical structure of a representative disk model obtained with the FEOSAD-PRODIMO interface. Figures 2 and 3 show the model WI-3a at 0.2 and one Myr, probing the disk as it progresses through the Class II stage. The first row of Fig. 2 shows the evolving surface density profile Σ_g(r). The early phases are featured by a relatively flat surface density between about one au and ten au, with a number of transient time-dependent pressure bumps, which are much more pronounced in dust than in gas. At later stages, only one pressure bump remains, whereas the disk outside of about three au follows a structure as expected from viscous spreading (see, e.g., McCaughrean & O’Dell 1996; Hartmann et al. 1998), with an approximate radial power-law surface density profile and an exponential decay beyond a critical radius of about 50 au. Inside of about one au, the surface density profile increases with the radius in both gas and dust at all evolutionary phases. The overall radial structure is similar to previous investigations with FEOSAD, where concentric gas and dust rings form inside of the innermost five au region of the dead zone due to viscous torques (Kadam et al. 2019, 2022). These regions are expected to host streaming instabilities and produce planetesimals (Yang et al. 2017).

A salient difference between an early (0.2 Myr) and an evolved Class II disk (one Myr), according to this model, is the distribution of the dust. At 0.2 Myr, the global dust-to-gas mass ratio only marginally deviates from the initial value of 1%, with the total dust and gas masses of 9 × 10⁻⁴ M_⊙ and 0.16 M_⊙, respectively. In the outer disk regions, the local dust-to-gas ratio is also maintained and is only marginally less than the initial value. However, at one Myr, most of the dust is found in the vicinity of the pressure maximum near 1 au. Both, the total dust and gas masses, 8.6 × 10⁻⁵ M_⊙ and 0.029 M_⊙, respectively, have diminished during disk evolution approximately by a factor of 10. The outer disk regions are severely depleted of dust; the local gas-to-dust ratio beyond three au is ≲ 1:10 000. This is a consequence of the dust dynamics via its growth in the low viscosity environment and subsequent inward radial drift.

In the second row of Fig. 2, we compare the scale heights calculated from FEOSAD’s midplane temperature and taking disk self-gravity into account (blue), to that calculated from PRODIMO’s midplane temperature (red dashed), and a more sophisticated variant thereof (black dotted). In case of the latter, the equation of hydrostatic equilibrium is integrated upward into the optically thin region using PRODIMO’s T_gas(r, ɀ) and then fitted by a Gaussian of single temperature. In the early phase, the disk is taller in the vicinity of inner rings, because of the combined effect of viscous heating and large vertical dust optical depths, which hinder the viscous heat to escape and create a temperature inversion (Fig. 3). As the disk evolves, the viscous heating diminishes and the scale heights in the inner regions decrease. Hence, the disk becomes increasingly more extended and flared due to the impact of direct irradiation from the star. Since FEOSAD and PRODIMO use different opacities (see Sect. 2.1 and 2.3), it is no surprise that the independent computation of the scale heights by PRODIMO is not identical with the values passed from FEOSAD. However, the overall magnitude and shape of the calculated H_g(r) matches reasonably well, especially in the outer disk beyond 10 au. Since the inner rim near the sink cell is directly illuminated by the star in PRODIMO, it is much warmer and the scale height is much larger as compared to that in FEOSAD, where this radiative effect is not included.

The gas density structure shown in the third row of Fig. 2 directly follows from Σ_g(r) and H_g(r). We obtain maximum gas particle densities of about 10¹⁵ cm⁻³ in the pressure bumps around one au. The last row of Fig. 2 shows the gas-to-dust mass ratio after settling. We obtain local dust-to-gas ratios ≳1 after dust settling in the pressure bumps, which is consistent with the midplane values in FEOSAD. However, the dust-to-gas ratio is <10⁻⁴ at the disk surface where the optical radial extinction A_V,rad = 1, which is an important result concerning the formation of mid-IR molecular emission lines observable with Spitzer and the James Webb Space Telescope (JWST) (Woitke et al. 2024).

In Fig. 3, the first row shows the dust temperature structure as calculated by PRODIMO. Due to viscous heating and close- to-diffusive radiative transfer, the midplane temperature features several smaller variations between one and five au at an age of 0.2 Myr. Consequently, as seen in the isothermal contours, there are multiple icelines (or snow-surfaces in 2D) around 150 K behind the pressure bumps, beyond five au. At the later stage, the icelines sit right behind the first and only pressure bump at about 1.5 au. Considering the second row of Fig. 3, the CO molecules populate a sandwich layer between disk midplane and its surface in the outer disk beyond the pressure bumps. This narrow layer is (i) warm enough to prevent CO freeze-out (T_d ≳ 25 K), and (ii) protected enough from UV radiation to prevent CO photo-dissociation. The latter indicated by χ/n ≲ 10⁻⁴ line, where χ is the UV field strength normalize to standard solar neighborhood, and n is the total hydrogen nuclei density.

The simulated continuum and line images are shown in the lower half of Fig. 3. For deriving synthetic observations of the disk sizes in dust (R_dust,Band6) and gas (R_{gas,(co2–1)}), we consider the radius that contains 90% of the spectral flux, which is applied to the simulated continuum images and the CO line maps. The CO gas radii found this way are always significantly larger than the continuum radii (see Table 3), because of two effects. Firstly, the CO lines are an optically thick tracer of the gas, whereas, excluding the innermost regions with dust rings, the 1.3 mm continuum (ALMA Band 6) is an optically thin tracer for the dust. The dust signal vanishes proportional to Σ_d. However, the CO signal features a slower decline with radius, reflecting the gas temperature radial profile, which is shallower than that of Σ_g, until it vanishes quickly as the line becomes optically thin. This first effect is always at work and it produces larger gas radii as compared to dust radii. Secondly, at the later stages of the disk evolution, the continuum disk radii are also affected by the dust growth followed by inward radial drift in the outer disk regions. This effect is evident from the much lower dust-to-gas ratios ≲10⁻⁴ after 1 Myr in the outer disk (see the first row of Fig. 2). In addition, those few dust grains remaining in the outer disk are the small grains, which have only little millimeter opacity. According to the FEOSAD model, the grains at 50 au have a maximum size of 500 µm after 0.2 Myr, but only 16 µm after one Myr. At late evolutionary stages, this second effect dominates, leading to an apparently shrinking dust radius while the gas radius is still expanding.

Table 3 shows how the derived continuum and CO line radii depend on the applied cut-off flux level, and on the choice of the line and ALMA band. The results are found to be rather independent of the choice of the ALMA band. However, the choice of flux threshold level is important, in particular for the continuum, which implies a considerable measurement uncertainty for the continuum radii. According to our model WI-3a, different threshold flux levels result in rather similar dust radii 50–70 au after 0.2 Myr, but in a large variety of values between 20 and 180 au after one Myr. With time, the dust grows and drifts, moving the carriers of the millimeter opacity inward. Despite the action of the winds, the outer disk still experiences gravitovis- cous spread, which contains trace abundances of larger grains. This explains the spread in radii with threshold flux at later times and implies that observations with high signal-to-noise ratio may result in much larger continuum radii. In summary, from our reference wind model WI-3a, we conclude that the disks tend to become fainter in the continuum with age, whereas they become brighter and larger in the CO lines due to gravitoviscous spreading.

Fig. 2 Disk structure of the WI-3a model after 0.2 Myr and one Myr. The upper plots show the gas and dust column densities, where the green dot marks the popular MMSN-value of 1700 g/cm2 at one au. The second row of plots shows the scale heights Hg(r) as passed from FEOSAD, and the values calculated from the gas temperatures calculated by PRODIMO, using either the midplane temperature only (red dashed) or the Tgas(ɀ) structure up to a height where the radial visual extinction is AV,rad = 0.001. The lower two rows show the hydrogen nuclei density n<H> and the gas-to-dust mass ratio. Additional contour lines show the radial and vertical optical extinctions, AV,rad and AV,ver, respectively.

Fig. 3 Continuation of Fig. 2 showing the dust temperature structure (top) and the CO-concentration (second row) calculated by PRODIMO based on the WI-3a model. Additional contour times show the ionization parameter χ/n (UV field strength divided by total particle density) to indicate where about the CO molecules photo-dissociate, and a dust temperature of 25 K to indicate where CO freezes out. The lower diagrams show the simulated 1.3 continuum and CO 2–1 line maps, from which we derive a continuum radius of 89 au and 109 au at 1.3 mm after 0.2 Myr and one Myr, respectively, and a disk radius in the CO 2–1 line of 198 au and 375 au, respectively. Each radius encircles 90% of the respective continuum or line flux.

3.3 Magnetic disk winds and observational constraints

The disk size and the mass are some of the most fundamental characteristics of a protoplanetary disk and some of the observable population-level differences between disks evolving with and without magnetic winds manifest in these global properties. In this section we compare our magnetic disk wind models against those evolving only with gravitoviscous forces, in the context of observational constraints. Here we briefly summarize the complications of observing protoplanetary disks (see Miotello et al. (2023) for a recent review). When measuring disk size, the dust radius is typically measured from surface brightness profile of dust continuum emission in millimeter wavelengths; the most favored for ALMA surveys is Band 6 at 1.3 mm. The radius is defined empirically as the radius which contains a certain fraction of the total disk emission (68%, 90%, or 95%). The gas emissions from a protoplanetary disk can be measured via molecular line traces such as CO isotopologs. For measuring the disk size in gas, the zero moment map of ¹²CO (J=2–1) is typically utilized, which is optically thick and easier to detect in small amounts at large radii. The gas radius is calculated similar to the methods as described above and it is observed to be larger than the dust radius, typically by a factor of a few (Ansdell et al. 2018). As discussed in Sect. 3.2, the larger gas disk radii can be attributed to optical depth effects as well as dust drift.

With the assumptions of optically thin emission along with a prescription of opacity and disk temperature, the continuum flux can be translated into the dust mass (Hildebrand 1983), although these assumptions do not always necessarily hold (Dunham et al. 2014; Woitke et al. 2019; Haworth 2021). Direct measurements of the gas mass, which is expected to contain the majority of the disk mass, are difficult due to the faintness and inaccessibility of H₂ line observations. However, the method of measuring gas masses using a tracer such as CO lines has severe limitations arising from the uncertainty in the CO/H ratio in disks (more specifically, gas-phase [CO]/[H2] ratio, see Fig. 3), which is inevitably connected to the disk thermochemistry (Ansdell et al. 2016; Bosman et al. 2021). Traditionally, the dust mass measured from continuum flux is translated into the total disk mass through a multiplicative factor of 100, as indicated by the interstellar gas-to-dust ratio. Predictions from hydrodynamic simulations of PPDs, however, point toward significant local variations from this value due to dust dynamics (Takeuchi & Lin 2002; Vorobyov & Elbakyan 2019) Thus, these different approaches of measuring the mass do not agree with each other for PPDs around low mass stars, resulting in significant uncertainties.

With these caveats in mind, we compare our results with the latest compilation of ALMA survey data of Class II proto- planetary disks. The star forming regions considered here are Ophiuchus, Taurus, Lupus and Chamaeleon, which have similar, relatively young (~0.8–3 Myr) ages and form a near complete sample (Ansdell et al. 2016, 2018; Manara et al. 2023; Miotello et al. 2023). Since not all data is available for each of the star forming regions, Table B.1 summarizes the data used for comparison with the simulated models. The theoretical masses and radii resulting from simulations are not directly observable, which motivates us to make synthetic observations with PRODIMO as elaborated in Sect. 3.2. Here we compare the synthetic observations in terms of radii R_dust,Band6 and $R_{{\text{gas,\hspace{0.17em}}}^{12}\text{CO}}$, both calculated at 90% flux threshold. The dust content of the disk M_dust,Band6 is also compared with the survey data, while for the gas component, we use the line fluxes of ¹²CO and C¹³O at a standard distance of 150 pc.

In Fig. 4, we show the temporal evolution of the centrifugal radius as well as synthetic gas and dust radii for the models. The vertical dotted lines in these plots indicate the average location of the Class 0/I and Class I/II boundaries, as defined in Sect. 3.1. The top panel of this figure shows the theoretical size evolution of the disks in terms of R_cf , which is smoothed with a box-car average with a box-width of 0.05 Myr, in order to get rid of stochastic variations caused by disk size pulsation inherit to young stellar objects with variable protostellar accretion (Vorobyov et al. 2020a). The significantly smaller sizes of low viscosity disks with winds in WI-3 and WI-4 are apparent in these plots. Even in the early Class 0 and I stages, these disks are much smaller than their gravitoviscous counterparts. As also seen in Fig. 1, the disks in models GV-4, WI-2 and WI-3a cluster around 150–250 au, while GV-2 produces a large viscously spreading disk. The exact values of disk radii as well as some of the global disk properties at one Myr are listed in Table 4 for comparison.

The second and third panels of Fig. 4 show synthetic measurements of gas and dust radii, respectively. As elaborated in Sect. 3.2, the gas radius is approximated from ¹²CO(J=2–1) line emissions, while the dust radius is estimated from simulated 1.3 mm continuum observations, each with a 90% flux threshold. The synthetic observations obtained via PRODIMO are performed after the embedded Class 0 and Class I stages are over, when the disk may be visible through the nascent cloud code. Starting at 0.2 Myr, a total of five PRODIMO calculations are performed per model, which are spaced uniformly in time so that the observational trends can be followed as a PPD evolves through the Class II stage. The shaded regions near the right axis of these panels depict half-violin plots of observed gas and dust radius data. The gas radii pertain to Lupus cluster only (Ansdell et al. 2018), while the dust radii are obtained from all four aforementioned star forming regions (Hendler et al. 2020; Manara et al. 2023). These violin plots show the probability distribution functions for the data, which are terminated at the extreme values denoted by dashes, with the central dash showing the median value.

In Fig. 4, we observe that for each model, R_cf lies between the synthetic gas and dust radii at all times, with the synthetic gas radius consistently larger than R_cf . As discussed before, the gravitoviscous disks spread rapidly to large sizes, both in dust and gas, with the disk spread increasing with the viscous α-parameter. The models with magnetic winds produce gas disks that are smaller than corresponding gravitoviscous models. At low values of α_MRI, the wind models WI-3 and WI-4 result in very small continuum sizes, which lie well outside of observational bounds. When the disk winds are active, these low alpha models produce significantly smaller initial Class 0 disks. The subsequent advective action of winds overpowers the gravitoviscous spread, which results in overall smaller gas disks.

Considering the evolution of dust radius in the third panel of Fig. 4, the GV disks show a monotonic increase. With the inclusion of winds, the dust disk radius decreases as a disk enters Class II stage and afterward, it may increase depending on the model. Recent survey of Orion molecular clouds indicates a similar decrease in the mean dust disk radius with evolution, which supports the scenario of disk winds (Tobin et al. 2020). The impact of dust evolution on the disk radius in the presence of winds can be summarized as follows. The dust size in the outer regions of a PPD is typically drift limited, while in the inner regions the dust may grow to reach the fragmentation barrier (Birnstiel et al. 2016). As a consequence of the dust partitioning, in the winds in the outer disk evacuate a relatively larger amount of (small) dust per unit mass. The remaining grown dust is subject to efficient radial drift toward increasing gas pressure in the low α environment, since a lower α increases dust fragmentation barrier, dust maximum size, and the Stokes number. PPDs subject to photoevaporative winds cause a reduction of dust reservoir at large radii by a similar mechanism (Sellek et al. 2020). The bulk of the disk is thus depleted of small dust and grown dust drifts to the inner disk regions where it gets trapped in the gaseous rings, with a large dust-to-gas ratio. As the result, the dust disk size in WI-3 and WI-4 shrinks to only a few au in Band 6 continuum flux. With the adjust wind model WI-3a, the effects of wind are attenuated by increasing C_ß parameter by a factor of 5. This modulates the inward transport of gas and consequently the dust drift so that the continuum disk in synthetic observations are congruent with observations even at low α_MRI = 10⁻³. The dust dynamics with disk winds will be studied in detail in a subsequent paper.

Figure 5 shows the evolution of disk mass for the models. The true disk gas mass in the model ( M_disk) is plotted in the first panel, which is calculated by summation over the extent of the centrifugally supported disk. All disks achieve the maximum disk mass at the beginning of the Class I stage and afterward in the Class II stage, the masses decrease monotonically at different rates. As a consequence, the disk mass is typically a substantial fraction of the stellar mass during the Class I stage and the disk is prone to gravitational instability (see second column of Fig. 1). In addition to being large in size, the GV models typically produce more massive disks, especially at low viscosity. With a low value of α_MRI = 10⁻⁴ in model GV-4, the transport of angular momentum due to viscous torques is diminished, which over time results in the accumulation of material and a more massive disk than GV-2. Both models WI-3 and WI-4 produce disks with masses in the vicinity of 10⁻⁴ M_⊙ or 0.1 M_J. Thus, these disks do not have sufficient mass to produce giant planets, which accrete gas over timescale of several million years in the core accretion scenario (e.g., Pollack et al. 1996). However, we note that the planet formation possibly begins before Class II stage and at that time, the mass available for this process is much more than typical assumptions of MMSN disk (Tychoniec et al. 2020; Schib et al. 2021). Models WI-2 and WI-3a produce disks that are of the order of 0.01 M_⊙ and are most consistent with our understanding of PPD masses in the T Tauri phase.

A key trend observed in the first panel of Fig. 5 is that the lower viscosity disks with gravitoviscous evolution are more massive. However, the trend is reversed with the inclusion of disk winds; that is, lower viscosity disks become progressively less massive over time. This can be explained with the strong dependence of disk wind mass loss and torque on the gas surface density. The low viscosity creates both mass and angular momentum transport bottleneck, and forces the disk to have a larger Σ_g, which in turn produces strong winds with more torque as well as mass loss. The net effect is that the winds remain strong, until the surface density is reduced, which ultimately results in a less massive disk. This should increase the mass of the central star in wind models relatively rapidly, which is seen in Table 4.

The second panel of Fig. 5 shows the theoretical dust mass of the disk, calculated as the sum of both small and grown dust components. This is done with an exclusion of the innermost five au of the disk, in order to avoid the region where the disk shows formation of dusty rings (Kadam et al. 2019, 2022). These rings form pressure maxima and tend to foster dust growth up to meter-sized boulders, which are not observable and contain significant amount of dust mass (Vorobyov et al. 2024). Additionally, the innermost region is typically optically thick, and does not reveal its true mass in observed disks. Thus, removing the innermost region from our theoretical dust mass estimates seems to be the best compromise, for comparison with the synthetic observations. For models WI-3 and WI-4, the dust radius collapses after about 0.5 Myr (last panel of Fig. 4), implying that the calculated dust masses outside of five au are correspondingly small and not representative of the dust content of the disk. However, for the rest of the models, the theoretical dust masses serve as a good proxy for the bulk of the observable dust mass in the disk.

The third panel of Fig. 5 shows the synthetically derived dust masses of the simulated disks in comparison with observations. The red violin plot in this panel shows the observed distribution of dust masses as calculated from ALMA Band 6 continuum data, where all four nearby star forming regions are considered (Miotello et al. 2023; Manara et al. 2023). The abrupt cutoff at the lower end of the observed kernel density distribution is caused by the detection limit of the surveys. Similar to the ALMA surveys, we assume that the dust emitting at submil-limeter wavelengths is optically thin and it emits isothermally at a characteristic temperature of T_d = 20 K (Hildebrand 1983; Andrews & Williams 2005; Manara et al. 2023). This allows us to estimate the disk mass from the Band 6 (1.3 mm) flux (F_ν) as $$M_{\text{dust,Band6\hspace{0.17em}}}=\frac{F_{\nu }{d}^2}{{\kappa }_{\nu }{B}_{\nu }\left(T_{\text{d}}\right)},$$(29)

where, B_ν(T_d) is the Planck function at a dust temperature of T_d = 20 K, and κ_ν is an average dust grain opacity, calculated as κ_ν = 3.5(870 µm/λ) cm²g⁻¹ (Beckwith et al. 1990; Miotello et al. 2023). Considering the masses at the Class 0/I boundary, the synthetic dust masses at this early times are consistently larger than those indicated by recent survey (≈10⁻⁵−10⁻⁴M_⊙ for Orion molecular clouds (Tobin et al. 2020)), except in the case of WI-3 and WI-4. When compared against the observed distribution of dust masses of Class II disks, GV-2 results in a massive disk at the upper end of the observed range. The significantly larger gravitoviscously spreading disk in GV-2 model contains the largest amount of dust mass. We find that the dust masses inferred from PRODIMO fluxes for all other disks are consistent with the observed range of the surveys, with low viscosity disks with winds producing a better fit. Comparing with the theoretical dust masses in the second panel, the synthetic disk dust masses show an excellent agreement for all models throughout the Class II evolution. As expected, WI-3 and WI-4 diverge after about 0.5 Myr, for aforementioned reasons.

As mentioned earlier, gas mass of a PPD may be inferred by multiplying its dust mass by a factor of 100. If we calculate the gas masses this way, the estimates will follow the exact same trends as the third panel of Fig. 5, shifted up by the same factor. As estimating disk gas mass comes with significant uncertainties, for the most transparent and direct comparison against the observations, the last two panels of Fig. 5 show the ¹³CO and C¹⁸O line fluxes at 150 au. The observed ranges are obtained from ALMA survey data for Lupus and Chamaeleon star forming regions (Ansdell et al. 2016; Long et al. 2017). In general, the line fluxes from the models are larger than the observed distributions and there may be several contributing factors. Firstly, we note that in PRODIMO, CO isotope-selective processes are not considered and the abundance of CO isotopologues is inferred with fixed ratios of [¹³C]/[¹²C] = 0.0140 and [¹⁸O]/[¹⁶O] = 0.0020. The less abundant CO isotopologues are thought to become optically thick much deeper in a PPD, thereby tracing the entire thickness of bulk of the disk. However, the increase in line flux with time for the higher mass models suggests that both of these traces are optically thick and therefore follow the increase in disk size with gravitoviscous spread. The simulated disks evolve only up to one Myr, while the considered star forming regions may be up three Myr old. As a PPD looses its mass with age, this would decrease the observed CO fluxes with respect to synthetic observations. The PRODIMO fluxes are calculated for a face-on disk with zero inclination, which may contribute to a higher CO flux as compared with the observed sample of disks. Finally, we note that CO observations imply a bulk gas-to-dust ratio between 1–10 for most disks, which is much lower than canonical MMSN or ISM value of 100 (Miotello et al. 2017). This discrepancy is attributed to carbon depletion in PPDs, with processes such as CO photodissociation in the upper layers and freeze-out near the midplane. If the CO emissions are optically thick and dust emission optically thin, this will also result in an underestimation of the measured bulk gas-to-dust ratio as compared to the actual value. The CO fluxes demonstrate that the adjusted disk wind model does not fit all observables simultaneously. Nonetheless, with the inclusion of disk winds, the simulated line fluxes decrease and move in the correct direction with respect to an agreement with the observed distributions.

We now test the disk evolution predicted by our models against observed spectral indices and the ratio of dust and gas disk sizes. Figure 6 compares the disk models with ALMA survey data of Class II PPDs in Lupus cluster (Ansdell et al. 2018). The top panel in Fig. 6 shows the total continuum flux versus far infrared spectral index between ALMA Band 6 (1.3 mm) and 7 (0.9 mm). The latter is the SED slope defined as $${\alpha }_{\text{mm}}=\frac{\log \left(F_{0.9\text{\hspace{0.17em}mm}}/F_{1.3\text{\hspace{0.17em}mm}}\right)}{\log \left({\nu }_{0.9\text{\hspace{0.17em}mm}}/{\nu }_{1.3\text{\hspace{0.17em}mm}}\right)},$$(30)

where F_ν is the continuum flux at the corresponding ALMA band. The α_mm = 2 line is characteristic of large particles represented by blackbody emission, since F_ν ∝ ν². The ISM value of 3.7 corresponding to small particles is also shown in the plot. Assuming optically thin emission, the flux in submillimeter continuum is a direct proxy of the total dust content in the disk. The spectral index probes the characteristic grain size and is hence often called the “growth index.” Although there is a higher degree of uncertainty, as additional factors such as optical depth and dust scattering complicate the picture (Draine 2006). PRODIMO calculates the fluxes at ALMA Bands at a distance of 150 pc and the spectral index is calculated from Eq. (30). We note that the observations have large error bars in spectral index, which are not shown in the figure. A disk starts at 0.2 Myr, when it is fully formed and typically evolves diagonally from lower right to upper left this phase space (diamond in the figure represents the end state at 1 Myr). Similar trends in the disk evolution have been previously shown in several studies (Birnstiel et al. 2010; Pinilla et al. 2012; Akimkin et al. 2020). The decrease in Band 6 flux indicates a decrease in dust mass or, in the case of optically thick emissions, an area effect due to inward dust drift. The evolution of spectral index is more complicated; the initial spectral index near 2 is a result of initial large optical depths and its value increases when the grain size approaches ≈λ/2π of the corresponding wavelength. The disks with large viscous a tend to show both larger flux and spectral index. The values of both spectral index and continuum flux obtained from the models are consistent with those in the ALMA surveys throughout the first Myr of evolution. The α_mm for WI-3 and WI-4 models first increases and then decreases with time. This decrease at later evolutionary stages is very likely a signature of dust clearing process, as all big grains drift inward to the pressure maxima (see Sect. 3.2). The α_mm is essentially disk-averaged spectral index and corresponds to the most luminous parts of the disk. As most of the dust mass in these low α models with winds eventually accumulates in the inner rings, the associated increase in the average optical depth may be responsible for the re-decrease in α_mm at later times.

The second panel of Fig. 6 compares the co-evolution of the gas and dust radii. These radii obtained from PRODIMO continuum emissions are same as those presented in Fig. 4. The observations suggest that a PPD has its gas radius larger than the dust radius by up to a factor of 4, as indicated by the diagonal dashed lines in the figure. Similarly, the disk radii inferred from simulations consistently show R_gas > R_dust. As elaborated in Sect. 3.2, a smaller dust radius can be attributed to dust growth and drift and it is also related to a drop in the dust opacity coefficient. The models with high α_MRI, GV-2 and WI-2 show evolution toward top right, as a consequence of gravitoviscous spreading. The low α models with winds, WI-3 and WI-4, deviate significantly from the observed trends. As discussed in the context of Fig. 4, the overpowering magnitude of the winds in this case results in small disks, and increased dust drift results in an increasingly compact disk in Band 6 continuum. We note that the only difference between models WI-2, WI-3 and WI-4 is viscous α. Thus, when the disk winds are involved, the exact magnitude of turbulence plays a significant role in overall disk evolution. The adjusted wind model WI-3a is much more consistent with the observational limits, as compared to WI-3. However, it still leaves some room for improvement, as the ratio of gas to dust radius is somewhat larger than R_g = 4R_d line. We note that the aim of model WI-3a is not to fit this wind model perfectly to the observations. Rather we adjust this model to show trends of varying the parameters associated with disk winds (C_β and C_H) and demonstrate that a disk with reasonable characteristics can be obtained at low alpha without fine tuning.

In Table 4, we list some of the global disk quantities that characterize the disks at 1 Myr. The numbers reinforce several conclusions made earlier and here we discuss some of the key points in terms of the mass budget of the PPDs. Model GV-4 produces a massive disk capable of forming envelopes of hundreds of gas giants, while WI-4 and WI-3 do not have enough mass to form a single gas giant. Most of the dust mass is in the grown dust, which is expected from dust dynamics occurring in a PPD. The low viscosity increases the fragmentation barrier for dust grains and the disk winds increase the inward dust drift, resulting in much more grown dust mass as compared to the small dust in these models. The global bulk dust to gas ratio (ζ_d2g = M_disk,dust/M_disk,gas) is approximately between 0.1 and 5% and in most cases, the dust is depleted below the canonical ISM value of 1%. The “efficiency” of star formation can be inferred from the stellar mass, as all models start from the same mass reservoir of 0.83 M_⊙. With gravitoviscous evolution, the lower α models produce a star with lower mass, as the diminished angular momentum transport maintains a relatively massive disk. With the inclusion of winds, however, we see that the efficiency of star formation is almost independent of disk α and varies with the strength of the winds. The total mass lost to the winds is inversely proportional to viscous α, as these disks produce stronger winds via its dependence on the gas surface density. However, a larger amount of net mass lost to the winds does not necessarily mean a less massive disk. For example, in the case of WI-3a, the wind stress is attenuated by increasing C_β, while the wind mass loss rate is increased by increasing C_H. Comparing against WI-3, the disk mass is larger, while the mass lost to the winds comes at the cost of the final mass of the star. Finally, the t_disk = M_gas/Ṁ_★, where Ṁ_★ is stellar accretion rate, is calculated at one Myr. This rough estimate of the disk lifetime implies that the disks with magnetic winds evolve faster with a noticeably shorter lifespan.

Fig. 4 Evolution of the disk size. Top panel: theoretical centrifugal disk radii. Bottom two panels: radii obtained from PRODIMO synthetic observations . The gas and dust radii for the models are calculated via synthetic 12CO and Band 6 continuum emissions with 90% flux threshold. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries. The shaded regions on the right depict half-violin plots of gas (red) and dust (red) radii of the Class II disks, as observed in a subset of ALMA surveys (see Table B.1).

Fig. 5 Evolution of the disk mass. Top panel: theoretical mass within the centrifugal region. Second panel: theoretical dust mass within the centrifugal disk, neglecting the innermost five au region. Third panel: Dust mass obtained from PRODIMO synthetic observations. Last two panels: synthetic CO line flux as a proxy for the disk gas content. The dust masses are calculated from continuum Band 6 fluxes, while the 13CO and C18O line fluxes are calculated at 150 pc. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries. The shaded regions on the right depict half-violin plots of dust masses (red) and line fluxes (blue) in Class II disks, as observed in a subset of ALMA surveys. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries.

Fig. 6 Comparison of the synthetic observations obtained from the disk models with the ALMA survey data with respect to the spectral indices (first panel) and dust-gas radii (second panel). The colored lines show evolutionary sequence of the models, starting at 0.2 Myr and ending at one Myr, with the latter point marked with a diamond. The blue scatter plots show observed quantities in Lupus cluster (Ansdell et al. 2018).

Fig. 7 Evolution of the disk models in the mass-radius plane showing how the wind models may diverge from the gravitoviscous models after the disk achieves its maximum mass. The solid lines depict theoretical disk size, Rcf, versus gas mass within the centrifugal disk. The dashed lines show $R_{\text{gas}{,}^{12}\text{CO}}$ obtained from synthetic emissions versus disk mass, the latter calculated as 100 times the dust mass obtained from the synthetic continuum emission. The diamonds mark the end points at one Myr and the blue scatter plot shows observed values from an ALMA survey of the Lupus cluster (Ansdell et al. 2016; Miotello et al. 2023).

3.4 Theoretical expectations from disk winds

In a recent review on protoplanetary disk demographics, Manara et al. (2023) illustrate some potential population level differences that can distinguish disks that evolve via predominantly viscous processes versus those driven primarily by disk winds. This difference arises because the action of magnetic winds tends to shrink the disk size as opposed to its viscous spreading, and it is best distinguished via the divergence of the evolutionary tracks in the mass-radius plane. Figure 7 depicts such a plot for the simulated models. Here the solid lines depict the theoretical centrifugal radius plotted against the disk mass, the latter corresponding to the gas mass enclosed within the centrifugal disk. R_cf is smoothed in time with a moving boxcar average with a box size of 0.05 Myr, hence these quantities correspond to the first rows of Figs. 4 and 5. The dashed lines show the disk sizes in ¹²CO, obtained from synthetic emissions and plotted against the synthetic disk mass at that time. The latter is calculated as 100 times the dust mass derived from the continuum emissions, since many studies still use this approximation. The diamonds mark endpoints of the simulations at one Myr. The scatter plot corresponds to the observed disks in Lupus cluster as seen in the ALMA survey data, where the disk masses are obtained from the Band 6 dust content (Ansdell et al. 2016; Miotello et al. 2023). The region on the extreme right represents the early and embedded GI-prone state when the disk is massive, while the trajectory will eventually move left toward lower mass and eventual dispersal, presumably via photoevaporation. The larger observed radius as compared to R_cf is seen in almost all models, which may be an artifact of a 90% flux threshold or additionally, of the criterion of surface density threshold in calculating R_cf. On the other hand, the observed disk masses are systematically smaller, which are a result of significant dust depletion as well as growth during the disk evolution.

The gravitoviscous models without disk winds spread monotonically, and end up with either a large disk radius (GV-2) or a large disk mass (GV-4). The low viscosity GV-4 model is able to hide most of its gas mass; however, such massive disks, which are a significant fraction of the stellar mass, typically exhibit GI spirals. Thus, GV models may be congruent with the observations of the largest or the most massive of the Class II disks. When both gravitoviscous and wind-driven accretion are present, the disk spreading is reduced and the radius evolution is essentially a result of the tug-of-war between these two processes. Hence, with the inclusion of disk winds, PPDs tend to evolve toward the diagonal in Fig. 7 so that they end up smaller and less massive. The disk wind models with low alpha, WI-3 and WI-4, produce compact disks that lie outside of observational expectations, although since $R_{\text{gas}{,}^{12}\text{CO}}>R_{\text{cf}}$, model WI-3 moves toward the observed values. Model WI-2 with a larger value of α as well as WI-3a are more or less congruent with the observed ranges. R_cf in WI-2 contracts near the end of its evolution, implying the size does not necessarily increase monotonically when disk winds are considered, although the synthetic disk size in $R_{\text{gas}{,}^{12}\text{CO}}$ increases with time. These results are qualitatively consistent with the expectations from previous models of disk winds, and imply that magnetic disk winds may be necessary to explain observed distributions (e.g., Bai et al. 2016; Tabone et al. 2022). Figure 7 demonstrates that in reality, both the measurements of disk sizes and masses contain several assumptions and significant uncertainties to draw decisive conclusions about the disk driving mechanism.

At this point, we discuss the dimensionless quantity of magnetic lever arm, λ, which is typically assumed as a parameter for wind-driven accretion models that are analogous to simple viscous α-disks (Kimmig et al. 2020; Alessi & Pudritz 2022; Tabone et al. 2022). This parameter characterizes the removal of angular momentum and it is essentially the ratio of specific angular momentum evacuated by the wind to that in the underlying Keplerian region that launches the wind. In the classical picture of magnetocentrifugal winds, the magnetic field lines rotate rigidly inside of the Alfvén surface so that the gas leaving as a wind increases its specific angular momentum (Blandford & Payne 1982). If a field line anchored at a radial distance r₀ crosses the Alfven surface at r_A, the magnetic lever arm can be written as $$\lambda ={\left(\frac{r_A}{r_0}\right)}^2=\frac{1}{2}\frac{{\dot{\Sigma }}_{\text{w},\text{tr}}}{{\dot{\Sigma }}_{\text{w},\text{\hspace{0.17em}g}}},$$(31)

where Σ_w,tr is the radial mass transport rate resulting from the magnetic winds. The lever arm decides how efficient the disk winds are in angular momentum removal, since with a larger λ, a small amount of mass lost carries a relatively large amount of angular momentum. In general, this parameter has to be greater than unity in order to extract angular momentum and brake the disk. Recent MHD simulations of weakly magnetized disks suggest λ to be between 1.1 and 5, with a typical value of 1.3 (Wang et al. 2019; Gressel et al. 2020; Lesur 2021a). Less efficient winds at λ ≈ 1 may primarily be driven by vertical gradients of magnetic pressure. Conical molecular outflows have been observed in a handful of sources, which give estimates for λ that are consistent with the simulations (de Valon et al. 2022; Louvet et al. 2018; Tabone et al. 2017).

In our disk wind model, both the mass loss and angular momentum loss rates depend on local disk quantities and stellar properties, while a lever arm parameter is not used explicitly. However, it is possible to infer the value of λ by calculating it a posteriori from Eq. (31). Assuming steady state accretion and vertically symmetrical wind geometry, $${\dot{\Sigma }}_{\text{w},\text{tr}}=\frac{8\pi }{\text{Ω}}r{T}_{ɀ\varphi }^{\text{Max}}.$$(32)

In the simulations, the local values of λ fluctuate severely. However, a global λ can be calculated, which is averaged over the entire centrifugal disk at a given time. Figure 8 shows the evolution of this inferred λ, which is smoothed with moving boxcar average with a box length of 0.05 Myr. For all models, λ decreases initially, especially during the Class I phase, while its value remains remarkably constant at later times. During Class 0 and I stages, the disk is massive and also relatively compact, which results in large average gas surface densities. This implies greater wind torques at early times and a correspondingly large values of λ. For GV-2 and GV-4, where disk winds are not included self-consistently, λ is relatively large at later times as well. This may be a result of a larger average Σ_g, as disk winds are not included self-consistently in the GV models. For the disk wind models, the lever arms converge to λ ≈ 1–2, which is in excellent agreement with both expectations from MHD simulations and observational constraints. The global values of λ are relatively steady during Class II stage in all cases. This suggests that simpler disk wind models that start with an MMSN disk and adopt a constant lever arm may also be suitable for long-term evolution.

Fig. 8 Evolution of the inferred value of the effective lever arm for the disk models. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries and blue, dashed line indicates lever arm length of unity. The shaded region on the right shows the expected range of values.

4 Discussions and conclusions

Magnetic disk winds play a crucial role in the evolution of pro-toplanetary disks, as suggested by both observational evidence and theoretical insights over the past two decades. In this study, we have introduced a simple semi-analytical formulation of magnetic disk winds in the thin-disk limit. The disk wind model is based on the insights gained from local shearing box simulations. The results of such local simulations can be stitched together to obtain analytical radial dependence of the wind properties throughout the disk (Bai 2013). These fitting formulae assume a power-law dependence of wind mass loss rate and surface stress as a function of local disk properties. We improved upon these results and took into account several corrections, for example, due to changing stellar mass and FUV luminosity. We incorporated this model of magnetic disk winds self-consistently in the FEOSAD code, which simulates the formation of protoplanetary disks and their evolution over the long term.

The wind mass loss rate and torque result in a sink term in mass continuity and momentum equations, respectively, while acting over the extent of the centrifugal disk. The corrections to other relevant equations, such as energy conservation and those for small dust, were also considered. The two free parameters in the disk wind model are C_β, which arises from the ideal MHD approximation used for the disk evolution, and C_H, which originates in the relatively unconstrained mass loss rate as predicted by the shearing box simulations. The simulation results depend sensitively on the factor C_β, which scales the disk plasma β₀ in our ideal MHD models. However, we note that this parameter is not arbitrary, as it is chosen so that the quantity β₀ × C_β over majority of the disk is of the order of 10⁴. This expected value is consistent with β₀ in protoplanetary disks from both observations and simulations, although these estimates vary by orders of magnitude. If our model is correct, only a narrow range of C_β may be suitable for disk evolution, implying that there could be an underlying mechanism – perhaps some self-regulation process during magnetic transport – that causes different protoplanetary disks to have similar magnetic β₀ values. Moreover, C_β may not be constant in space and time, and this uncertainty makes it challenging to compare models against observations. For the model to work without C_β, solving the non-ideal MHD equations is necessary. An alternative approach could be capping B_ɀ to a reasonable upper limit corresponding to observational constraints, for example, from paleomagnetism of meteorites in the early Solar System (Fu et al. 2014; Borlina et al. 2021). This approximation handles the two limiting cases of no and strong diffusivity correctly, as B_ɀ will be approximately constant in the limit of strong magnetic diffusion. Despite the aforementioned drawbacks, our model gives us a reasonable insight into the long-term evolution of wind-driven disks.

Previously, FEOSAD considered gravitoviscous evolution, wherein a PPD forms from the collapsing parent cloud and subsequently evolves via combined action of turbulent viscosity and gravitational torques. In this framework, we conducted numerical experiments that incorporated the effects of disk winds. In our limited parameter space study, three values of viscous α were considered (0.01, 10⁻³, and 10⁻⁴). Synthetic observations were obtained by post-processing the simulation results with the radiation thermo-chemical code PRODIMO. For comparison with observations, we aggregated available ALMA survey data for nearby star forming regions. We validated our magnetic disk wind model with respect to observational constraints, while the disk properties and characteristics of the magnetic winds will be investigated in detail in a subsequent paper. Without discussing the individual models in detail, our salient findings can be summarized as follows:

1. Protoplanetary disks evolving under gravitoviscous torques only (without disk winds) typically result in Class II disks that are too large in mass, and often in size, in comparison with observational expectations.

2. With the inclusion of magnetic disk winds, the resulting disks are smaller and less massive. This is expected, as advective evolution caused by vertical angular momentum transport from disk winds tends to suppress the outward gravitoviscous spread and promote inward mass transport.

3. Some of the models with winds produce Class II disks that are unreasonably small in terms of dust disk size as well as the true disk mass. Obtaining a reasonable fit to the observations is possible with a few adjustments in the two disk wind model parameters C_β and C_H.

4. The synthetic observations of global disk properties – gas radius in ¹²CO, dust radius in Band 6 continuum, dust content as calculated by Band 6 flux, and ¹³CO and C¹⁸O line fluxes – move in the direction of better agreement with the ALMA survey data when disk winds are included.

5. Similar to previous gravitoviscous simulations, models including winds also show the formation of a dead zone, pressure maxima, and dusty rings with an enhanced dust-to-gas ratio in the inner disk region. These can assist rapid planetesimal formation and point toward planet formation starting early in the Class I stage.

6. The divergence of evolutionary tracks with respect to disk winds may be seen on the disk mass versus radius plane, where careful choice of disk wind parameters leads to models that are more consistent with the observed ranges.

7. The ratio of synthetic gas to dust radii stays within the observed ranges of one to four for the gravitoviscous models. At a low viscous α, the models that include disk winds may produce anomalously small dust disks due to the loss of small dust and efficient migration of grown dust. By changing model parameters, the magnetic winds can be attenuated so that the resulting disk fits the observational constraints.

8. The synthetic spectral indices for all models are broadly consistent with the survey data throughout their evolution. As a disk evolves, an increasing α_mm indicates grain growth, and a decrease in Band 6 flux indicates a decrease in the dust content.

9. In agreement with simulations and observations, the inferred effective magnetic lever arm for an entire disk is between approximately one and three for the models that include the action of disk winds, and it stays remarkably steady during the Class II stage.

In conclusion, our magnetic disk wind model is able to reproduce the observational properties while simultaneously satisfying theoretical expectations over the long term (over one million years), wherein the disk forms and evolves through the Class 0, I, and II stages. Even now, most models of planet formation start with a fully formed MMSN protoplanetary disk and consider only its viscous evolution. Recent observations of PPDs put severe constraints on planet formation, as Class II disks do not contain enough dust mass to support the observed population of exoplanets (Greaves & Rice 2010; Tychoniec et al. 2020). This suggests that planetesimals and planetary embryos start forming earlier in the embedded stages, when the gravitational torques are significant, while disk winds control the fundamental disk properties, such as its size and the mass budget, affecting their subsequent evolution. Our model marks a step forward in the direction of representing a more complete picture of such disks, wherein the concurrent effects from viscous, gravitational, and magnetic disk wind processes guide the disk evolution and ultimately decide the fate of the planetary system that emerges.

Acknowledgements

We are thankful to the anonymous referee for constructive comments and suggestions. We thank Aleksandr Skliarevskii and Indrani Das for useful discussions. This work utilized computing resources provided by the Digital Research Alliance of Canada. E.I.V. acknowledges support by the FWF project I4311-N27.

References

Appendix A Escape velocity of disk winds

The escape velocity of the wind, v_w, constrains the relative strength of the wind mass loss rate and wind stress. It can be estimated to the first order by assuming that the work done by the wind torque equals the kinetic energy gained by the mass launched, resulting in $$v_w={\left(\frac{2{\tau }_{w}d\theta }{d{M}_w}\right)}^{0.5},$$(A.1)

where dM_w is the total mass lost to the wind (small dust and gas) during a time step. The angular distance traveled by the gas parcel during the time step can be calculated as $$d\theta =\arctan \left(\frac{v_{\varphi ,𝑔}dt}{R}\right),$$(A.2)

where v_ϕ,ɡ is the velocity of the gas parcel in the azimuthal direction. The wind torque is given by $${\tau }_w=\frac{1}{4\pi }{T}_{ɀ\varphi }^{\mathrm{Max}}RdA,$$(A.3)

where dA is the surface area of the grid cell. Thus, a moment arm length equal to the local radius is assumed. With the small angle approximation, Eq. (A.1) simplifies to $$v_w={\left(\frac{1}{2\pi }\frac{T_{ɀ\varphi }^{\text{Max}}}{{\dot{\Sigma }}_{w,𝑔}}{v}_{\varphi }\right)}^{0.5}$$(A.4)

Since the ratio $T_{z\varphi }^{Max}/{\dot{\Sigma }}_{w,ɡ}\approx c_s$, the wind velocity should also typically be of the order of the local sound speed. Wind velocity constraints the relative strengths of wind mass and angular momentum loss and since it compares favorably with observations (~ a few tens of kilometers per second), it provides an additional sanity check for the disk wind model.

Appendix B Summary of observational data

All Tables

All Figures

	Fig. 1 Evolution of the gas surface density for the protoplanetary disk models showing the global picture. The white contours show the extent of the centrifugal disk. We note that the time intervals are not uniform and are chosen to depict viscous spread of the disks. (See Table 4 for a summary of the properties of the resulting systems.)
In the text

Fig. 2 Disk structure of the WI-3a model after 0.2 Myr and one Myr. The upper plots show the gas and dust column densities, where the green dot marks the popular MMSN-value of 1700 g/cm2 at one au. The second row of plots shows the scale heights Hg(r) as passed from FEOSAD, and the values calculated from the gas temperatures calculated by PRODIMO, using either the midplane temperature only (red dashed) or the Tgas(ɀ) structure up to a height where the radial visual extinction is AV,rad = 0.001. The lower two rows show the hydrogen nuclei density n<H> and the gas-to-dust mass ratio. Additional contour lines show the radial and vertical optical extinctions, AV,rad and AV,ver, respectively.

In the text

Fig. 3 Continuation of Fig. 2 showing the dust temperature structure (top) and the CO-concentration (second row) calculated by PRODIMO based on the WI-3a model. Additional contour times show the ionization parameter χ/n (UV field strength divided by total particle density) to indicate where about the CO molecules photo-dissociate, and a dust temperature of 25 K to indicate where CO freezes out. The lower diagrams show the simulated 1.3 continuum and CO 2–1 line maps, from which we derive a continuum radius of 89 au and 109 au at 1.3 mm after 0.2 Myr and one Myr, respectively, and a disk radius in the CO 2–1 line of 198 au and 375 au, respectively. Each radius encircles 90% of the respective continuum or line flux.

In the text

Fig. 4 Evolution of the disk size. Top panel: theoretical centrifugal disk radii. Bottom two panels: radii obtained from PRODIMO synthetic observations . The gas and dust radii for the models are calculated via synthetic 12CO and Band 6 continuum emissions with 90% flux threshold. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries. The shaded regions on the right depict half-violin plots of gas (red) and dust (red) radii of the Class II disks, as observed in a subset of ALMA surveys (see Table B.1).

In the text

Fig. 5 Evolution of the disk mass. Top panel: theoretical mass within the centrifugal region. Second panel: theoretical dust mass within the centrifugal disk, neglecting the innermost five au region. Third panel: Dust mass obtained from PRODIMO synthetic observations. Last two panels: synthetic CO line flux as a proxy for the disk gas content. The dust masses are calculated from continuum Band 6 fluxes, while the 13CO and C18O line fluxes are calculated at 150 pc. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries. The shaded regions on the right depict half-violin plots of dust masses (red) and line fluxes (blue) in Class II disks, as observed in a subset of ALMA surveys. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries.

In the text

Fig. 6 Comparison of the synthetic observations obtained from the disk models with the ALMA survey data with respect to the spectral indices (first panel) and dust-gas radii (second panel). The colored lines show evolutionary sequence of the models, starting at 0.2 Myr and ending at one Myr, with the latter point marked with a diamond. The blue scatter plots show observed quantities in Lupus cluster (Ansdell et al. 2018).

In the text

Fig. 7 Evolution of the disk models in the mass-radius plane showing how the wind models may diverge from the gravitoviscous models after the disk achieves its maximum mass. The solid lines depict theoretical disk size, Rcf, versus gas mass within the centrifugal disk. The dashed lines show $R_{\text{gas}{,}^{12}\text{CO}}$ obtained from synthetic emissions versus disk mass, the latter calculated as 100 times the dust mass obtained from the synthetic continuum emission. The diamonds mark the end points at one Myr and the blue scatter plot shows observed values from an ALMA survey of the Lupus cluster (Ansdell et al. 2016; Miotello et al. 2023).

In the text

	Fig. 8 Evolution of the inferred value of the effective lever arm for the disk models. The vertical, dashed gray lines mark the average Class 0/I and Class I/II boundaries and blue, dashed line indicates lever arm length of unity. The shaded region on the right shows the expected range of values.
In the text