Protoplanetary disk formation from the collapse of a prestellar core
AAstronomy & Astrophysics manuscript no. disk_formation_2020Feb © ESO 2021February 17, 2021
Protoplanetary disk formation from collapse of prestellar core
Yueh-Ning Lee , , , , Sébastien Charnoz , and Patrick Hennebelle , , Department of Earth Sciences, National Taiwan Normal University, 11677 Taipei, Taiwane-mail: [email protected] Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, UMR 7154 CNRS, F-75005 Paris, France Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS, F-91191 Gif-sur-Yvette, France IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France LERMA (UMR CNRS 8112), Ecole Normale Supérieure, 75231 Paris Cedex, FranceAccepted 14th Feb. 2020
ABSTRACT
Context.
Between the two communities of star formation and protoplanetary disk evolution, there has only been a few e ff orts tryingto bridge the gap between a collapsing prestellar core and a developed disk. While it is generally accepted since about a decade thatthe magnetic field and its non-ideal e ff ects play important roles during the stellar formation, simple models of pure hydrodynamicsand angular momentum conservation are still widely employed in the studies of disk assemblage in the framework of the so-called“alpha-disk” model due to their simplicity. Aims.
The goal of the present work is an attempt to revisit the assemblage phase of the protoplanetary disk, with more up-to-dateknowledge of the prestellar core collapse.
Methods.
We perform 3D magnetohydrodynamic (MHD) simulations with the ambipolar di ff usion and full radiative transfer to followthe formation of the protoplanetary disk within a collapsing prestellar core. The global evolution of the disk and its internal propertiesare analyzed to understand how the infalling envelope regulates the buildup and evolution of the disk. We follow the global evolutionof the protoplanetary disk from the prestellar core collapse during 100 kyr with a reasonable resolution of one AU. Two snapshotsfrom this reference run are extracted and rerun with significantly increased resolution to resolve the disk interior. Results.
The disk formed under our simulation setup is more realistic and are in agreement with recent observations of disks aroundclass-0 young stellar objects. The source function of mass flux arriving onto the disk and the radial mass accretion rate within the diskare measured and compared to analytical self-similar models based on angular momentum conservation. The source function is verycentrally peaked compared to classical hydrodynamical models, implying that most of the mass falling onto the star does not transitthrough the mid-plane of the disk. We also found that the disk mid-plane is almost dead to turbulence, whereas upper layers and thedisk outer edge are very turbulent, and this is where the accretion happens. The snow-line, situating at around 5-10 AU during theinfall phase, is significantly further away from the center than in a passive disk, while this result could be of possible numerical origin.
Conclusions.
We have conducted self-consistent protoplanetary disk formation from prestellar core collapse, taking into accountnon-ideal MHD e ff ects. We developed a zoomed rerun technique to quickly obtain a reasonable disk that is highly stratified, weaklymagnetized inside, and strongly magnetized outside. During the class-0 phase of protoplanetary disk formation, the interaction be-tween the disk and the infalling envelope is important and ought not be neglected. We measured the complex flow pattern, andcompared this to the classical models of pure hydrodynamical infall. Accretion onto the star is found to mostly depend on dynamicsat large scales, i.e., the collapsing envelope, rather than the details of the disk structure.
1. Introduction
Whereas it is commonly accepted that the planets of our SolarSystem (S. S.) formed within 1 to 100 million years (Myr), sev-eral lines of evidence show that planetary formation processesmay have started on a much shorter timescale. The Calcium-Aluminum-rich inclusions (CAIs) that define the time origin ofthe S. S. most likely have formed in a hot environment close tothe Sun. To explain their presence in carbonaceous chondrites farfrom the Sun, it is proposed that they formed in the disk duringthe collapse stage of the prestellar core, where there was signif-icant infall of the di ff use medium, and were subsequently trans-ported outward during the disk viscous expansion (see e.g. Yang& Ciesla 2012; Desch et al. 2018) or transported in disk winds orstellar outflow (Shu et al. 2004; Yang & Ciesla 2012; Pignataleet al. 2018). Moreover, simulations of stellar cluster formationwithin a molecular cloud showed that the typical infall time ontoindividual star-disk systems is about a fraction of a Myr (Padoanet al. 2014). The low abundance of water in the inner Solar System (Mor-bidelli et al. 2016; Hyodo et al. 2019) and the identification ofseveral isotopic reservoirs (Kruijer & Kleine 2017; Van Kootenet al. 2016) in the S. S. argue for an early formation of a dynam-ical barrier, like Jupiter, at about one Myr, while the buildingblocks of this latter must form much earlier. Finally, a recentstudy showed that disks do not contain enough detectable dustto form the exoplanet populations (Manara et al. 2018), imply-ing that either (1) most planets form in a few 100 thousand years(kyr) at the class-0 / I phase, or (2) a continuous dust replenish-ment process feeds planets over the disk lifetime (Manara et al.2018). All these elements raise a fundamental question: did theplanetary accretion start during the assemblage of the protoplan-etary disks or after, as is often assumed? If accretion did start inthe disk during the prestellar core infall, this would be a majorchange in our understanding of planet formation, as most stud-ies still presume initial conditions from the
Minimum Mass So-lar Nebula model (MMSN, Hayashi et al. 1979; Weidenschilling1977), which is pertinent only for an isolated disk.
Article number, page 1 of 27 a r X i v : . [ a s t r o - ph . E P ] F e b & A proofs: manuscript no. disk_formation_2020Feb
Protoplanetary disks are classically described using the “alpha-disk model” (Shakura & Sunyaev 1973), which considers an iso-lated disk. This model, whereas very idealized, is a useful tool tostudy planet formation due to its flexibility and is still the basisof most dust growth studies (see e.g. Birnstiel et al. 2012; Testiet al. 2014, for a review) and planet population synthesis relevantfor exoplanets (Alibert et al. 2013). Accretion and transport pro-cesses across the disk are described with the e ff ective viscosity ν = α c / Ω , where α is the dimensionless local shear stress co-e ffi cient that will be defined further in the text, c s is the thermalsound speed, and Ω is the Keplerian frequency.The accretion rate onto the star in steady state, ˙ M = πν Σ ,is also controlled by α . The disk is heated by both the proto-stellar radiation and the viscous dissipation (Hueso & Guillot2005). The value of α is linked to the intensity and the struc-ture of turbulence in the disk. It was thought during the last 15years that the source of turbulence is the magneto-rotational-instability (MRI), giving α > − (see e.g. Balbus & Hawley1991; Fromang & Nelson 2006) in ideal magneto-hydrodynamic(MHD) cases. However, it was later realized that non-ideal MHDe ff ects (ambipolar di ff usion, ohmic di ff usion) may preclude diskmagnetization, such that α may be in fact (cid:28) − in the planet-forming region (see e.g. Turner et al. 2014, for a review).Whereas such a low turbulence state may favor dust coagulationand planetesimal formation (Bai & Stone 2010; Youdin 2011), alow value of α cannot account for observed accretion rates ontothe star, especially for young disks, with 10 − M (cid:12) / yr < ˙ M < − M (cid:12) / yr (Hartmann et al. 1998).Studies of disks with non-ideal MHD show that the disk maybe, in fact, vertically stratified, with a low- α mid-plane, high- α upper layers that drives accretion (Turner & Sano 2008; Turneret al. 2014), and possibly additional MHD structures above a fewpressure scale-heights ( H ) such, as disk winds (see e.g. Bai &Stone 2013; Riols et al. 2016). These new understandings leadto descriptions that are more complex than the α -prescriptionand might be time- and position-dependent (Suzuki et al. 2016).Nevertheless, all these studies still consider isolated disks. While many works have addressed the issue of disk formation,only a handful of works have attempted to study early planetaryformation processes (dust coagulation, transport, planetesimalformation) during the collapse of the prestellar core that feedsthe growing disk. The main di ffi culties are (1) to model the infallof the envelope and (2) to connect this latter properly to the cir-cumstellar disk. Classically, the self-similar singular isothermalsphere collapse solution has been used (SIS, Shu 1977; Ulrich1976). Angular momentum conservation inside the system is al-ways assumed (a critical approximation that neglects magneticforces). Several works showed (Nakamoto & Nakagawa 1994;Hueso & Guillot 2005) that, under the above hypothesis, disksmay form with characteristics similar to the MMSN (Nakamoto& Nakagawa 1994) or to some disks observed in young stellarclusters (Hueso & Guillot 2005). Using the SIS infall solution,it was shown that CAIs could form close to the star, from con-densing gas, and then be transported away due to viscous diskrelaxation (Yang & Ciesla 2012; Pignatale et al. 2018).Dr ˛a˙zkowska & Dullemond (2018) showed that planetesimalsmay form during the gas infall. In such models, the infall fromthe envelope is represented as a source term for the disk surfacedensity, and the disk transport processes (controlled by α ) are set independently. For example, Hueso & Guillot (2005) considereda fully active disk, whereas Pignatale et al. (2018) considered adisk with a dead-zone. However these two critical aspects are infact linked. For example, as the gas falls onto the disk, it changesthe ionization state of the disk, and thus α . Secondly, the turbu-lent motions inside the disk may be inherited from the turbulenceor the asymmetry in the infalling envelope. In addition, due toMHD e ff ects, some outflow may be present, and part of the gasmay be driven above the mid-plane. Finally, the accretion shockof the infalling envelope may heat up the disk upper layers, thuschanging the thermal structure of the disk. Therefore, it is indis-pensable to properly study the gas infall, that is coupled to thetransport properties of the disk, and the resulting disk structurein terms of the surface density, temperature profile, and α . Vorobyov & Basu (2010, 2015) studied 2D thin-disk modelswith radiative cooling, mechanical heating, as well as irradia-tion from the protostar and the background. The magnetic tensorwas treated as a dilution factor of gravity and magnetic pres-sure a multiple of the thermal pressure for assumed frozen-infields. They showed that the accretion processes within the diskare intrinsically variable, while external infall and radiation doa ff ect the fragmentation within the disk. Vorobyov et al. (2015)studied the influence of environment on an embedded disk. Thesize of the disk is highly subject to the properties, the angularmomentum in particular, of the infalling material. The disk frag-mentation is a ff ected in consequence. The global mass accretionhistory onto the star is consistent with that onto the disk, whilethe burst modes are largely regulated by fragmentation withinthe disk. It is to be noted that this 2D model does not describethe vertical stratification and thus does not have a source termfor the surface density. Whether mass transits through the entiredisk before reaching the star or fall through upper layers cannotbe distinguished in such setup.Ku ff meier et al. (2017) have developed zoom-in techniqueto follow the evolution of nine sink particles at 2 AU in a gi-ant molecular cloud simulation of 126 AU resolution with idealMHD prescription, and showed that the sink mass accretion ishighly sensitive to large-scale environment and it not as simpleas a SIS collapse. Ku ff meier et al. (2018) further increased theresolution of one of those zooms to 0.06 AU and demonstratedthat the numerical convergence is unlikely reached. The short-term variation (bursts) of the mass accretion is much attenuatedwith increased resolution.The above works have demonstrated the challenges in thestudy of protoplanetary disk evolution during the embeddedphase. Compromises, among which thin-disk approximation andideal MHD, were made to alleviate computational cost. The purpose of the present paper is to revisit the classical purehydrodynamical model of disk assemblage (Nakamoto & Naka-gawa 1994; Hueso & Guillot 2005) with three-dimensional non-ideal MHD simulations, during the first 100 kyr of the buildingphase of the disk. We wish to investigate the following questions: – What is the mass flux onto the disk as a function of positionand time? – What is the mass accretion rate from the core onto the disk,and that from the disk onto the star?
Article number, page 2 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core – What controls the accretion rate onto the star? Is it the turbu-lent stress inside the disk, or is it the envelope infall rate? – How is the angular momentum lost during the infall and theaccretion? – What is the disk size and surface density? – What is the temperature profile within the disk? Where is thesnow-line? – How does α behave inside the disk?To avoid confusion, we will refer to infall as mass fallingfrom the envelope to the disk, and accretion as mass falling fromthe disk onto the star throughout this manuscript. The di ffi cultyin this kind of study lies in the fact that long-term evolution anddetailed structures need to be followed simultaneously, which isa computational challenge. The present work is an attempt inthis direction where we present a detailed study of the disk andinflow structure of an assembling protoplanetary disk. The pa-per is structured as follows: the numerical simulations are pre-sented in Section 2. In Section 3, we briefly review the classicaldisk models applied in planetary science. Detailed analysis ofthe simulation results are presented in Section 4. The results andcomparison with the α -disk model are discussed in Section 5.Comparison to other numerical works is discussed in Section 6.We conclude the paper in Section 7.
2. The numerical simulations
To study self-consistently the formation of the protoplanetarydisk, we simulate the collapse from the scale of a prestellarcore. The code RAMSES (Teyssier 2002; Fromang et al. 2006),which treats the MHD with an adaptive mesh refinement (AMR)scheme, is used. Since the initial prestellar core is of sub-parsecsize and the disk forms as a high density region of a few tensof AU, this collapse problem implies a wide dynamical range ofscales, equivalent to four orders of magnitude. Such type of sim-ulations has been proven to be computationally di ffi cult due tothe required high resolution, which means both high computa-tional power and large memory. Even with an AMR scheme, thehighly structured mesh configuration makes the parallelizationof the computation ine ffi cient. For this reason, the simulations inthis study are typically performed on about a hundred cores atmost. The initial condition is a prestellar core with flat density profile,0.02 pc (3800 AU) radius, and total mass of 1 M (cid:12) . The coreinitially has a weak solid-body rotation around its center. Theradial profile of the rotation is not a major concerns, since theinitial amount of rotation does not have strong impacts on theoutcome of the disk (Hennebelle et al. 2020a, paper I hereafter).Uniform magnetic field threads the core and is inclined by θ mag = ◦ with respect to the rotation axis. Magnetic field strength isdescribed with the mass-to-flux ratio, µ , which is the normalizedvalue with respect to its critical value. The value µ = /µ = . α vir = . β rot = .
04 the ratio ofrotational to gravitational energy. A weak m = ff ects ofthese physical parameters are discussed in paper I.We solve the complete MHD equations considering the am-bipolar di ff usion (ion-neutral friction, Masson et al. 2012, 2016;Marchand et al. 2016). The temperature is calculated with full ra-diative transfer with a flux-limited di ff usion scheme (FLD, Com-merçon et al. 2011), and a temperature floor is set at 20 K. Thethermal irradiation from the protostar is also considered usingthe prescription of Kuiper & Yorke (2013). The accretion lumi-nosity is not included here because of the large uncertainties re-garding its actual value (see the discussion in Hennebelle et al.2020b).In this study, the Ohmic dissipation and the Hall e ff ect arenot taken into account when treating the non-ideal MHD. TheOhmic dissipation is important at higher densities and inducesthe formation of a small secondary disk around the protostar(Tomida et al. 2015). Including this e ff ect might not a ff ect ourgeneral conclusions since the central high density part is highlydemagnetized as will be shown in Sect. 4.4.2. On the other hand,the Hall e ff ect seems to indeed play a role during the formationof the disk (Tsukamoto et al. 2015; Marchand et al. 2018; Zhaoet al. 2020), and this remains to be considered in future studies.It is however hampered by large uncertainties due to the limitedknowledge on the grain properties. The AMR scheme requires that the local Jeans length is resolvedby 20 cells. The base grid has 2 =
64 cells in each dimension,and the canonical run has maximal refinement level (cid:96) max = dx = .
93 AU. With the canonicalsetup, the interior of the disk is not very well resolved, and thedynamics within the disk is not fully physical. Therefore, theanalyses of this run focus more on the infall of the prestellarcore envelope onto the disk and the formation of this latter. Suchcanonical choice of resolution allows to follow the disk evolu-tion during almost 100 kyr (150 kyr after the simulation starts).Increasing the resolution requires much longer computationaltime, and thus we only rerun the simulation from some snap-shots with increased resolution to examine the consequences andto resolve the structure of the disk.At the center of the collapse, where the high density is nolonger resolved by the fluid description, we insert a sink parti-cle to trap the collapsed mass (Bleuler & Teyssier 2014). Thisis necessary because otherwise the warm and dense gas accu-mulated at the center will expand in a way that is not physical.The sink particle forms when the gas molecular number densityexceeds n thres = × cm − . It accretes from its surroundingwithin 4 dx . A threshold accretion scheme is used such that everycell having density higher than n acc = n thres / c acc , is setto 0.1 in this study. In the canonical run, the entire disk is refined at the highest reso-lution (2 ). The vertical structure of the disk is however poorlyresolved. To study the detailed dynamics of the disk interior,two restarts are performed at roughly 40 (R_40ky_ (cid:96)
18) and 80
Article number, page 3 of 27 & A proofs: manuscript no. disk_formation_2020Feb
Table 1.
Simulation parameters.
Columns from left to right: label, initial time, final time, restart snapshot origin, maximum level of refinement,physical resolution, and sink accretion density threshold. The canonical simulation is evolved during roughly 100 kyr after the sink formation.This corresponds to 600 orbits at 1 AU distance and 10 orbits at 20 AU for a star of 0 . M (cid:12) . Two high-resolution restarts and one run withoutstellar irradiation are performed. With increased resolution, the time steps are significantly reduced and the time span of the restarts are very short. Label t int (kyr) t fin (kyr) Origin (cid:96) max Resolution (AU) n acc (cm − )R_ (cid:96)
14 0. 151. - 14 0.93 1 × R_40ky_ (cid:96)
18 102.532 103.144 R_ (cid:96)
14 18 0.06 1 × R_80ky_ (cid:96)
18 138.203 138.519 R_ (cid:96)
14 18 0.06 2 × R_ (cid:96) (cid:96)
14 14 0.93 1 × (R_80ky_ (cid:96)
18) kyr after the sink particle formation (at simulationtime 61 kyr) with increased resolution. The model parameters aresummarized in Table 1.Figure 1 shows four column density snapshots of the diskwith edge-on and face-on views. The first snapshot correspondsto the time right after the sink particle formation. Shortly beforethe sink particle forms, the dense core grows and flattens dueto the infall of the envelope. At the same time, the accreted gaswith angular momentum leads to the formation of a rotating pro-toplanetary disk around the central dense core. This core-diskstructure is initially thick, and reaches mass of almost 0 . M (cid:12) .Once the central star forms, the majority of mass of this densestructure is quickly accreted, and the disk mass drops and staysat around 2 −
3% of the stellar mass during the following evolu-tion.Despite the inclined initial magnetic field, the disk is almostaligned with the prescribed rotation. In the long-term evolution,the disk gradually grows in size from a few tens to almost a hun-dred AU at simulation time 120 kyr, while there are fluctuatingbehaviors due to the turbulent nature of the infalling envelope.The disk sometimes decreases in size due to sudden strong ac-cretions, that are often asymmetric, i.e., coming from one sideof the disk, and due to the interchange instability. We also seethe expulsion of magnetic field as an outflow in form of loops orbubbles (see e.g. second row of Fig. 1), which is highly dynamic.Figure 2 shows the restart at 40 kyr along with the canonicalrun. With increased resolution, the disk structure is better re-solved, while the global profile stays similar. The refinement isincreased by steps, passing progressively from level 15 to level18, in order to avoid numerical artifacts due to abruptly increasedresolution. The final resolution at level 18 equals to 0.06 AU.When increasing refinement levels, we had to adjust n acc corre-spondingly to avoid artifacts. As explained in paper I, the diskmass turns out to be very sensitive to n acc . Estimating with asimple analytical model, n acc should scale roughly as dx − .We caution here that the value of n acc likely sets an innerboundary condition for the disk and determines the density pro-file throughout the disk, and the temperature is a ff ected in turn(see discussions in Hennebelle et al. 2020a). Machida et al.(2014); Vorobyov et al. (2019) have also shown similar results.One should therefore be very careful when making physical in-terpretations of disk simulations. As we will show later, the e ff ectof the choice of n acc is seen in the disk global density, evolutionof the snow-line, and mass accretion rate onto the central star.
3. Existing models of disk formation
The collapse of the prestellar core is usually treated with ax-isymmetric and unmagnetized assumptions due to their simplic-ity and lack of observation evidences to support the contrary.Before analyzing our simulations, we briefly review what havebeen learned from the existing models.
One classical example of the collapse of a prestellar core is thesingular isothermal sphere (SIS, Shu 1977), which describes aninside-out collapse of mass shells that accrete onto a central starwith a time-independant mass accretion rate, ˙ M = . c / G ,where G the gravitational constant. While it should be kept inmind that this solution is just a special case among a wholefamily of self-similar solutions (e.g. Whitworth & Summers1985), it has been applied to many models for its simplicity.Furthermore, it is to be noted that, in reality, a prestellar corehas finite extent and the mass accretion rate onto the centralobject actually diminishes during the collapse and eventuallystops, giving an isolated disk. Indeed, observations have mea-sured mass accretion rate a few times larger than this canonicalvalue ( ∼ × − M (cid:12) / yr for isothermal gas at 10 K) for class-0 objects, that are at the beginning of their protostellar collapse(e.g., Ohashi et al. 1999; Hirano et al. 2002; Ward-Thompsonet al. 2007). Moreover, a declining mass accretion curve has beenobserved for young stellar objects between class-0 and class-Istages and simple accretion models have been proposed (e.g.,Henriksen et al. 1997; Whitworth & Ward-Thompson 2001).An easy, though not unique (see e.g. Verliat et al. 2020), wayto form a disk is to introduce a small amount of rotation into theprestellar core, under the assumption that the collapse solutionstays una ff ected at large scales. Due to conservation of angularmomentum during the collapse (with at least about a hundredtimes contrast in size), the rotation is significantly amplified andthis leads to the formation of a flattened structure in the center-most region around the star. While some numerical works exist(e.g. Li et al. 2014), very few studies have examined analyticallyhow a prestellar core collapses to form a star surrounded by aprotoplanetary disk because the breakup of the spherical sym-metry complicates the problem.This simplified picture consists of prescribing each collaps-ing mass shell with a certain angular velocity and deriving howthis mass is distributed radially when it arrives onto the equato-rial plane. The accretion onto the disk is expressed as the sourceterm S ( r , t ), that describes the mass flux as function of the radialposition and time. We discuss two widely accepted propositions.These models provide, in a deterministic manner, (1) the totalmass accretion rate onto the disk as a function of time, (2) theradius of the accreting zone of the disk as a function of time, and(3) the spatial distribution of mass flux onto the disk surface. As proposed by Ulrich (1976), in spherical coordinate (radius r ,latitude θ , and longitude φ ) all particles passing through radialposition r are assumed to have angular velocity ˙ φ and movein the plane instantaneously defined with the vectors r and φ . Article number, page 4 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core
100 50 0 50 100y (AU)10050050100 z ( A U ) Time = 61.096 kyr l o g ( k g m )
100 50 0 50 100x (AU)10050050100 y ( A U ) Time = 61.096 kyr l o g ( k g m )
150 100 50 0 50y (AU)050100150 z ( A U ) Time = 90.054 kyr l o g ( k g m )
500 450 400 350 300x (AU)15010050050 y ( A U ) Time = 90.054 kyr l o g ( k g m ) z ( A U ) Time = 120.034 kyr l o g ( k g m ) y ( A U ) Time = 120.034 kyr l o g ( k g m )
50 100 150 200 250y (AU)50100150200250 z ( A U ) Time = 142.264 kyr l o g ( k g m ) y ( A U ) Time = 142.264 kyr l o g ( k g m ) Fig. 1.
Column density snapshots of the canonical simulation R_ (cid:96)
Left: edge-on.
Right: face-on. The coordinates are in AU, and the originis at the box center. The images are centered on the sink particle andthe white dashed circle indicates the numerical radius of sink accretion(4 dx ). After the formation of the sink particle, the originally roundishfirst-Larson-core-disk structure is quickly accreted and a flat disk forms,with sightly smaller size. The protoplanetary disk does not stay at thecenter of the computation box and might experience a complex accre-tion history. The size of the disk grows slowly from its initial radiusaround 20 AU and presents fluctuations. At around 120 kyr, the diskexperiences a more rapid growth and reaches almost 100 AU in size.Although only a few snapshots are presented here, the actual picture ishighly variant in time. A parabolic trajectory, with the central star being the focus, isuniquely defined with these two criteria, while the non-zero ve-locity in the two other directions are not explicitly specified. Thefact that this spherical shell is not a rotating solid body did notcause any concern to the author, probably because the introducedamount of rotation is small. The particle is considered as ac-creted onto the disk when its trajectory intercepts the equatorialplane. The radial extent, or the centrifugal radius, of this region z ( A U ) Time = 103.161 kyr l o g ( k g m )
650 640 630 620 610x (AU)0102030 y ( A U ) Time = 103.161 kyr l o g ( k g m ) z ( A U ) Time = 103.144 kyr l o g ( k g m )
650 640 630 620 610x (AU)0102030 y ( A U ) Time = 103.144 kyr l o g ( k g m ) Fig. 2.
Zoomed views of the disk at age (103.1 kyr in simulationtime) from edge-on (left column) and face-on (right column).
Top row: R_ (cid:96) Bottom row:
R_40ky_ (cid:96)
18. The column density is shown as inFig. 1. The disk interior is clearly more resolved in R_40ky_ (cid:96)
18 with0.06 AU resolution, particularly in the vertical direction. of accretion is derived as r c = ˙ φ r GM , (1)which is the pericenter of the particle that arrives from the equa-torial plane, with M being the mass of the star-disk system,which is assumed to be dominated by the stellar mass.Assuming initial uniform distribution of mass on the shell,the source function can be derived as: S p ( r , t ) = ˙ M ( t )4 π r ( t ) (cid:32) rr c ( t ) (cid:33) − (cid:32) − rr c ( t ) (cid:33) − / , (2)where ˙ M ( t ) is the mass accretion rate of the mass that falls fromthe radius r , and is assumed to be a constant of time in most ofthe existing models as previously explained. On the other hand, Nakamoto & Nakagawa (1994) and Hueso& Guillot (2005) assumed that the particle, conserving its an-gular momentum, arrives at the radius of Keplerian rotation, re-gardlessly of the trajectory. With the exact setup of initial massand angular momentum distribution as the previous model, thesource function becomes: S k ( r , t ) = ˙ M ( t )8 π r ( t ) (cid:32) rr c ( t ) (cid:33) − / (cid:32) − (cid:114) rr c ( t ) (cid:33) − / . (3) In the model by Ulrich (1976), the particle has sub-Keplerianrotation when it reaches the disk and will fall inwards. The ex-change of angular momentum by friction of particles in the diskis not described, and a disk evolution model is further required tofollow the re-distribution of matter inside the disk. In the model
Article number, page 5 of 27 & A proofs: manuscript no. disk_formation_2020Feb r / r c S ( M t o t / ( r c )) Disk surface mass flux
Ulrich 1976Hueso & Guillot 2005 r / r c M ( < r , M t o t ) Cumulative mass accretion rate
Ulrich 1976Hueso & Guillot 2005
Fig. 3.
Upper panel:
Source function as proposed by Ulrich (1976) andHueso & Guillot (2005). The radius is normalized to the centrifugalradius r c , and the source term is normalized to the mean surface massflux ( ˙ M tot / ( π r )) such that its surface integral gives unity. Lower panel:
Cumulated mass accretion rate (normalized by the total mass accretionrate), ˙ M , integrated within radius r . The values of ˙ M tot and r c are to bespecified as functions of time according to the initial conditions of theprestellar core envelope. by Nakamoto & Nakagawa (1994) and Hueso & Guillot (2005),on the other hand, the trajectory is not described and the parti-cle reaches directly its centrifugal radius, giving a slightly morecentrally concentrated source function.The two models are shown in Fig. 3. For more intuitiveunderstanding, we also show the accumulated mass accretionwithin radius r :˙ M ( < r ) = (cid:90) r S ( r (cid:48) )2 π r (cid:48) dr (cid:48) , (4)which equals to zero at the origin and the total ˙ M at r c . The dif-ference between the two prescriptions do not seem to have strongimpact on the evolution of the disk (Hueso & Guillot 2005). Bothmodels rely on pure hydrodynamical arguments and the initialcloud is mapped to the disk deterministically, in a shell-by-shellmanner. The instantaneous source function depends completelyon the assumption of uniform mass distribution on the shell andconstant angular velocity all over the shell. While widely ac-cepted due to their convenience, these models are probably toosimplistic.The mass accretion rate ˙ M ( t ) and centrifugal radius r c ( t ) areboth determined by the initial conditions (density and angularmomentum profiles) of the core. The inside-out collapse solu-tion of Shu (1977) suggests a constant accretion rate ˙ M ( t ). Byassuming a rotation profile of the prestellar core (often solid-body which leads to r c ( t ) ∝ t ), the full accretion history of thedisk can be constructed. Simple variations of this class of mod- els can be prescribed with di ff erent initial density and rotationprofiles, leading to di ff erent histories of ˙ M ( t ) and r c ( t ).The source function is of particular interest for the studies ofearly S. S. evolution because it implies the temperature-pressurehistory that the matter inside the S. S. can experience. In a purelyhydrodynamic model, a non-negligible fraction of mass arrivesat the inner part of the disk, where the temperature is high dueto viscous heating and protostellar irradiation. The centrifugalradius grows rapidly in time such that later-on accretion arrivesat the cooler part of the disk. Such model has been used to ex-plain the composition diversity of chondritic meteorites (Pigna-tale et al. 2018). More and more observational evidences show that the star-forming gas is likely magnetized, and so is the disk since it isunlikely to completely get rid of the magnetic field at this scale.In the last decade, numerical simulations started to take into ac-count the non-ideal MHD e ff ects in the studies of early disk for-mation. Those authors found that the size of the disk is regulatedby the magnetic field and is small compared to what would havebeen if pure hydrodynamics is considered (see e.g. Dapp & Basu2010; Inutsuka 2012; Li et al. 2014; Tsukamoto et al. 2015; To-mida et al. 2015; Hennebelle et al. 2016; Machida et al. 2016;Zhao et al. 2018; Wurster & Li 2018; Lam et al. 2019). Thisis also in line with the observations suggesting that most class-0 disks are small ( <
50 AU, Andrews et al. 2018; Maury et al.2019).In this work, we aim to provide a revised model for the pro-toplanetary disk assemblage by taking into account the non-idealMHD e ff ects and the finite extent of the prestellar core. This canserve as more realistic framework for the future studies of theearly protoplanetary disk evolution.
4. Disk characteristics evaluated from simulations
Before performing the canonical run, we have simulated the col-lapse of a non-magnetized core. The outcome is very similar towhat is suggested by classical models: a disk forms due to an-gular momentum conservation and it size grows quickly in time.However, this disk is highly unstable and fragments very quickly(see e.g. Matsumoto & Hanawa 2003). Therefore, it is unlikelythat this scenario can lead to an extended and evolved disk. Ob-servations (e.g. ALMA Partnership et al. 2015) also suggest thatdisks are probably small at the early times. This is what moti-vated the inclusion of the magnetic field along with its non-ideale ff ects.In the presence of ambipolar di ff usion, the disk size is regu-lated by the competition between the magnetic braking that de-creases the angular momentum of the gas and the ambipolar dif-fusion that leaks the field lines outwards and reduces the braking(Hennebelle et al. 2016). Given the very di ff erent picture with re-spect to ideal MHD simulations, the non-ideal MHD e ff ects areindeed necessary for understanding the evolution of disk struc-ture. The magnetic flux is allowed to leak outwards, reducing thebraking, and a rotation-supported, demagnetized disk can formin consequence.The simulation starts with a core in free collapse, and thusit takes some time for the mass to accumulate before the pro-tostar and the disk form. When the central density increases,the gas heats up due to the increase of dust opacity, and thepressure locally supports against the self-gravitating collapse. A Article number, page 6 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core quasi-stationary first hydrostatic dense core (or first Larson core,Larson 1969; Masunaga et al. 1998; Vaytet & Haugbølle 2017),forms in consequence, and this core has a slightly flattened shapedue to the rotation. Further mass accumulation leads to the for-mation of a sink particle, that represents the protostar, at the cen-ter roughly 61 kyr after the beginning of the simulation.We evaluate several disk properties from the simulations. Allquantities are averaged in the azimuthal direction and betweenupper and lower halves of the disk.
The first step is to clearly identify the disk region and infer itsgeometrical properties. For this purpose, we analyze the distri-bution of mass as function of radius and altitude to derive thescale height, surface density, and radial extent of the disk. Theanalyses are done inside a selected cylindrical region with radius r cyl =
100 AU and half height z cyl =
100 AU, of which the cen-ter and the axis are first evaluated by selecting dense cells over acertain threshold as a first step of disk identification (see detailsin Appendix A). At later times when the disk has grown in size,the selected region is increased to r cyl = z cyl = H We calculate the azimuthally averaged vertical density profile forall radii. The binning size in the radial direction is taken to be dx ,the smallest cell size (or several times dx if the disk is large withrespect to the resolution). We evaluate the disk scale height withseveral approaches, which should all give identical results if thedisk is vertically isothermal and follows exact Gaussian densityprofile in pressure equilibrium. The thermal scale height fromvertical hydrostatic equilibrium (assuming a dominating stellarmass) is h p ( r ) ≈ (cid:115) c ( r ) r G [ M ∗ + M d ( r )] , (5)where c s ( r ) is the local sound speed, r the cylindrical radius, M ∗ the stellar (sink) mass, and M d ( r ) the disk mass inside r . Notethat the disk is not a point mass and the complete gravitationalfield should be calculated by integrating over the whole disk.We use this approximative formula for simplicity since the diskmass is only a few percent that of the star. Because the disk is notexactly isothermal, the thermal sound speed is calculated as c = (cid:80) PdV / (cid:80) dm , where P , dV , dm are the pressure, volume, massof each cell and the summation is performed vertically withincylindrical shells. Here again, we have verified that the verticalupper limit of the sum has no significant impact and we simplysum over the whole cylinder.In the case of an exact vertical Gaussian profile, ρ ( z ) = Σ tot √ π H exp (cid:34) − z H (cid:35) (6)the scale height can be recovered from the first or second mo-ment of mass: H = (cid:114) π Σ tot (cid:90) ∞−∞ ρ zdz = (cid:34) Σ tot (cid:90) ∞−∞ ρ z dz (cid:35) / . (7) r(AU)10 H ( A U ) h p r h h r(AU)10 H ( A U ) h p r h h r(AU)10 H ( A U ) h p r h h r(AU)10 H ( A U ) h p r h h R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 4.
Scale heights from vertical hydrostatic equilibrium, h p (blue),from the first ( h , orange) and second ( h , green) moment of mass, plot-ted against the radius. Top two panels: R_ (cid:96)
14 and R_40ky_ (cid:96)
18 at 103kyr. . Bottom two panels: R_ (cid:96)
14 and R_80ky_ (cid:96)
18 at 138 kyr. The hori-zontal and the left vertical gray lines mark the numerical resolution ( dx ,outside the figure extent thus not shown for the high resolution runs).The values are averaged within rings of thickness dx and plotted againstthe mean radius of the rings. The second vertical gray line marks theradius of the numerical sink accretion scheme (4 dx ). The two verticaldotted lines mark the characteristic radii of the disk ( R kep and R mag ), aswill be defined in Sect. 4.1.2. Within R kep the disk is close to verticalthermal equilibrium and the three definitions of scale height coincidevery well. Beyond this radius, magnetic support becomes important anddeviation occurs. The thermal scale height, h p ( r ), is fitted to a powerlaw,which is shown in the legend (in unit of AU). In practice, we calculate h ( r ) = (cid:114) π (cid:80) zdm (cid:80) dm and h ( r ) = (cid:34) (cid:80) z dm (cid:80) dm (cid:35) / , (8) Article number, page 7 of 27 & A proofs: manuscript no. disk_formation_2020Feb where the summation is performed vertically within cylindricalshells. These quantities are used as estimates of the thermal scaleheight.Figure 4 shows snapshots of the scale height calculated forthe two high resolution restarts compared to the canonical run.These results of h p , h , and h indeed reasonably coincide withone another inside the disk, and the significant incoherencemarks the radius at which the disk ends (left dotted vertical line).This also confirms that our choice of the cylindrical region sizeis su ffi cient to avoid errors, since the Gaussian density profiledrops very quickly at high altitudes.The two vertical dotted line denote the disk characteristicradii, that will be defined in Sect. 4.1.2 and further discussedin Sect. 4.1.3. The first radius, R kep , corresponds to the innerpart of the disk that is centrifugally supported in the radial direc-tion and thermally supported in the vertical direction (see Sect.4.4.1). The scale heights determined from three di ff erent defini-tions thus coincide very well with one another. The second ra-dius, R mag , corresponds to an outer region of the disk where themagnetic support is relatively important, and deviations start toappear. Beyond R mag , the moments of mass no longer give muchinformation as the infalling envelope does not have a flattenedgeometry. For a vertically uniform density profile (which is al-most the case in the outer part of the disk and the envelope), h = √ π/ z cyl ≈ . z cyl and h = √ / z cy ≈ . z cyl .They are almost indistinguishable either inside or outside thedisk, thus it is their discrepancy with h p that exhibits the clearchange to a non-thermally supported region.The overall behavior is similar between the high and low res-olution runs. However, the high resolution runs provide moreprecise descriptions of the disk height profile and allow to probethe high density interior of the disk. The horizontal gray line inFig. 4 marks the resolution dx . The disk in the canonical run con-tains only one cell in the vertical direction within 10 AU radius,while the high-resolution restarts resolve the disk down to radiusof 1 AU. Run R_40ky_ (cid:96)
18 shows a larger and better delineatedzone of transition between the two characteristic radii, while inthe canonical run this region is barely resolved by a few cells.The disk has nearly constant aspect ratio H / R ≈ . (cid:96)
18 resolvesbetter the disk center but does not give much extra informationat large radius since it is already well resolved in the canonicalrun. Σ The surface density of the disk is obtained as Σ cyl ( r ) = (cid:80) dm π rdr = Σ tot ( r ) erf (cid:32) z cyl √ (cid:33) ≈ Σ tot (r) (9)where Σ cyl is the column density within the cylindrical region,and Σ tot is the total surface density if the vertical density profilefollowed an exact Gaussian, with the mass contained within ± z cyl equaling to the numerically measured value. The error functionis a correction for the mass outside the cylinder vertical extentthat is not taken into account. However, this correction is almostnegligible for since H (cid:28) z cyl .Figure 5 shows snapshots of the Σ profile. At 138 kyr, thesurface density profile is close to a powerlaw with a sharp dropnear 13 AU, and becomes flat again outside 30 AU. The rapidlydecreasing part is due to increased magnetic support towards theouter disk (see Sect. 4.4.2 for discussions of the magnetic field). r(AU)10 ( k g / m ) M = 0.357( M ) M d = 0.0059 + 0.0008( M ) cyl r r(AU)10 ( k g / m ) M = 0.356( M ) M d = 0.0055 + 0.0016( M ) cyl r r(AU)10 M d ( M ) r(AU)10 ( k g / m ) M = 0.366( M ) M d = 0.0083 + 0.0019( M ) cyl r r(AU)10 ( k g / m ) M = 0.366( M ) M d = 0.0077 + 0.0024( M ) cyl r R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 5.
Disk surface density profile Σ cyl ( r ) measured inside a selectedcylindrical region (blue). Top two panels: R_ (cid:96)
14 and R_40ky_ (cid:96)
18 at103 kyr.
Middle: radially accumulated mass of the disk of R_40ky_ (cid:96) . Bottom two panels: R_ (cid:96)
14 and R_80ky_ (cid:96)
18 at 138 kyr.A 3-segment powerlaw fit is over-plotted (orange), and the powerlawexponents are shown in the legend. Two vertical dotted lines show thecharacteristic radii R kep and R mag , that correspond to the transition of thepowerlaw slope. The mass of the star M ∗ and that of the disk M d (inside R kep plus between R kep and R mag ) are also displayed. Most of the diskmass is contained within R kep .Article number, page 8 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core D i s k r a d i u s ( A U ) R Keplerian R magnetized R_40ky_l18R_80ky_l18
Fig. 6.
Radius evolution of the disk. The radii measured from simula-tions as described in Sect. 4.1.3 are plotted with blue solid ( R kep ) anddashed ( R mag ) lines for R_ (cid:96)
14 in blue. The disk has initially a Keplerianradius of about 20 AU for almost 50 kyr, and a magnetized envelopewhich is roughly 1.5 times larger. At later times, the disk experiences agrowth to nearly 50 AU in radius, and the magnetized region grows evenmore significantly to more than twice the size of the demagnetized Ke-plerian disk. Runs R_40ky_ (cid:96)
18 and R_80ky_ (cid:96)
18 are over-plotted (withsame line styles), and the size of the disk is consistent with di ff erentresolutions. Most of the disk mass is within the inner region, while some-times the region with sharply decreasing surface density can con-tain up to a quarter of the total disk mass. The outer flat part,on the other hand, is not a real surface density, but a truncatedvalue of the di ff use envelope that is more vertically extended.We perform a 3-segment powerlaw fit to the measured Σ cyl ( r ),and the transition points define the two characteristic radii, R kep and R mag , which will be discussed in the following paragraph. R Two characteristic radii, R kep and R mag , are defined from a 3-segment powerlaw fit of the surface density profile (Sect. 4.1.2).This is shown in Fig. 5 with the vertical dotted lines. The innerpart of the disk has a shallow powerlaw profile and is almostnon-magnetized. The outer part, on the other hand, is supportedby the magnetic field and the surface density drops quickly (seeSect. 4.4.2). These two characteristic radii will be marked in allthe following radial profile figures to show that they indeed cor-respond to some transitions of the properties inside the disk.It is to be kept in mind that the there does not exist a clearcutdisk radius, while our choice of characteristic values allows mak-ing reasonable evaluations of other disk properties in our analy-ses. With the non-ideal MHD e ff ects, the disk becomes self-regulatedand its size evolves slowly. This is due to the balance be-tween the rotation and the magnetic braking as well as that be-tween the magnetic field induction through rotation and its lossthrough di ff usion (Hennebelle et al. 2016). The radius reachesvery quickly ∼
20 AU and becomes quasi-stationary. The diskreceives mass from the infalling envelope, and at the same timesome of its mass is lost through accretion onto the central star.The mass inside the protoplanetary disk stays at ∼ . M (cid:12) .Figure 6 shows the disk radius evolution for the canonicalsimulation at 1AU resolution, with R kep plotted in solid blue lineand R mag in dashed blue line. The sink particle forms at simu-lation time 61 kyr, and we only show measurements from ∼ M a ss ( M ) MM Keplerian M magnetized R_40ky_l18R_80ky_l18
Fig. 7.
Mass evolution of the star-disk system. Run R_ (cid:96)
14 is plottedin blue, thick line for the sink mass, thin line for the mass within R kep and dashed line for the mass within R mag . A flattened core-disk struc-ture forms at the beginning. Once the sink forms, it quickly accretes themajority of the mass inside the disk, and the disk mass drops to aboutone percent that of the star. The magnetized region around the Keple-rian disk, though significant in size, contains much less mass. The highresolution restarts are over-plotted with same line styles, and there is anoverall good agreement. sink particle and a flattened disk is well defined. As can be seenin Fig. 1, the star-disk system is more like a flattened ellipsoidat 61 kyr when the sink particle has just formed. The size, R kep in particular, does not vary strongly in a time span of more than50 kyr. At around 120 kyr, the disk starts to grow in size up to ∼
50 AU. The magnetized part grows significantly, reaching acouple of times the size of the Keplerian disk. The growth of thedisk is possibly a consequence of the evolution of the surround-ing medium properties and, in particular, of the magnetic field.It may also indicate that the coupling with the surrounding en-velope, with which angular momentum is exchanged, becomesless tight, and thus the e ffi ciency of magnetic braking is reduced.This particular point certainly requires further investigation. Thehigh resolution restarts are over-plotted and the disk sizes areconsistent. With increased resolution, the time steps are signif-icantly reduced such that the evolution is not followed duringvery long time spans.Figure 7 shows the mass evolution of the star (thick blueline), the disk mass within R kep (thin blue line), and that within R mag (dashed blue line). Before the sink formation, a flatteneddense structure is formed. Quickly after the sink particle for-mation, most of the mass is accreted by the star and the diskmass drops to stably ∼ −
3% of the stellar mass. The majorityof the mass is contained within R kep , where the disk is essen-tially Keplerian. When the disk starts to expand at later times,the magnetized region grows significantly and its mass can reachsometimes a third of that of the Keplerian part. The over-plottedhigh-resolution runs are essentially identical.Figure 8 shows the mass accretion rate of the sink particleand that across the disk surface (at min(3 h p , dx )). The accre-tion rate of the sink is directly derived from the sink mass history.There are some bursty behaviors, which is probably due to themass accumulating at the inner disk boundary before the densityexceeds n acc . Ku ff meier et al. (2018) also saw similar behaviorin their zoomed-in systems at 2 AU resolution, while such fea-ture diminishes with increased resolution. Therefore such phe-nomenon should be interpreted with care. The two values arebroadly comparable. The mass accretion onto the disk surface isslowly decreasing in time. On the other hand, the sink accretionrate experiences a more serious drop. This is probably becausethe accretion becomes numerically more di ffi cult as the disk den-sity slowly decreases while n acc is kept a constant. At later times Article number, page 9 of 27 & A proofs: manuscript no. disk_formation_2020Feb M ( M / y r ) MM Disk , 3 H R_40ky_l18R_80ky_l18
Fig. 8.
Mass accretion rate evolution. The canonical run at 1 AU resolu-tion is plotted in blue. The mass accretion rate of the sink particle (thickline) is derived directly from its mass evolution. The accretion rate ontothe disk is the integral of the measured source function (dashed line,as described in Sect. 4.5.3). The mass accretion rate is initially highan decreases slowly with time. The mass accretion onto the sink parti-cle, though generally decreasing in time, shows some bursty behavior.This is due to the accumulation at the inner disk border shortly beforereaching the accretion threshold. After 130 kyr, the density of the diskdrops so much such that the threshold density of the sink accretion algo-rithm is no longer met and the sink hardly accretes. The high resolutionrestarts are over-plotted (solid line for sink accretion and stars for disksurface accretion). The mass flux onto the disk surface is in generalunchanged, while the sink accretion is sensitive to the inner boundarycondition of the disk (i.e., the choice of n acc , see Sect. 2.3). These resultsactually suggest that the threshold is too high in the restarts at 40 kyr,giving zero sink accretion, and slightly too low in the 80 kyr restart. ˙ M Disk > ˙ M sink . This is also reflected in the growth of mass andsize of the disk (Figs. 6 and 7).The mass accretion rate across the disk surface in the high-resolution runs (marked with stars) is compatible with the canon-ical run. However, R_40ky_ (cid:96)
18 presents zero sink accretion.This reflects the sensitivity of sink evolution to the numericalparameters, and n acc was probably not very well adjusted forthis case (see Sect. 2.3). On the other hand, R_80ky_ (cid:96)
18 showshigher sink accretion rate (solid line) than the canonical run. Thismight be due to the better description of the central high-densityregion, or it might as well suggest that the value of n acc is slightlytoo low. One of the important properties of the protoplanetary disk is thetemperature distribution, as it is important for the condensationof di ff erent minerals, as well as water. Knowing the tempera-ture conditions will allow to study the formation of precursorsof di ff erent mineral components of the building blocks of the S.S.. The water content, on the other hand, plays crucial roles notonly in the conditions of chemical reactions but also but also forthe growth of solid particles, since water ice represent the mostmassive reservoir of solid material in the S. S.. In this session,we discuss the thermal structure of the disk and ts evolution. Figure 9 shows snapshots of radial temperature profile. The pro-files averaged (mass-weighted) within one and three times thedisk scale height H are shown, and they are indeed almost iden-tical, confirming the vertical isothermal simplification. Note thatthere is no properly defined mid-plane temperature, because we r(AU)10 T ( K ) H H r(AU)10 T ( K ) H H r(AU)10 T ( K ) H H r(AU)10 T ( K ) H H R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 9.
Radial temperature profiles of the disk. Same snapshots as inFig. 4 are shown. The temperature is averaged (mass-weighted) ver-tically within H (orange) and 3 H (green). The snow-line ( R ) andthe CAI condensation line ( R , if present) are marked with verticalcyan dashed lines. In the low resolution run, the snow-line R is veryclose to the resolution limit and thus not well resolved. High resolutionrestarts are indeed necessary for studying the interior structure of thedisk correctly. The CAI condensation line R appears in R_40ky_ (cid:96) (cid:96)
18 this line has migrated inward due to the decreasedoverall density and temperature. only have access to the temperature at the resolution limit, thatis, within one cell height. The value of h p is used for H , thus itis more reliable within R kep . Meanwhile, the gas is relatively dif-fuse outside and vertical temperature variation is considerablysmall. As a verification for the thermal sound speed evaluatedin Sect. 4.1.1, we examined the vertical temperature profiles.The temperature only deviates significantly at rather high alti-tude compared to the disk scale height, and we refer readers toAppendix B for more details. Article number, page 10 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core
In general, the temperature decreases with increasing dis-tance to the central star. This is mostly due to the heating fromthe central star, since the temperature of the disk is significantlylower when the stellar irradiation is not considered (see Figs. 10and 11). We caution that the kink seen near 4 AU in the canonicalrun is not physical. It actually corresponds to the sink accretionradius and where gas behavior is no longer physical. The tem-perature profile is close to what would be expected for a flareddisk illuminated by the star, where T ∝ r − / (Kenyon & Hart-mann 1987). Such temperature profile is derived by balancingthe flux received at the disk surface ( ∝ r − ) to the local blackbody radiation ( ∝ T ). More specifically, Chiang & Goldreich(1997) derived T ∝ r − / α / , (10)where α g = dh p / dr − h p / r is the grazing angle at which the stel-lar irradiation strikes the disk. Combined with our measurementsof h p ∝ r . − . (see Fig. 4), one can derive a temperature profileslightly shallower than r − / . However, our disk is surrounded byan envelope which is not optically thin, and the geometry of thedisk is not very flared (see Fig. 4). The photons emitted from theproto-star is unlikely able to penetrate freely through the upperlayers and strike the disk surface at large radii directly. It is thusmore appropriate to describe the temperature distribution in thedi ff usion regime, where the radial temperature gradient is estab-lished through the equilibrium of radiative di ff usion, such that(cf. Eq. (5) of Hennebelle et al. 2020b) − π r σ κ ( T ) ρ ( r ) ∂ T ∂ r = L ∗ , (11)where L ∗ is the luminosity from the proto-star, κ the opacity,and σ the Stefan-Boltzmann constant. The density profile ρ ( r ) isslightly shallower than r − in general (see Fig. 3 of Hennebelleet al. 2020a). For roughly constant opacity (Semenov et al. 2003)and ρ ∼ r − ∼− , one can derive T ∝ r − / ∼− / . In conclusion, thedisk temperature profile is similar to what would be expectedfrom a irradiated disk, while the physical origin of such profileis actually di ff erent.For a passive disk that is heated only by viscous dissipation,Bitsch et al. (2013) derived T ∝ r − / − s / (constant viscosity) or T ∝ r − / − s / ( α viscosity), where s is the powerlaw exponentof the surface density profile Σ ∝ r − s , depending on di ff erentviscosity models. Again, these profiles are very similar to theabove-derived ones. Therefore, the temperature profile is not avery sensitive proxy for distinguishing among di ff erent heatingmechanisms of the disk.Here we are particularly interested in two temperature val-ues: the water snow-line at 150 K and the CAI condensationtemperature at 1500 K. We use 1500 K as a representative value,while it should be kept in mind that the CAIs are a family ofrefractory minerals and form in a temperature range around thisvalue. We identify the two radii, R ∼ R < dx = .
12 AU) that R is marginally resolved. Figure 11 shows the evolution of R , the water snow-line. Wedo not discuss R since it is only marginally resolved in theruns with the highest resolution. R starts at around 10 AU and decreases slowly in time, asthe global density of the disk decreases due to accretion onto thestar and redistribution of material in the envelope surrounding r(AU)10 T ( K ) H H R_ (cid:96)
14 vs. R_ (cid:96)
Fig. 10.
Radial temperature profiles of the disk with and without stellarirradiation at time 120 kyr. Values of R_ (cid:96) H (orange) and 3 H (green). The snow-line ( R ) is marked withvertical cyan dashed line. The disk temperature is significantly lowerwhen no stellar irradiation is introduced, while the slope of radial profiledoes not seem to di ff er significantly. S n o w li n e ( A U )
150 KR_40ky_l18R_80ky_l18R_l14_nofeed
Fig. 11.
Evolution of the water snow-line, R . The snow-line situ-ates initially around 10 AU and migrates inward in time as the diskevolves and loses its mass, finally stabilizing at ∼ (cid:96)
18 has consistent temperature, while inR_40ky_ (cid:96)
18 the snow-line migrates to larger radius, possibly due to themass accumulation in the disk and the increase in opacity. This reflectsthe importance of properly choosing n acc with increased resolution (seeSect. 2.3). The run with no stellar irradiation is also plotted for compar-ison, where the temperature is significantly lower. the disk. After about 80 kyrs the snow-line stabilizes at about5 AU, close to the present location of Jupiter, which may haveconsequences for its formation (see Sect. 5). R from high-resolution restarts are over-plotted. Run R_80ky_ (cid:96)
18 shows con-sistent snow-line position. On the other hand, the R_40ky_ (cid:96) n acc is not entirely con-sistent when we increase the resolution (see Sect. 2.3). Recallthat the sink does not accrete in this high-resolution restart (Sect.4.2). The slight increase in surface density, which is barely per-ceptible if one looks at the Σ profiles, leads to increased opacityand thus increased temperature. Due to the same issue, the evo-lution trend of the snow-line we obtained should be regarded assuggestive and its exact position depends on the detailed physicsof the stellar accretion.To evaluate the significance of the stellar irradiation, we per-formed run R_ (cid:96) ff at simulation time 119 kyr. As shownin Fig. 11, the temperature of the disk turns out to be signifi-cantly lower right away. This implies that in our simulations, thedisk temperature is primarily set by the stellar irradiation. How- Article number, page 11 of 27 & A proofs: manuscript no. disk_formation_2020Feb r(AU)10 u ( m / s ) u kep u ( H ) u (3 H ) u r ( H ) u r (3 H ) r(AU)10 u ( m / s ) u kep u ( H ) u (3 H ) u r ( H ) u r (3 H ) R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)
Fig. 12.
Velocity profiles of the disk at two resolutions: rotational veloc-ity within H (orange) / H (green) and infall velocity within H (red) / H (violet). Both positive and negative values with | u | > are plotted inlogarithmic scale, and this is connected with a linear range in between.Keplerian rotational velocity is plotted in blue for reference. The innerdisk is very close to Keplerian rotation, and the deviation becomes re-markable in the outer magnetized part of the disk. The envelope is farfrom Keplerian rotation, while the infalling gas has roughly half thefree-fall velocity. Below the disk accretion radius (vertical gray line),the velocity is not physically meaningful. The disk characteristic radii R kep and R mag (the dotted vertical lines) mark the transition between dif-ferent behaviors. ever, since this run is only performed for a short time span, it isunclear whether the disk will reach a di ff erent equilibrium statewith slightly higher temperature by increasing the density. Thee ff ect of irradiation should be investigated in detail, along withthe luminosity resulting from mass accretion, while we leave itfor future studies. Some disk properties, such as the growth and radial drift ofdust and pebbles, as well as planetesimal formation, are closelylinked to its dynamics and evolution. As it is always assumedthat the protoplanetary disk has close-to-Keplerian rotation, wefirst check the validity of this common assumption. We then con-tinue to discuss the magnetization and the transport parametersthat are essential in controlling dust-growth and transport withinthe disk.The azimuthally averaged quantities are presented, and weshow both the average within one and three times the disk scaleheight H . Since respectively 68.3 % and 99.7 % of mass arecontained within these extents for a vertical Gaussian profile,these quantities should be representative of the disk. Note thatall quantities are calculated throughout the selected radius, r cyl ,while the notion of scale height is only valid within the diskradius R kep . The quantities inferred in the envelope region aretherefore only suggestive of the values around mid-plane. We show the rotation velocity, u φ , and the radial velocity (par-allel to the disk plane), u r , in Fig. 12. The Keplerian rotationvelocity is plotted for reference: u kep ( r ) ≈ (cid:114) G [ M ∗ + M d ( r )] r . (12)This approximation is valid for disk mass dominated by that ofthe star.The mass weighted average is calculated within one andthree disk scale heights. The azimuthal velocity is plotted in or-ange ( H ) and green (3 H ). Di ff erent behaviors are clearly seen inregions separated by the disk characteristic radii R kep and R mag (see definition in Sect. 4.1.2 and here we see the reason for giv-ing those names). The velocity is very close to the Keplerianrotation between ∼ ∼
13 AU( R kep ). Inside the sink accretion radius, the flow property is notphysical and should not be considered. Beyond R kep , the rota-tion starts to deviate from the Keplerian value since significantsupport comes from the magnetic field. Beyond R mag it is clearthat the gas inside the infalling envelope is far from Keplerianrotation.On the other hand, the infall velocity (negative radial veloc-ity) is plotted in red ( H ) and purple (3 H ). In the envelope, thegas falls radially on to the disk at almost half the free-fall veloc-ity, which has the same expression as the Keplerian velocity. Inthe interior of the disk, the radial velocity is significantly lowerthan the rotational velocity (by at least two orders of magnitude).The radial velocity flips between inward and outward directionswith no clear pattern, and also varies with time. The di ff erencebetween the two runs at di ff erent resolutions indicates that theinternal dynamics is highly variant and has not yet reached nu-merical convergence. Nonetheless, the innermost disk (below 1AU) always shows accretion towards the central star. As will bediscussed in Sect. 4.5.1, most of the mass accretion onto the staroccurs directly through a high altitudes channel and does notpass through the bulk of the disk. Therefore, the disk is not nec-essarily always in accretion. β parameter: magnetization The ratio between the thermal and magnetic pressures, β = P therm / P mag = πρ c / B , indicates which supporting agent isdominant inside and outside the disk region. This is calculatedas β ( r ) = π (cid:80) PdV (cid:80) B dV , (13)Figure 13 shows β values within H and 3 H . Due to the hightemperature inside the disk, the thermal pressure largely domi-nates with β decreasing outwards from 10 at 1 AU. The value of β crosses unity inside the magnetized zone of the disk (between R kep and R mag ), marking a transition to the magnetized infallingenvelope. As the global collapse proceeds, the central high den-sity part loses the field lines through ambipolar di ff usion and themagnetic pressure is therefore weak in the inner part of the disk.The disk is more magnetized at higher altitudes, giving lower β values at 3 H than at H (since the temperature is almost identicalfrom Fig. 9). Although significantly demagnetized, the magneticfield is still stronger than what is typically assumed for the stud-ies of disk instabilities. For example, Béthune & Latter (2020)used β ∈ [10 , ]. Article number, page 12 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core r(AU)10 H H r(AU)10 H H r(AU)10 H H r(AU)10 H H R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (138 kyr)R_ (cid:96)
14 (103 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 13.
Plasma beta, the ratio between thermal and magnetic pressures, β = P therm / P mag , evaluated within one (orange) and three (green) timesthe disk scale height. The disk is highly demagnetized due to the am-bipolar di ff usion that leaves the field lines outside during the collapse.The β value passes through unity around the junction between the diskand the envelope. The di ff erence between the two curves comes fromthe lower magnetization level close to the disk mid-plane. α parameter: accretion and diffusion The α parameter describes the transport in the disk as a resultof torques from the turbulence or magnetic field. The cross-correlation of di ff erent components of a vector field generate atensor term that implies the transportation of mass and angu-lar momentum (Balbus & Papaloizou 1999). It has contributionsfrom the turbulent fluctuation (Reynolds tensor), the magneticfield (Maxwell tensor), and the self-gravity of the disk. The tur-bulent di ff usion is not very well known within the disk, and thus ν = α Rey c / Ω is often used as the turbulent viscosity for describ-ing the di ff usion within the disk. We evaluate radial dependence of these terms inside the disk and discuss their e ff ect on the trans-port of mass and angular momentum.The stress tensor, T p φ ( p , φ ) ≡ (cid:104) ρδ u p δ u φ (cid:105) φ , is an azimuthalaverage along a circle perpendicular to the axis, where the sub-script p stands for the poloidal direction that has two compo-nents r and z and δ denotes the deviation from the mean value.The dimensionless α expression is defined as T p φ / (cid:104) ρ c (cid:105) φ , wherethe denominator is the azimuthally averaged thermal pressure.The radial mass accretion rate of a stationary disk is linkedto the stress tensor with the transport equation (See details inAppendix C)˙ M ( r ) = − π r Σ (cid:104) u R (cid:105) r (14) = π ( r Ω ) (cid:48) (cid:34) ddr (cid:32) r (cid:90) H − H T r φ dz (cid:33) + r T z φ | H − H (cid:35) = π ( r Ω ) (cid:48) (cid:34) ( r Σ c α r ) (cid:48) + r H Σ c α z (cid:35) ≈ (cid:114) π Σ rG c (cid:18) α r + r H α z (cid:19) , where (cid:48) signifies the derivative with respect to r . The approx-imation is obtained by making several assumptions such that d / dr ≈ / r and Ω ≈ G (cid:82) π r Σ dr ≈ π G Σ r , which are not toofar from reality (see e.g. Fromang et al. 2013), with the purposeto give an appreciation of how the two terms in the parenthesiscompare to one another. More precisely, α r is measured in thebulk of the disk, while α z is measured on the surfaces, and theire ff ects on the mass accretion rate di ff er by a factor r / (2 H ). Thesetwo values should therefore not be compared directly.While the stress tensor is defined along a thin circular line, inpractice the average is taken within a finite volume. The α valuesare thus calculated as α Rey , r = ˆ T Rey , r φ / ˆ P , α Rey , z = ( ˆ T H Rey , z φ − ˆ T − H Rey , z φ ) / ˆ P , (15) α Max , r = ˆ T Max , r φ / ˆ P , α Max , z = ( ˆ T H Max , z φ − ˆ T − H Max , z φ ) / ˆ P , (16) α Grav , r = ˆ T Grav , r φ / ˆ P , α Grav , z = ( ˆ T H Grav , z φ − ˆ T − H Grav , z φ ) / ˆ P , (17)where the hut symbol, ˆ , signifies the volume average. The re-gion of the summations is an annulus of radial thickness dx ,while its vertical thickness is to be specified.The volume averaged stress tensors and pressure are definedas:ˆ T p φ, Rey (cid:88) dV = (cid:88) δ u p δ u φ ρ dV = (cid:88) δ u p δ u φ dm (18)ˆ T p φ, Max (cid:88) dV = π (cid:88) B p B φ dV = (cid:88) u A , p u A ,φ dm , (19)ˆ T p φ, Grav (cid:88) dV = π G (cid:88) g p g φ dV = (cid:88) u G , p u G ,φ dm , (20)ˆ P (cid:88) dV = (cid:88) P therm dV = (cid:88) c dm , (21)where u A = B (cid:112) πρ , (22) u G = ∇ Φ (cid:112) π G ρ = − g (cid:112) π G ρ . (23) dV and dm represent the volume and mass of each cell.The three radial terms α Rey , r , α Max , r , and α Grav , r are relatedto the radial transport within the disk. Their values are there-fore calculated as the mean across the disk vertical extent. The Article number, page 13 of 27 & A proofs: manuscript no. disk_formation_2020Feb r ( AU )10 r , d x Rey ( H ) Rey (3 H ) Max ( H ) Max (3 H ) r ( AU )10 z , d x Rey ( H ) Rey (3 H ) Max ( H ) Max (3 H ) Fig. 14.
Disk α for R_40ky_ (cid:96)
18 at 103 kyr. Top panel presents theReynolds and Maxwell stresses in r direction, and bottom panel thetwo components in z direction. Both positive and negative values with | α | > − are plotted in logarithmic scale, and this is connected with alinear range in between. summations are made within H or 3 H . The total α r is the sumof all the three components. On the other hand, α Rey , z , α Max , z ,and α Grav , z are the vertical counterparts that are linked to the ac-cretion or angular momentum transport across the disk surfaces.The sum of the three terms, α z , also leads to radial mass accre-tion. The summation is performed within a vertical thickness of dx around H and 3 H in this case for the stress tensors, while thepressure is still averaged within the disk vertical extent. With thedefinition in Eq. (14), it naturally justifies the choice of normal-ization by the bulk pressure, whether the stress tensor is averagedin the bulk ( T r φ ) or only on the surface ( T z φ ).Figure 14 shows the r -(top panel) and z -(lower panel) com-ponents of α . The gravitational term is mostly noise and we donot show it for conciseness (see Appendix D for more details).The solid lines are the values measured at H and dashed lines arethose measured at 3 H . For the α r component at H , the Maxwellterm at inner part of the disk is almost zero since the disk isstrongly demagnetized as a consequence of the ambipolar di ff u-sion (see also Sect. 4.4.2). It grows steadily when approachingthe magnetized outer part of the disk and reaches ∼ − at R kep .The Reynolds term fluctuates around zero with a maximal abso-lute value lower than 10 − all along the disk. Both componentsgrow significantly in the outer zone of the disk between R kep and R mag . The Reynolds component flips between positive and neg-ative values, showing a complex transport pattern in the innerdisk, which is also seen in the radial velocity (Sect. 4.4.1). Thevalues at H and 3 H are generally comparable, implying that theradial stress tensors does not cause significant stratified behav-ior within the vertical extent of the disk. On the other hand, thevalues of α z show more significant variation across di ff erent al-titudes, and this is the dominant component that leads to radialaccretion within the disk. The Maxwell term at 3 H is higher thatat H above r ∼ H , the Reynolds component is significantly negative, implying possibly a viscous outwardtransport near the disk mid-plane.Equation (14) implies that α z has a stronger e ff ect than α r byroughly a factor 5 in the case of aspect ratio 0.1 if α z (cid:39) α r . Withthe slightly higher value of α z compared to α r , one can infer thatmost of the radial mass accretion is due to the angular momen-tum transport along the field lines on the surface of the disk, andthis accretion happens mostly in the upper layers (between H and 3 H , or even above) of the disk.The mass accretion rate derived from α is shown in Fig. 25with dotted curves and compared to the directly measured valueslater. Large fluctuations are present, while the overall magnitudeis comparable with the directly measured values of radial accre-tion (see Sect. 4.5.4).The above measurements of α along given surfaces reveal animportant information, that is, the transport onto and within thedisk is highly structured. To give a better idea of the transportpattern and the behavior of di ff erent stress tensors, we show α values on the r − z plane for several snapshots of the canoni-cal run in Fig. 15. We also examine the e ff ect of high resolutionrestarts in Figs. 16 and 17. Positive values are shown in red, cor-responding to radial accretion (or angular momentum removal,see Eq. (14)), and negative values in blue signify radial expan-sion. The total e ff ect has contribution from all the four termsshown in the figure.Figure 15 shows that at early times α Rey usually dominateswithin the disk but does not show regular pattern (first two snap-shots). This is probably linked to the strong accretion that con-stantly perturbs the disk, and it is to be kept in mind that the disksize fluctuates in time under the e ff ect of irregular infall of theenvelope. Moreover, the disk size is very small at the beginning(10 AU for one AU resolution), and a good determination of theinterior structure of the disk is in principle not possible. At latertimes, the disk grow to a few tens of AU in size, giving moreresolved internal structure.Indeed, the middle layers of the disk does not have very wellbehaved values of α and some parts of the disk mid-plane canactually expand. Most of the accretion is channeled through ahigh-altitude layer of the disk, moving radially inwards beforereaching the surface of the disk. This can be seen most clearly inthe figures of α Rey , z and α Max , z as a dark red region slightly abovethe line of 3 H .Figures 16 and 17 show the comparisons with the canonicalrun of R_40ky_ (cid:96)
18 and R_80ky_ (cid:96)
18. At 40 kyr, the restart withincreased resolution allows us to study the internal structure ofthe disk without having to evolve the whole simulation at highresolution since the beginning. The di ff erence in α Rey , r is moreremarkable but it is not the dominant term. Accretion happensmostly at high altitude (see α Max , z ) and the middle layers of thedisk is highly structured. At 80 kyr, the disk grows much larger,and thus even the canonical run allows satisfying evaluation of α . The similarity of the α patterns between the two runs with afactor 16 resolution di ff erence is a validation of numerical con-vergence. Moreover, the high resolution restart also allows usto probe much smaller radii and show that the similar disk andaccretion structure continue as we go to smaller radii. Q parameter: stability and fragmentation The Toomre parameter, Q , indicates whether the disk is subjectto gravitational fragmentation. It is defined as Q ( r ) = κ ( r ) c s ( r ) π G Σ ( r ) ≈ Ω ( r ) c s ( r ) π G Σ ( r ) , (24) Article number, page 14 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core
Rey , r ( t =110 kyr , M =0.36 M ) Rey , r ( t =118 kyr , M =0.36 M ) Rey , r ( t =126 kyr , M =0.37 M ) Rey , r ( t =134 kyr , M =0.37 M ) Max , r ( t =110 kyr , M =0.36 M ) Max , r ( t =118 kyr , M =0.36 M ) Max , r ( t =126 kyr , M =0.37 M ) Max , r ( t =134 kyr , M =0.37 M ) Rey , z ( t =110 kyr , M =0.36 M ) Rey , z ( t =118 kyr , M =0.36 M ) Rey , z ( t =126 kyr , M =0.37 M ) Rey , z ( t =134 kyr , M =0.37 M ) Max , z ( t =110 kyr , M =0.36 M ) Max , z ( t =118 kyr , M =0.36 M ) Max , z ( t =126 kyr , M =0.37 M ) Max , z ( t =134 kyr , M =0.37 M ) Fig. 15. α evolution of R_ (cid:96)
14 ( increasing time from left to right ). The four components α Rey , r , α Max , r , α Rey , z , α Max , z ( from top to bottom ). The r components are vertically averaged values between zero and given altitude, and the z components are local values. The white curves trace thethermal scale height H and 3 H . Positive values ( red ) can be roughly interpreted as a radial accretion and negative values ( blue ) for radial expansion.The resolution is one AU. As the disk grows in size, the interior of the disk is described with more cells and there is a clearer pattern of α . Theearly disk has relatively strong α Rey compared to its magnetic counter part. It is probably under constant perturbation by the important infall fromthe envelope. The first two chosen snapshots are not necessarily representative of this epoch since the behavior of α Rey is very irregular in time,leading to a disk that fluctuates in size around 10 AU. Later on, the infall from the envelope weakens, and the disk becomes more regular, with itsinner parts readily accreting. This accretion is mostly caused by vertical magnetic transport of angular momentum within the envelope. where κ = r d ( r Ω ) dr = Ω r d ( r Ω ) dr (25)is the epicyclic frequency. The approximation, κ ≈ Ω , is validfor a disk in Keplerian rotation, therefore valid only within ∼ R mag . The Toomre parameter compares the stabilizing radialshear from di ff erential rotation to its self-gravity that induces lo-cal collapse. For a value of Q >
1, the disk is stable againstself-gravitating fragmentation.Figure 18 shows the Toomre Q profile for some snapshots.The value is very high near the center due to the hight tempera-ture, and decreases to reach a lowest point near R kep . The valueof Q increases again in the exterior magnetized part of the disk,while the approximation in Eq. (24) starts to be invalid. Withthe large value of Q , the whole disk is very stable against self-gravitating fragmentation. All values of c s , Ω , and Σ are evalu-ated locally within H or 3 H . The di ff erence in the vertical direc-tion is mostly because all mass is not taken into account whenonly the region within H is considered.The Toomre Parameter Q ∼
10 near the edge of the Ke-plerian disk, and it reaches very high values at smaller radii(10 − ). Such extreme stability, from high temperature and / orlow surface density, should be regarded with care and we pointout some possible origins. First, the disk could be potentiallyover-heated with our prescription of proto-stellar irradiation (seeSect. 4.3 and 5). Secondly, the disk global profile is likely tobe governed by the star-disk boundary condition, or the sink ac-cretion algorithm from a numerical aspect. Our simulation haspercent-level disk-star mass ratio, which is somehow low com-pared to some other studies (e.g. 30-40% was found by Tomidaet al. 2017). Such result is directly linked to the density thresh- old for the sink accretion (see Sect. 2.3). This has a direct conse-quence on the surface density of the disk, and therefore large Q values could occur if Σ is too low. The primary goal of this study is to follow the building phase ofthe protoplanetary disk: how mass arrives onto the disk and howmass is transported within the disk. We therefore start by lookingat some general properties of the global flow, and then continuewith more detailed quantification of the mass flow properties.
Figure 19 shows R_ (cid:96)
14 (top) and R_40ky_ (cid:96)
18 (bottom) at103 kyr. Colored streamlines present the azimuthally averagedpoloidal mass flux, (cid:104) ρ u p (cid:105) φ , on top of the specific angular mo-mentum, (cid:104) ρ ru φ (cid:105) φ / (cid:104) ρ (cid:105) φ , shown in gray scale. White curves markthe disk scale height H and 3 H up to R kep . The general behav-ior is not strongly a ff ected by the resolution. The flow patternvaries slightly at di ff erent instants. Nonetheless, they all showsome similar characteristics that we present below: infall, out-flow, shock, and accretion (or sometimes expansion / decretion).Figure 20 presents a schematic view of the various zones. – Infall: The envelope is globally infalling. The infall is mostlyradial near the equatorial plane. Slightly above the disksurface, the vertical mass flux becomes more significant.Most importantly, the flux intensity varies (change of color)with the altitude due to the small radial component thatbrings the mass flux to smaller radii as it approaches the disk
Article number, page 15 of 27 & A proofs: manuscript no. disk_formation_2020Feb
Rey , r ( t =103 kyr , M =0.36 M ) Rey , r ( t =103 kyr , M =0.36 M ) Max , r ( t =103 kyr , M =0.36 M ) Max , r ( t =103 kyr , M =0.36 M ) Rey , z ( t =103 kyr , M =0.36 M ) Rey , z ( t =103 kyr , M =0.36 M ) Max , z ( t =103 kyr , M =0.36 M ) Max , z ( t =103 kyr , M =0.36 M ) Fig. 16.
Measured α values of R_ (cid:96)
14 ( left ) and R_40ky_ (cid:96)
18 ( right )at 103 kyr. The disk radius is about 15 AU. The bins for averagingare chosen to be larger than dx in the high resolution run in order toget rid of numerical noises. These figures clearly show that the tech-nique of restart using higher resolutions is needed to study the interiorproperties of the disk. Note that the z components are more importantsince they dominate in governing the mass accretion inside the disk.The middle layers of the disk are more irregular with some partially ex-panding parts, while most of the accretion happens through the upperlayer slightly above 3 H and reaches small radius before penetrating thesurface of the disk. Rey , r ( t =138 kyr , M =0.37 M ) Rey , r ( t =138 kyr , M =0.37 M ) Max , r ( t =138 kyr , M =0.37 M ) Max , r ( t =138 kyr , M =0.37 M ) Rey , z ( t =138 kyr , M =0.37 M ) Rey , z ( t =138 kyr , M =0.37 M ) Max , z ( t =138 kyr , M =0.37 M ) Max , z ( t =138 kyr , M =0.37 M ) Fig. 17.
Measured α values of R_ (cid:96) (left) and R_80ky_ (cid:96) (right) at138 kyr. Disk radius ( R kep ) grows to 50 AU at this time. The values areevaluated on grid of 1 AU resolution for both runs. The global patternsare compatible for the two runs with 16 times di ff erence in resolution.Disk middle layers are radially structured and partially expending attimes, while the accretion happens through a channel slightly above 3 H and brings significant amount of mass from the infalling envelope tosmall radius without passing though the disk. r(AU)10 Q H H r(AU)10 Q H H r(AU)10 Q H H r(AU)10 Q H H R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 18.
Toomre Q ≈ Ω c s / ( π G Σ ) parameter evaluated within one (or-ange) and three (green) times the disk scale heights. Same snapshots aretaken as in Fig. 4. The disk is very stable against self-gravitating frag-mentation. The di ff erence between the two curves comes mainly fromthe total mass not being taken into account for the value within H and3 H . surface from the top. Therefore, most of the infalling massreaches small radius of within a few AU before penetratingthe disk surface. Let us recall the source function, which isthe mass flux crossing the disk surface, that we aimed tomeasure. This kind of infalling pattern gives a surface fluxthat is centrally concentrated and strongly depends on thereference surface (altitude). – Outflow: In the axial direction, there is a low densityoutflow cavity that is launched at 10-20 AU above the diskmid-plane, and the radial extent increases with the verticaldistance from the disk. The transition between infall andoutflow happens at high altitudes above the disk, roughly
Article number, page 16 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core u p ( ru ) u p ( ru ) Fig. 19.
Poloidal mass flux ( (cid:104) ρ u p (cid:105) φ ) and specific angular momentum( r (cid:104) u φ (cid:105) φ ) around the disk region (color streamlines on gray image, inSI units), for R_ (cid:96)
14 ( top ) and R_40ky_ (cid:96)
18 ( bottom ) at simulation time103 kyr. All quantities are azimuthally averaged and averaged betweenthe upper and lower planes. White curves outline the disk scale height H and 3 H at r < R kep . An outflow is launched at high altitudes at z ≈
20 AU in R_ (cid:96)
14. In the high resolution run, the flow is slightlymore messy but qualitatively similar. The outflow is launched closer tothe disk ( z ≈
10 AU) and the disk size is slightly smaller. The outflowhas relatively low density and manifests itself as a cavity. Significantinfall from the envelope is channeled through a zone of low angularmomentum at several scale-heights (darker gray region near z ≈ r ) ofthe disk. The horizontal infall from the envelope reaches the outer edgeof the disk (as can be seen in Fig. 2 at r ≈
15 AU) and the densitycontrast generates a shock, that pushes the gas to higher altitude. Thiscreates a counterclockwise poloidal circulation near the disk edge be-tween 10 and 15 AU in R_40ky_ (cid:96)
18 since the outer disk is decreting atthis instant, while it is not seen in the upper panel.
Fig. 20.
Illustration for the zones showing major properties of the disk-envelope system in Sect. 4.5.1 and 4.5.2. near the cone z ≈ r . Right below this is a low-angular-momentum zone (darker on the gray image), through whichsignificant infall is channeled. There might be some circula-tion and recycling of the wind material across the cavity wall. – Shock: Recalling from Fig. 12, the radial velocity in theequatorial plane reaches a fraction of the free-fall velocityand quickly drops when arriving at ∼ R mag . In the post-shockregion, the rotation replaces the infall as the dominant ki-netic energy. A slowly expanding zone appears very often inthe magnetized outer region of the disk and meets the radialinfall at the shock front. This shock situates at r ∼
15 AU inboth panels of Fig. 19, where the radial velocity apparentlydiminishes. In R_ (cid:96)
14, it is more evident that the gas ispushed upwards due to the increased pressure and falls backto the disk surface at smaller radius, while the streamlinesin R_40ky_ (cid:96)
18 are more regular and there is only very mildupward motion when the gas reaches the shock front. Theradial flow inside the disk is sometime outwards, causinga circulation at the edge of the disk: matter goes upwardsat the shock and falls back from the top. This behaviorhas been suggested in observations (Sakai et al. 2017), andthe shock can be possibly traced with chemical molecules.However, in our case, the disk is probably too thin and theshock-generated heat is quickly radiated away. – Accretion: Within the disk ( z (cid:46) H ), the behavior of the flowis less clear. R_ (cid:96)
14 shows an almost vertical flux downwards,while in R_40ky_ (cid:96)
18 the flux is upwards. The disk internaldynamics is not numerically converged and conclusions can-not be drawn. Globally speaking, the innermost parts of thedisk (a few AU) is always in radial accretion, which leads tothe mass growth of the central star. However, the radial ve-locity within the rest of the disk is not strictly negative allthe time. Outward expansion sometimes occurs in the mid-dle layers of the disk. As mentioned earlier, this happens veryoften in the magnetized zone of the disk (between R kep and R mag ). However, the disk is confined by the magnetized zoneand the infalling envelope, and thus can not expand freely.At 103 kyr, the disk radius is roughly 15 AU, while it growsto 50 at 138 kyr (Fig. 21). The general behavior is the same as de-scribed above, and the low-angular-momentum region becomesa clear channel of mass infall. The expanding magnetized zone isvery large and almost comparable to the disk radius, R kep . Hereagain, the disk internal dynamics is not identical in the two runswith di ff erent resolutions. Nonetheless, the two major remarksstay valid: matter might circulate near the wind-infall wall andat the disk edge. The commonly accepted scenario describes the loss of angularmomentum due to various transport mechanisms in the envelope,which allows the matter to be accreted within the disk. The equa-tion of angular momentum transport is obtained by multiplyingthe φ -component of the momentum conservation equation by r and averaging azimuthally, which gives: ∂∂ t (cid:104) ρ ru φ (cid:105) φ + ∇ · (cid:34) (cid:104) r ρ u p u φ − r B p B φ π (cid:105) φ (cid:35) = . (26)The second term is the divergence of the poloidal flux of an-gular momentum. By separating the toroidal velocity into the Article number, page 17 of 27 & A proofs: manuscript no. disk_formation_2020Feb u p ( ru ) u p ( ru ) Fig. 21.
Poloidal mass flux ( (cid:104) ρ u p (cid:105) φ ) and specific angular momentum( r (cid:104) u φ (cid:105) φ ) around the disk region (color streamlines on gray image, in SIunits), for R_ (cid:96)
14 ( top ) and R_80ky_ (cid:96)
18 ( bottom ) at simulation time 138kyr. All quantities are azimuthally averaged and averaged between theupper and lower planes. The white curves outline the disk scale height H and 3 H at r < R kep . An outflow is launched at z ≈
20 AU. The diskhas grown in size while the flow pattern stay very similar to that in Fig.19. At this moment, the magnetized zone is almost twice the size ofthe Keplerian disk ( R mag ≈ R kep ), and the counterclockwise poloidalcirculation around R kep is evident in both panels and occupies almostthe totality of the magnetized zone and the outer Keplerian disk. circular and fluctuating components u φ = r Ω + δ u φ , where r Ω = (cid:104) ρ u φ (cid:105) φ / (cid:104) ρ (cid:105) φ , we can rewrite the angular momentum fluxas F j = (cid:104) r ρ u p u φ − r B p B φ π (cid:105) φ (27) = r Ω (cid:104) ρ u p (cid:105) φ + r (cid:104) ρδ u p δ u φ (cid:105) φ − r (cid:104) B p B φ (cid:105) φ π . The first term is the laminar transport, which simply describesthe the angular momentum flux associated with the mass fluxif the momentum is to be frozen within the mass. The sec-ond term is the turbulent transport terms that originated fromthe cross-correlation between the fluctuations. This is the majortransport mechanism in the α -disk model, that allows the diskto accrete. The last term is the magnetic transports associated tothe Maxwell stress tensor which tends to deviate the movementalong the field line direction, thus accelerating or braking therotation.In the case of an isolated disk at later phases, one is more fo-cused on the angular momentum within the disk and there is noimportant net mass flux. The tensor terms are therefore of ma-jor interests. In this study however, the disk is actively receiving r u u p ( ru ) r u u p ( ru ) rB B p ( ru ) r u u p + r u u p rB B p ( ru ) r u u p ( ru ) r u u p ( ru ) rB B p ( ru ) r u u p + r u u p rB B p ( ru ) Fig. 22.
Angular momentum flux and angular momentum ( r (cid:104) ρ u φ (cid:105) φ )around the disk region (color streamlines on gray image, in SI units),for R_ (cid:96)
14 ( left ) and R_40ky_ (cid:96)
18 ( right ) at 103 kyr. All quantitiesare azimuthally averaged and averaged between the upper and lowerplanes.
From top to bottom: laminar transport ( r (cid:104) u φ (cid:105) φ (cid:104) ρ u p (cid:105) φ ), turbu-lent transport ( r (cid:104) ρδ u φ δ u p (cid:105) φ ), magnetic transport ( − r (cid:104) B φ B p (cid:105) φ / (4 π )), andtheir sum. White curves outline the disk scale height H . The laminarflow brings in mass with positive angular momentum, while the tensorterms (Reynolds and Maxwell) transport the angular momentum out-wards. In particular, the magnetic tensor dominates in the envelope andthis allows the mass to fall onto the disk through the loss of angular mo-mentum. The general pattern persists with improved resolution. Whitepatches in the gray image is due to locally negative angular momentum.This is simply an artifact of the averaging, because the selected axisdoes not align exactly with the asymmetric outflow. mass from the infalling envelope and the laminar transport is notto be neglected. Figure 22 shows the laminar, turbulent, mag-netic, and total angular momentum fluxes (from top to bottom)at 103 kyr at two resolutions (left and right). From the pattern ofangular momentum transport, the disk-envelope system can bedivided into four zones (see Fig. 20 for an illustration). – The braked infalling layer: This region is dominated by thevertical upward transport of angular momentum by the mag-netic field, while the velocity is downwards. This means thatduring the infall, the gas is braked by the magnetic field thatis not aligned with the flow. This is the dark zone in Fig. 19where we see the specific angular momentum is the lowestand the mass is strongly infalling.
Article number, page 18 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core r u u p ( ru ) r u u p ( ru ) rB B p ( ru ) r u u p + r u u p rB B p ( ru ) r u u p ( ru ) r u u p ( ru ) rB B p ( ru ) r u u p + r u u p rB B p ( ru ) Fig. 23.
Angular momentum flux and angular momentum ( r (cid:104) ρ u φ (cid:105) φ )around the disk region (color streamlines on gray image, in SI units), forR_ (cid:96)
14 ( left ) and R_80ky_ (cid:96)
18 ( right ) at 138 kyr. The disk scale height H is shown in white. All quantities are azimuthally averaged and averagedbetween the upper and lower planes. From top to bottom: laminar trans-port ( r (cid:104) u φ (cid:105) φ (cid:104) ρ u p (cid:105) φ ), turbulent transport ( r (cid:104) ρδ u φ δ u p (cid:105) φ ), magnetic trans-port ( − r (cid:104) B φ B p (cid:105) φ / (4 π )), and their sum. – The free infalling layer: The region contained within a fewscale heights above the disk surface is dominated by thefree infall of mass. The turbulent and magnetic stress tensorstransport the angular momentum upwards, while this is dom-inated by the laminar flow that brings in directly the angularmomentum associated with the gas. In this region, the mag-netic field is mostly lost through ambipolar di ff usion. Thiszone is clearly seen by comparing Figs. 19 and 22 where themass flux and the total angular momentum flux are almost inthe same direction. – The wind cavity: This low density region is dominated bythe vertical upward transport of angular momentum by themagnetic field. Matter moves outwards along field lines andits angular momentum increases. – The disk interior: The disk interior (below a few H ) is dom-inated by the laminar transport of the angular momentum.This is linked to the global behavior of the disk and has tobe investigated with great care. While there are also radi-ally outward magnetic transport and vertical Reynolds trans-port, these two contributions are weak. We restrain ourselvesfrom over-interpreting the results since the low- and high-resolution runs are not consistent at this scale. Throughout the braked layer and the wind cavity, the magneticangular momentum flux is almost constant. This means that thefield lines channel out the angular momentum in an almost sta-tionary way and launches it away eventually through the wind.Figure 23 shows the same plots at 138 kyr, where the diskgrows to larger size. The general behavior is the same, whilewe have better resolved disk internal dynamics, especially inR_80ky_ (cid:96)
18. The disk mid-plane is under expansion at this mo-ment. It is clear that the angular momentum transport is largelydominated by the laminar transport that simply advects the an-gular momentum along with the mass.From the results shown above, the α -disk model is probablyunsuitable for an embedded disk, because the disk is not in a sta-tionary state and is very sensitive to irregularities in the infallingenvelope. The source function, S ( r ), in unit of mass per unit time per unitsurface, describes the rate and the spatial distribution of massarrival from the envelope onto the disk surface. Note that thesurface we refer to when discussing the mass flux is always theprojection of the disk’s surface onto the equatorial plane, thatis equal to the actual (inclined) disk surface multiplied by thecos( θ ) of the local inclination angle. This choice is made suchthat the description is consistent with the classical model thatregards the disk as infinitely thin (e.g. Hueso & Guillot 2005).The disk scale height h p (see Sect. 4.1.1) is used to define thesurface where the source function is evaluated. To avoid numeri-cal errors due to the small scale fluctuations of h p , we first fit thescale height to a powerlaw expression, H = η (cid:18) r (cid:19) p (28)such that it becomes a monotonically increasing smooth func-tion. The fitted parameters are shown in Fig. 4. The disk as-pect ratio is almost constant at 10 %. We choose to measure thesource function at h s = H and h s = H . The source function istherefore measured as S ( r ) = − ρ ( r , z ) (cid:2) u z ( r , z ) − u r ( r , z ) tan( θ ( r )) (cid:3) (29) = − ρ ( r , h s ( r )) (cid:34) u z ( r , h s ( r )) − u r ( r , h s ( r )) dh s dr ( r ) (cid:35) . Figure 24 shows snapshots of the measured S ( r ) at H (or-ange) and 3 H (green) at two resolutions. Closer to the equatorialplane, the source function is more centrally concentrated, as wehave discussed earlier due to the converging flow pattern (seeSect. 4.4.3 and 4.5.2). If we compare with the source functionsfrom Ulrich (1976) and Hueso & Guillot (2005), the measured S ( r ) is much steeper, meaning that the infall onto the disk ismore centrally concentrated for our simulated disk. Particularly,it means that most of the disk at r (cid:38) − ffi ciently far away from the mid-planeto avoid being a ff ected by the disk internal motions.Moreover, due to the finite size of sink particle accretion ra-dius, the flux within the disk scale height at small radius is notphysical. We therefore measure the flux at the horizontal sur-face z = dx (shown in gray). It is evident that the extreme high Article number, page 19 of 27 & A proofs: manuscript no. disk_formation_2020Feb r(AU)10 S ( k g / m / s ) H H dx r(AU)10 S ( k g / m / s ) H H dx r(AU)10 S ( k g / m / s ) H H dx r(AU)10 S ( k g / m / s ) H H dx R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 24.
The source function S ( r ) measured at H (orange) and 3 H (green) in R_ (cid:96)
14 and R_40ky_ (cid:96)
18 at 103 ky, R_ (cid:96)
14 and R_80ky_ (cid:96)
18 at138 ky ( from top to bottom ). The flux measured at the horizontal plane z = dx is shown in gray. The curves in lighter colors are the analyti-cal functions from pure hydrodynamic predictions which have the sameamount of total flux (see Fig. 3). The source function measured in thesimulations is much steeper than that of the classical model, with largeamount of infall happening at the central part of the disk. At a given ra-dius, the flux is generally greater at lower altitudes, indicating that thereis a strong radial flux in the interior of the disk. At small radii (below thesink accretion radius), the measured flux at H / H is no longer realisticand we used the measurement at z = dx to calculate the total flux. flux values at small radii are not physical. When getting veryclose to the disk center, the measurement of the vertical and hor-izontal flux becomes sensitive to the exact selection of the axisand there might be some confusion between the two. Therefore,when calculating the total mass flux, we use min(4 dx , H ) for thedisk surface. This introduces some uncertainties in the measure-ment, but it is the best we can do at this given resolution with theconstraints from the algorithm. r(AU)10 M ( M / y r ) M s ( H ) M s (3 H ) M r ( H ) M r (3 H ) M ( H ) M (3 H ) r(AU)10 M ( M / y r ) M s ( H ) M s (3 H ) M r ( H ) M r (3 H ) M ( H ) M (3 H ) r(AU)10 M ( M / y r ) M s ( H ) M s (3 H ) M r ( H ) M r (3 H ) M ( H ) M (3 H ) r(AU)10 M ( M / y r ) M s ( H ) M s (3 H ) M r ( H ) M r (3 H ) M ( H ) M (3 H ) R_ (cid:96)
14 (103 kyr)R_40ky_ (cid:96)
18 (103 kyr)R_ (cid:96)
14 (138 kyr)R_80ky_ (cid:96)
18 (138 kyr)
Fig. 25.
Mass accretion rate measured across the disk surface ( ˙ M s ( < r )),measured radially inside the disk ( ˙ M r ( r )), and inferred from the mea-sured α ( ˙ M α ( r )) for R_ (cid:96)
14 and R_40ky_ (cid:96)
18 at 103 ky, R_ (cid:96)
14 andR_80ky_ (cid:96)
18 at 138 ky ( from top to bottom ). Both positive and nega-tive values with | ˙ M | > − are plotted in logarithmic scale, and thisis connected with a linear range in between. Most of the surface fluxarrives where the slope of ˙ M s ( < r ) is the steepest. Below this radius,the disk in readily accreting in the radial direction only at very small ra-dius ( < α is not exactly the same as ˙ M r ( r ), although theamplitudes are comparable. The straightforward measure of the mass accretion rate, ˙ M r , isthe mass that crosses radius r inwards per unit time. This is cal-culated as˙ M r ( r ) = − π r h s (cid:88) − h s (cid:104) ρ u r ( r , z ) (cid:105) φ dz , (30) Article number, page 20 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core where ρ u r is the local radial mass flux. This radial accretion rateis shown in Fig. 25 with dashed curves for the integrated valueswithin H (orange) and 3 H (green). At small radii, the radial fluxwithin ± dx , which is larger than the scale height, is shown, andtherefore the green curves superimpose the orange curves. Asdiscussed earlier, the disk can be sometimes in partial expansion(negative ˙ M r ( r )), while in general the inner part of the disk (be-low ∼ H and 3 H curves shows that theradial accretion within the disk happens mostly in upper layers.This can be seen in the two lower panels of Fig 25. Beyond ∼ H ) curve is often above the orange ( H ) curve,meaning more accretion (positive) or less decretion (negative).In the region slightly beyond the sink accretion radius (secondvertical grey line), the value of ˙ M r ( r ) increases toward the centerof the disk in all four panels, meaning either the disk is beingemptied or there is a surface source.To examine this, we also calculated the surface mass accre-tion rate, that is, the source function integrated over the disk sur-face. We choose to perform this integration between 0 and r ,giving˙ M s ( r ) = r (cid:88) S ( r )2 π rdr . (31)This is shown on the same figure with solid curves for sourcefunctions measured at H (orange) and 3 H (green). Due to thenumerical limits described above, we use the flux at the surface z = dx whenever H is smaller than this value. Note that forR_ (cid:96)
14 at 103 kyr, the disk is so small that H < dx for almostthe entire disk.Most of the mass penetrates the disk surface where the slopeof ˙ M s ( r ) is the steepest. The total mass that falls on to thedisk, with significant contribution from this zone at a few AU,is broadly comparable with the mass that falls radially inwardsonto the star ( ˙ M r ( r ) at small radii). As a general remark, the in-falling material from the envelope can reach quite small radiusof a few AU before penetrating the disk surface. Inside the disk,the inward accretion is well established only for the innermost ofthe disk, while some of the mass that falls from the surface ontothe disk might move outwards. This expansion reaches a shockgenerated by the radial infall form the envelope near the equa-torial plane and brings the material to high altitudes, possiblygenerating a circulation in the outer part of the disk.Radial accretion behavior is present in all panels in themarginally-resolved region, that is, a few times the sink ac-cretion radius. Nonetheless, we would like to caution that thiszone is primarily dominated by the numerical viscosity, andthis is probably why the accretion is readily established. Thisshould not have strong consequence on the long-term stellarmass growth, since the infall from large-scales is the dominantfactor. However, such numerical e ff ect does a ff ect whether one isable to properly describe the evolution of the star-disk system onshort time-scales, which remains a computational challenge (e.g.Ku ff meier et al. 2018). In the outer part of the high-resolutionruns, both restarts show oscillating radial motions. This againsupports our proposition that the disk is highly sensitive to theenvironment, while the role of the disk as a mass bu ff er betweenthe envelope and the star is less evident.Figure 25 also shows ˙ M α , the mass accretion rate estimatedfrom the α parameter (see Sect. 4.4.3, Eq. (14)). The values of˙ M r and ˙ M α do not match perfectly with each other. However,they are similar in magnitude.
5. Comparison to standard alpha-disk models
The present paper shows that the disk formed during the infallof the dense core has a number of similarities to, but also dif-ferences with classical isolated disks assumed in standard alpha-disk models. Among the most notable similarities we report:1. The Keplerian disk structure is established rapidly in lessthan 61 kyr, a timescale shorter than the envelope infalltimescale.2. The disk mass is about 1% of the stellar mass after about 80kyr (Fig. 7). This is a typical disk / star mass ratio (like forthe MMSN, see e.g. Hayashi 1981). However, note that thestellar mass is only 0 . M (cid:12) at the end of the simulation sothat we cannot confirm for the moment that this results isvalid up to one Solar Mass. It is therefore premature to statethat there is enough mass to build the Solar System. More-over, we caution again that the disk mass is (unfortunately)sensitive to the disk inner boundary condition ( n acc ).3. The mid-plane of the disk is almost devoid of MHD turbu-lence with α r < − . This is in agreement with previousstudies arguing that a protoplanetary disk’s mid-plane shouldbe dead because of lack of ionization (see Turner et al. 2014,for a review). Several authors showed that low-values of α inthe disk mid-plane can lead to e ffi cient planetesimal forma-tion (Dr ˛a˙zkowska & Dullemond 2014, 2018; Charnoz et al.2019). The resulting disk could potentially be the birthplaceof a massive planetesimal population. However, we cautionthat the Toomre Q parameter is high in our simulation. Plan-etesimal formation can be allowed if Σ is higher in reality(see Sect. 4.4.4), or if the metallicity is enhanced towards thedisk mid-plane with dust enrichment of as high as a factor (cid:39)
2, as indeed suggested by the recent work of Lebreuillyet al. (2019, 2020).Among the notable di ff erences with classical alpha-diskmodels we report:1. The disk in our simulation is truncated at about 20 AUdue to magnetic braking, consistent with Hennebelle et al.(2016). This may imply that the classical picture of a vis-cously spreading disk, up to several tens, or one hundred AU,is maybe not the general rule, whereas such disk does exist(see e.g. Isella et al. 2010). However, recent ALMA obser-vations showed that the majority of class-0 disk are likelysmall (Andrews et al. 2018; Maury et al. 2019). The absenceof viscous spreading may preclude the outward transport ofCAIs due to disk outward viscous relaxation (Yang & Ciesla2012; Pignatale et al. 2018). However, alternative solutionscould exist as since transient outward directed flow are de-tected in the disk (see Sect. 4.5) and could potentially driveCAIs outward from the proto-star. Alternatively, the disk canpossibly expand in a later phase when the envelope is almostexhausted.Interestingly, the truncation of the disk at 20 AU is in agree-ment with some versions of the Nice Model, as well as themore generalized Giant Planets instability model which ap-peared later. In the recent view of Solar System formation,the giant planets formed in a much more compact config-uration than today, within about 15-20 AU of the Sun (seee.g. Gomes et al. 2005; Tsiganis et al. 2005; Morbidelli &Nesvorný 2020). The main argument for this is the neces-sity to form a low mass Kuiper Belt with a correct orbitalarchitecture in which an accumulation of bodies is observed Article number, page 21 of 27 & A proofs: manuscript no. disk_formation_2020Feb in mean-motion resonance with Neptune. To reach the ob-served state, it is necessary to invoke an outward migra-tion of Uranus and Neptune starting from 15-20 AU. Whenthe gas disperses, the four giant planets are still in a com-pact configuration with Jupiter and Saturn close to their 3:2mean motion resonance. This scenario is in agreement withdisk / planets evolutionary models (Bitsch et al. 2015) at leastfor Jupiter and Saturn.The Nice model was originally aimed to explain the lunarbombardment (Gomes et al. 2005). The key idea is that theorbits of giant planets were unstable during the S. S. evolu-tion. If the giant planets were formed in a more compact con-figuration than today (as suggested by the Nice Model, or asa result of evolution within the protoplanetary disk, see Kley& Crida 2008; Walsh et al. 2011), then they should have, atsome moment, su ff ered from a global dynamical instability,which lead to a very short and violent reorganization of theirorbits (Deienno et al. 2017; Clement et al. 2018). During thisre-arrangement, small bodies belts were destabilized, result-ing in massive bombardment of nearby planets. For example,the outward migration of Neptune may have pushed outwardthe original Kuiper belt, that could have been initially con-tained within about 20 AU, rather than at the present posi-tion of about 50 AU. Therefore, it is possible that the youngS. S. was more compact than today, and reached its presentdimensions after such an instability (Levison et al. 2008).In the original version of the Nice Model, the giant-planetsinstability happened at 800 Myr afters CAIs. More recentmodels of Giant Planets instability, using configuration ofthe four giant planets at the end of the solar nebula obtainedfrom hydrodynamical simulations of planet migration, alsoconverge to a formation of Neptune within 25 AU (Nesvorný& Morbidelli 2012). de Sousa Ribeiro et al. (2020) suggestedthat this event could have happened as early as in the firstmillions of years. Forming Uranus and Neptune at 30 AUand beyond would have disrupted the Kuiper belt and in par-ticular the population of dynamically cold objects (low or-bital eccentricity and inclination). This suggests that the diskitself was initially in a compact configuration, and was notextending presumably to 100 AU.Nonetheless, the disk in our simulation also experienced arapid grow in size after ∼
120 kyr. This suggest that the disksize is sensitive to the properties of the infalling material,which can be highly variable in time (see also Vorobyov et al.2015; Ku ff meier et al. 2017). Therefore, more constraint isneeded if one wishes to reconstruct the history of the SolarSystem.2. The formed disk is significantly hotter than classical isolateddisks. The snow-line starts at about 10 AU, and shifts in-wards down to 5 AU, where it remains constant after (untilthe end of the simulation) while the accretion drops to about10 − M (cid:12) / yr . This temperature profile is notably higher thanin passive radiatively heated disks (see e.g. Chiang & Gol-dreich 1997), where the snow-line falls around 2 AU (in-dependently of the accretion rate, and assuming a 2 . R (cid:12) star with T = − M (cid:12) / yr, consistent with the de-tailed model of D’Alessio et al. (1998) that finds a snow-line around 2 AU for an accretion rate of ∼ − M (cid:12) / yr. Asnow-line located at 5 AU would be the ideal position totrigger fast planetesimal formation of the Jupiter core, as ithas been demonstrated that streaming-instability and plan- etesimals formation is strongly enhanced at the snow-line(Dr ˛a˙zkowska & Alibert 2017; Charnoz et al. 2019).Nonetheless, the temperature of the disk should be inter-preted with caution since it is subject to some model un-certainties. First, the prescription by Kuiper & Yorke (2013)was used to describe the proto-stellar irradiation on theHosowaka track. This heating could possibly be excessiveand overheats the surrounding of the proto-star. It can be seenfrom Fig. 11 that the snow-line is shifted to 2-3 AU in the ab-sence of this heating term. Secondly, the the radiative trans-fer was treated with a flux-limited di ff usion scheme, whichis mostly valid for infrared radiation, while optical photonsmay escape more easily. In reality, the absorption of opticalradiation by dust and the re-emission in infrared wavelengthsshould not be complete and should be further investigatedwith care.3. The disk surface density does not follow the classical -1.5 ex-ponent of the Minimum Mass Solar Nebula (Hayashi 1981)but is closer to -0.5, with a sharp drop beyond R kep . Thisis remindful of observed protoplanetary disks, see e.g. Isellaet al. (2010) who suggested that the surface densities have anexponent > -1 and that the outer edge of the disks is between23 and 27 AU.
6. Comparison with similar numerical works
When looking back into the previous numerical studies, someessential conclusions could be drawn along with our currentwork. Primarily, the disk sensitively responds to its environment,that is, the collapsing prestellar core. Therefore, class-0 / I disksshould not be studied in the same way as isolated disks. Sec-ondly, whether the disk has significant impacts on the evolutionof the star remains debatable. Bursty mass accretion may resultfrom fragmentation within the disk, while the numerical originof such gravitational instabilities is not yet completely ruled out.An alternative way of regarding the disk as an accessory of theprotostar and an inevitable existence due to residual angular mo-mentum is emerging. Evolution of the star is determined mostlyfrom the core scale, and partly mediated by the disk. This devi-ates from the widely accepted view that mass first falls onto thedisk and then the star accretes from the disk.Vorobyov & Basu (2010, 2015); Vorobyov et al. (2015) em-ployed 2D thin-disk approximation. Such models cannot de-scribe the mass accretion that reaches the star through upperlayers and everything seems to transit through the disk. In theembedded phase, the stability conditions in the disk is mostlyreflected in the fluctuations (bursts) of the mass accretion rate,while the accretion rate from the envelope to the disk and thatfrom the disk to the star are broadly consistent (for class-0 / I,where envelope infall is significant).Ku ff meier et al. (2018) re-simulated small regions around sixsink particles selected from a giant molecular cloud simulation atincreased resolution of 2 AU for 200 kyr. One of those runs wasfurther evolved for 1000 yr at 0.06 AU resolution. These param-eters are similar to ours and we share some common findings.First, the overall mass accretion history is likely determined bylarge-scale flow properties (see also Ku ff meier et al. 2017), in-stead of the viscous accretion process in the disk. Second, thesink mass accretion mediated through the disk is bursty, thoughit is possibly of numerical origin. Finally, significant mass infallarrives onto the star without transiting radially through the disk. Article number, page 22 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core
7. Conclusions
We have simulated the collapse of a magnetized prestellar core,considering the non-ideal MHD e ff ects of the ambipolar di ff u-sion that rids the gas of its magnetic field. The temperature wascalculated with a flux-limited di ff usion scheme for the radiativetransfer. We followed the evolution during the first 100 kyr af-ter the protoplanetary disk formation, and we performed high-resolution restarts in parallel to study the numerical convergenceas well as the detailed interior structure of the disk. This work isamongst the first works that try to self-consistently form a proto-planetary disk and resolve the interior structure at the same time.Nonetheless, it should be kept in mind that the numerical con-vergence is not strictly demonstrated in this work and challengesremain.The collapse problem is known to be computationally chal-lenging due to the large dynamical ranges. For such reason, theformation of the protoplanetary disk has rarely been followedfor long duration. On the other hand, the detailed structure ofthe disk interior is usually studied with a prescribed disk, ora shearing box that contains a small region of the disk. In thisstudy, we have followed the evolution of the disk during a rela-tively long period of time at 1 AU resolution. The disk is formedself-consistently from a prestellar core that is a few thousands ofAU in radius. Such compromise of resolution allows us to reachmore evolved stages within reasonable amount of computationtime. In order to study the detailed interior dynamics of the disk,we selected snapshots from the canonical simulation and rerunwith increased resolution down to 0.06 AU.From this work, we were able to follow the time evolutionof disk global properties such as its mass and size. The disk isestablished roughly 60 kyr after the beginning of the simulationand grows slowly in size. The mass accretion rate is a decreas-ing function of time and the disk gradually approaches an iso-lated state. The global appearance (density structure) of disk isquasi-stationary, and has mass of 1 . α − disk model that describes the turbulenttransport, linked to the local surface density and temperature, isunlikely to be applicable for a disk in the embedded phase.Due to the small size of the disk ( <
25 AU) and to the pres-ence of strong MHD turbulence in the outer regions, all planetsmust form within ∼
30 AU ( R kep ). This is qualitatively consis-tent with the modern view of Solar System formation, where allplanets are suggested to have formed well inside 20 AU (e.g.Raymond & Izidoro 2017; Morbidelli et al. 2018) and most plan-etesimals may have accreted up to 30 AU (Gomes et al. 2004;Nesvorny et al. 2013; Nesvorný et al. 2018). The temperatureprofile evolution can also be extracted and allows the determi-nation of the migration of the snow-line. The radial temperatureprofile is consistent with a powerlaw exponent between -1 and-0.5. This is somewhat between a purely irradiated disk and aviscously heated disk. However, the physical origin of this pro-file is likely the radiative di ff usion in the disk radial directioninstead of the aforementioned mechanisms. The snow-line startsat about 10 AU, and ends at about 5 AU after about 100 kyrevolution. This is much further away than commonly assumedfor isolated and passively irradiated disks. This result should beregarded with care because it might have some numerical ori-gins. First of all, the FLD method treats only infrared radiation and thus could trap more energy than it should in reality. Sec-ondly, the temperature is highly sensitive to the density (opti-cal depth) of the disk, which in turn depends on the disk innerboundary conditions. Moreover, it should also be kept in mindthat the stellar mass has only reached ∼ . M (cid:12) in this study.Initial conditions with higher mass should be further explored.The disk has an inner part with a shallow powerlaw surfacedensity profile ( (cid:38) − (cid:38) H , reaching the interior of the disk withoutpenetrating the disk surface at large radius , which is also sug-gested by Ku ff meier et al. (2018). The mid-plane of the disk,on the other hand, can be sometimes expanding. The expandingmaterial from the disk mid-plane sometimes exceed the Keple-rian radius, R kep , of the disk and confronts the infalling materialin the magnetized zone between R kep and R mag . Some material ispushed up vertically and fall down again through the disk surfaceat smaller radius. This can potentially generate some circulationof material, which remains to be verified using tracer particles orpassive scalar fields. In particular, the source function evaluatedfrom our simulations is much more centrally peaked comparedto classical models, and most of the mass that ends up inside thestar does not transit through the whole disk. Acknowledgements.
We thank the referee, Enrique Vazquez-Semadeni, forthoutough comments that have significantly improved the manuscript. Thiswork was granted access to HPC resources of CINES under the allocationx2014047023 made by GENCI (Grand Equipement National de Calcul Inten-sif). Y.-N. Lee acknowledges the financial support of the UnivEarthS Labexprogram at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), MOST (108-2636-M-003-001, 109-2636-M-003-001), and the MOEYushan Young Scholar program.
References
Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109ALMA Partnership, Vlahakis, C., Hunter, T. R., et al. 2015, ApJ, 808, L4Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41Bai, X.-N. & Stone, J. M. 2010, ApJ, 722, 1437Bai, X.-N. & Stone, J. M. 2013, ApJ, 769, 76Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214Balbus, S. A. & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650Béthune, W. & Latter, H. 2020, MNRAS, 494, 6103Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A,549, A124Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575,A28Bleuler, A. & Teyssier, R. 2014, MNRAS, 445, 4015Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50Chiang, E. I. & Goldreich, P. 1997, ApJ, 490, 368Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icarus, 311,340Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011,A&A, 529, A35Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ,600, 279D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411Dapp, W. B. & Basu, S. 2010, A&A, 521, L56de Sousa Ribeiro, R., Morbidelli, A., Raymond, S. N., et al. 2020, Icarus, 339,113605Deienno, R., Morbidelli, A., Gomes, R. S., & Nesvorný, D. 2017, AJ, 153, 153Desch, S. J., Kalyaan, A., & O’D. Alexander, C. M. 2018, ApJS, 238, 11Dr˛a˙zkowska, J. & Alibert, Y. 2017, A&A, 608, A92
Article number, page 23 of 27 & A proofs: manuscript no. disk_formation_2020Feb
Dr˛a˙zkowska, J. & Dullemond, C. P. 2014, A&A, 572, A78Dr˛a˙zkowska, J. & Dullemond, C. P. 2018, A&A, 614, A62Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71Fromang, S. & Nelson, R. P. 2006, A&A, 457, 343Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466Gomes, R. S., Morbidelli, A., & Levison, H. F. 2004, Icarus, 170, 492Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35Hayashi, C., Nakazawa, K., & Mizuno, H. 1979, Earth and Planetary ScienceLetters, 43, 22Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830,L8Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020a, A&A, 635,A67Hennebelle, P., Commercon, B., Lee, Y.-N., & Chabrier, G. 2020b, arXiv e-prints, arXiv:2010.03539Henriksen, R., Andre, P., & Bontemps, S. 1997, A&A, 323, 549Hirano, N., Ohashi, N., Dobashi, K., Shinnaga, H., & Hayashi, M. 2002, in 8thAsian-Pacific Regional Meeting, Volume II, ed. S. Ikeuchi, J. Hearnshaw, &T. Hanawa, 141–142Hueso, R. & Guillot, T. 2005, A&A, 442, 703Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90Inutsuka, S.-i. 2012, Progress of Theoretical and Experimental Physics, 2012,01A307Isella, A., Carpenter, J. M., & Sargent, A. I. 2010, ApJ, 714, 1746Kenyon, S. J. & Hartmann, L. 1987, ApJ, 323, 714Kley, W. & Crida, A. 2008, A&A, 487, L9Kruijer, T. S. & Kleine, T. 2017, Earth and Planetary Science Letters, 475, 15Ku ff meier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475,2642Ku ff meier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7Kuiper, R. & Yorke, H. W. 2013, ApJ, 772, 61Lam, K. H., Li, Z.-Y., Chen, C.-Y., Tomida, K., & Zhao, B. 2019, MNRAS, 489,5326Larson, R. B. 1969, MNRAS, 145, 271Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96Lebreuilly, U., Commerçon, B., & Laibe, G. 2020, A&A, 641, A112Levison, H. F., Morbidelli, A., Van Laerhoven, C., Gomes, R., & Tsiganis, K.2008, Icarus, 196, 258Li, Z. Y., Banerjee, R., Pudritz, R. E., et al. 2014, in Protostars and Planets VI,ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 173Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2014, MNRAS, 438, 2278Machida, M. N., Matsumoto, T., & Inutsuka, S.-i. 2016, MNRAS, 463, 4246Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3Marchand, P., Commerçon, B., & Chabrier, G. 2018, A&A, 619, A37Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016,A&A, 587, A32Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G.2012, ApJS, 201, 24Masunaga, H., Miyama, S. M., & Inutsuka, S.-i. 1998, ApJ, 495, 346Matsumoto, T. & Hanawa, T. 2003, ApJ, 595, 913Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477, 2760Morbidelli, A., Bitsch, B., Crida, A., et al. 2016, Icarus, 267, 368Morbidelli, A. & Nesvorný, D. 2020, Kuiper belt: formation and evolution, ed.D. Prialnik, M. A. Barucci, & L. Young, 25–59Morbidelli, A., Nevorny, D., Laurenz, V., et al. 2018, LPI Contributions, 2107,2005Nakamoto, T. & Nakagawa, Y. 1994, ApJ, 421, 640Nesvorný, D. & Morbidelli, A. 2012, AJ, 144, 117Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, NatureAstronomy, 2, 878Nesvorny, D., Vokrouhlicky, D., & Morbidelli, A. 2013, in AAS / Division forPlanetary Sciences Meeting Abstracts / Division for Planetary Sci-ences Meeting Abstracts, 508.03Ohashi, N., Lee, S. W., Wilner, D. J., & Hayashi, M. 1999, ApJ, 518, L41Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E. 2018, ApJ, 867, L23Raymond, S. N. & Izidoro, A. 2017, Icarus, 297, 134Riols, A., Ogilvie, G. I., Latter, H., & Ross, J. P. 2016, MNRAS, 463, 3096Sakai, N., Oya, Y., Higuchi, A. E., et al. 2017, MNRAS, 467, L76Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A,410, 611Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33Shu, F. H. 1977, ApJ, 214, 488Shu, F. H., Li, Z.-Y., & Allen, A. 2004, ApJ, 601, 930Suzuki, T. K., Ogihara, M., Morbidelli, A. r., Crida, A., & Guillot, T. 2016, A&A,596, A74 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed.H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339Teyssier, R. 2002, A&A, 385, 337Tomida, K., Machida, M. N., Hosokawa, T., Sakurai, Y., & Lin, C. H. 2017, ApJ,835, L11Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015,ApJ, 810, L26Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI,ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 411Turner, N. J. & Sano, T. 2008, ApJ, 679, L131Ulrich, R. K. 1976, ApJ, 210, 377Van Kooten, E. M. M. E., Wielandt, D., Schiller, M., et al. 2016, Proceedings ofthe National Academy of Science, 113, 2011Vaytet, N. & Haugbølle, T. 2017, A&A, 598, A116Verliat, A., Hennebelle, P., Maury, A. J., & Gaudel, M. 2020, A&A, 635, A130Vorobyov, E. I. & Basu, S. 2010, ApJ, 719, 1896Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115Vorobyov, E. I., Lin, D. N. C., & Guedel, M. 2015, A&A, 573, A5Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627,A154Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M.2011, Nature, 475, 206Ward-Thompson, D., André, P., Crutcher, R., et al. 2007, in Protostars and Plan-ets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 33Weidenschilling, S. J. 1977, Ap&SS, 51, 153Whitworth, A. & Summers, D. 1985, MNRAS, 214, 1Whitworth, A. P. & Ward-Thompson, D. 2001, ApJ, 547, 317Wurster, J. & Li, Z.-Y. 2018, Frontiers in Astronomy and Space Sciences, 5, 39Yang, L. & Ciesla, F. J. 2012, Meteoritics and Planetary Science, 47, 99Youdin, A. N. 2011, ApJ, 731, 99Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868Zhao, B., Caselli, P., Li, Z.-Y., et al. 2020, MNRAS, 492, 3375
Article number, page 24 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core
Appendix A: Identification of the disk
Because we are interested in the interior properties of the diskas well as its accretion onto the central star, we choose the sinkparticle as the center of the star-disk system. Despite some non-axisymmetric oscillation of the disk, the gravitational e ff ectsfrom the central star cancels out when calculating azimuthallyaveraged quantities, α in particular. This allows a discussionof the disk properties in coherence with the axisymmetric diskmodels (which is the case for almost all existing models). Thedisk is identified in two steps. The first step is based on puredensity thresholding, and the second step is more refined andtakes into account the radial and vertical mass distributions ofthe disk (see Sect. 4.1).The disk is significantly marked by a density jump, which isused as a simple first identification. The rotating dense disk issurrounded by a collapsing di ff use envelope. The procedure is asfollows:1. Calculate the inverse cumulative density function (CDF) ofmass against density, that is, the amount of mass that is atdensity above a certain value, M ( ρ ) = (cid:90) ∞ ρ dm . (A.1)This results in a decreasing function of ρ that has a plateau,which corresponds to the density jump between the disk andthe envelope.2. Define the disk threshold ρ th density as the point where thederivative of the CDF is highest, i.e. where the CDF is flat-test. This writes: ρ th = arg max d log M ( ρ ) / d log ρ. (A.2)3. Select all the cells that have density ρ ≥ ρ th .4. Calculate the angular momentum of these cells with respectto the sink particle and define this as the disk axis.A simple density thresholding is robust to the turbulent fluc-tuations for, while not always, most of the time. As the high den-sity regions dominate the mass inside the disk, the rotation axishere-defined is robust to slightly di ff erent definitions of the disk.The rotation in the initial condition is aligned to the − z axis, andthe axis of the disk is very close to this direction. Appendix B: Disk vertical temperature distribution
Figure B.1 shows and example of disk vertical temperature dis-tribution, where the temperature from each cell is plotted di-rectly. Within a few times the scale height, the disk is basicallyisothermal, and decrease in temperature is only significant faraway from the disk mid-plane. Vertical temperature variation ismore clearly seen at small radii, i.e., within 7 AU, where thevertical structure of the disk is, however, not well resolved. Thetemperature dispersion amongst cells at same radial distance isactually sometimes more prominent than the vertical variation.
Appendix C: The stress tensors linked to thetransport
Here we define self-coherently the stress tensors used in thiswork such that there should be no confusion. Following Balbus z (AU) T ( K ) Time = 89 kyr, M* = 0.34 Ms r=1.0 AUr=3.0 AUr=5.0 AUr=7.0 AUr=9.0 AUr=11.0 AUr=13.0 AUr=15.0 AUr=17.0 AUr=19.0 AUr=21.0 AU
Fig. B.1.
Temperature plotted against altitude, z , at various distance tothe central star, r , from R_ (cid:96)
14 (recall that the disk aspect ratio is ∼ . ff erence between the grid direction and the selected disk axis. & Hawley (1998), the angular momentum conservation averagedin azimuthal direction becomes ∂∂ t (cid:104) ρ r Ω (cid:105) φ = ∇ · r [ (cid:104) ρ r Ω u p (cid:105) φ + T ] = , (C.1)where T is the poloidal stress tensor T = (cid:104) ρ u φ u p − B φ B p / π (cid:105) φ , (C.2)and W = T / (cid:104) ρ (cid:105) φ . Instead of integrating over the whole verticaldomain, we integrated over the disk scale height [ − H , H ], anddefine a mass weighted mean (cid:104) X (cid:105) ρ = π Σ H ∆ r (cid:90) H − H (cid:90) r +∆ r / r − ∆ r / (cid:90) π ρ Xd φ drdz , (C.3)where ∆ r is a radial averaging scale chosen such that the pro-file is not significantly a ff ected by fluctuations. The choice of ∆ r = dx is reasonable for most variables, while for the stresstensors larger values are used and this is discussed in the Ap-pendix D. It is easily shown from this definition that the disksurface density within a the thickness ± H , Σ H , can be found bycalculating (cid:104) (cid:105) ρ = ∂ Σ H ∂ t + r ∂ ( r Σ H (cid:104) u r (cid:105) ρ ) ∂ r + (cid:104) ρ u z (cid:105) H φ − (cid:104) ρ u z (cid:105) − H φ = , (C.4)where the superscript ± H denotes the values at ± H . The secondterm is the radial mass flux, and the last two terms are the fluxacross the disk surfaces.Likewise, Eq. (C.1) becomes: ∂∂ t (cid:16) Σ H r (cid:104) Ω (cid:105) ρ (cid:17) + r ∂∂ r (cid:16) r Σ H (cid:104) Ω u r (cid:105) ρ + r Σ H (cid:104) W r φ (cid:105) H (cid:17) (C.5) + r (cid:16) Ω H (cid:104) ρ u z (cid:105) H φ − Ω − H (cid:104) ρ u z (cid:105) − H φ (cid:17) + r (cid:16) (cid:104) ρ W z φ (cid:105) H φ − (cid:104) ρ W z φ (cid:105) − H φ (cid:17) = , where we also define a vertical density weighted average as (cid:104) X (cid:105) H = Σ H ∆ r (cid:90) H − H (cid:90) r +∆ r / r − ∆ r / (cid:104) ρ (cid:105) φ Xdrdz . (C.6)The second term on the first line contains the the radial transportdue to the laminar flow and the stress tensors, while the secondline corresponds to the same contributions across the disk sur-face. Article number, page 25 of 27 & A proofs: manuscript no. disk_formation_2020Feb
Assuming a stationary disk, we ignore the time derivatives.Combining eqs. (C.4) and (C.5) and assuming that the disk issymmetric with respect to the mid-plane ( Ω H = Ω − H , rememberthat Ω = (cid:104) u φ (cid:105) φ / r by definition), we can derive Σ H (cid:104) u r (cid:105) ρ ∂∂ r (cid:16) r Ω H (cid:17) + r ∂∂ r (cid:16) r Σ H (cid:104) W r φ (cid:105) H (cid:17) (C.7) + r (cid:16) (cid:104) ρ W z φ (cid:105) H φ − (cid:104) ρ W z φ (cid:105) − H φ (cid:17) − r ∂∂ r r Σ H (cid:42) ∂ Ω ∂ z (cid:90) ρ u r dz (cid:43) ρ = . The last two terms disappear when the integration is performedover the whole vertical extent with vanishing boundary condi-tions and the rotation profile is assumed vertically invariant.The radial mass accretion rate within the disk˙ M ( r ) = − π R Σ H (cid:104) u r (cid:105) ρ (C.8)no longer has a simple expression. In this work, we limit ourdiscussion to the vertical extent where the disk is practically Ke-plerian and thus the last term in Eq. (C.7) can be neglected. Wecan then infer the Eq. (14).There are two contributions to the mass accretion due to an-gular momentum removal. The first one is in the radial directionand the contribution is integrated across the disk thickness, forwhich we can define a mean α value such that α r = (cid:104) W r φ (cid:105) H (cid:104) c (cid:105) ρ = ˆ T Hr φ ρ H (cid:104) c (cid:105) ρ , (C.9)where ρ H = Σ H / (2 H ) is the disk vertical mean density. The sec-ond contribution comes from the angular momentum lost due tovertical transport across both surfaces of the disk, and we definethe corresponding α to be α z = ( T Hz φ − T − Hz φ ) ρ H (cid:104) c (cid:105) ρ = ( (cid:104) ρ (cid:105) H φ W Hz φ − (cid:104) ρ (cid:105) − H φ W − Hz φ ) ρ H (cid:104) c (cid:105) ρ . (C.10)The normalization is always done with respect to the bulk aver-aged pressure. Appendix D: Numerical discretizing effects on theevaluation of the stress tensor
The azimuthally-averaged stress tensor of the velocity is definedlocally in a narrow range of radius ∆ r around position r , suchthat rapid radial fluctuations are smoothed out: W p φ, loc = (cid:15) (cid:90) r + (cid:15)/ r − (cid:15)/ (cid:104) δ u p δ u φ (cid:105) φ dr , (D.1)where (cid:15) is a minimum value such that the function is smooth.When numerically evaluating the stress tensor, this average isperformed within a practically larger region. We discuss herethe e ff ects of discretized calculation on the stress tensor. Whencalculating the stress tensor in a larger finite region, there is anextra velocity di ff erence with respect to the mean velocity in thewhole region. We have thus W p φ, fin = ∆ r (cid:90) r +∆ r / r − ∆ r / (cid:104) ( ∆ u p + δ u p )( ∆ u φ + δ u φ ) (cid:105) φ dr (D.2) = W p φ, loc + ∆ r (cid:90) r +∆ r / r − ∆ r / (cid:104) ∆ u p ∆ u φ (cid:105) φ dr = W p φ, loc + W p φ, cor , where ∆ specifies the di ff erence between the mean value at r andthe mean value averaged in [ − ∆ r / , ∆ r / · stands for theradially averaged value. Below we performed a calculation withsome simplifying assumptions to estimate the e ff ect of averagingin a finite region on the radial component of the tensor.First of all, we assume Keplerian rotation such that the az-imuthal average gives (cid:104) u φ (cid:105) φ ( r ) = r Ω ( r ) = (cid:114) GM ∗ r . (D.3)In the radial direction, the accretion velocity is given by (cid:104) u r (cid:105) φ ( r ) = − ( r Σ W r φ, loc ) (cid:48) r Σ ( r Ω ) (cid:48) ∼ − η W r φ, loc r Ω = − η W r φ, loc (cid:114) rGM ∗ , (D.4)where η is a factor of order unity from the unknown profiles of Σ and W r φ, loc . By assuming the average inside ∆ r is simply thevalue at r , we can calculate ∆ u r ∆ u φ = − η W r φ, loc (cid:16) (cid:112) r + ∆ r / − √ r (cid:17) (cid:32) √ r + ∆ r / − √ r (cid:33) (D.5) = − η W r φ, loc (cid:32) − r + ∆ r / √ ( r + ∆ r / r (cid:33) To calculate the average in the whole region, we assume that W r φ, loc ( r ) and Σ ( r ) are slowly varying and can be locally approx-imated as a constant. We also neglect the disk geometry and per-form simply a linear integration along r . The correction term foraveraged stress tensor calculated in a non-negligible finite regionof [ − ∆ r / , ∆ r /
2] is thus W r φ, cor = ∆ r (cid:90) r +∆ r / r − ∆ r / (cid:104) ∆ u r ∆ u φ (cid:105) φ dr (D.6) = η ∆ r ∆ r / (cid:90) − ∆ r / (cid:104) ∆ (cid:48) u r ∆ (cid:48) u φ (cid:105) d ∆ (cid:48) r = η η W r φ, loc (cid:32) ∆ rr (cid:33) − (cid:32) ∆ rr (cid:33) + O (cid:32) ∆ rr (cid:33) , where η is also a factor of order unity from the simplifications.The stress tensor that we practically measure from the simula-tions is therefore W r φ, fin = W r φ, loc + W r φ, cor ≈ W r φ, loc + η η (cid:32) ∆ rr (cid:33) . (D.7)As long as ∆ r / r is small and W r φ, loc slowly varying, the correc-tion term remains negligible and does not a ff ect the estimationof the stress tensor. Nonetheless, the choice of ∆ r / r still has tobe reasonable such that the simplifying assumptions stay valid.Indeed, when measuring the α values from our simulations,the size of the binning has to be carefully chosen. As shownin Fig. D.1 (top panel), the value of the Reynolds α fluctuatesstrongly when the bin size in r is taken to be the smallest cellsize. Nonetheless, by increasing the bin size, the measured α converges towards a smooth function. This function does notseem to be dependent of the bin size as long as the bin size ≥ dx , confirming the above discussions.On the other hand, this does not seem to be true for the grav-itational α . The measurement remains noisy irrespectively of theincreased bin size. This might be a result of noises introduced bythe cartesian discretization for the disk system. Article number, page 26 of 27ueh-Ning Lee et al.: Protoplanetary disk formation from collapse of prestellar core r ( AU )10 R e y dx2dx4dx8dx16dx r ( AU )10 M a x dx2dx4dx8dx16dx r ( AU )10 G r a v dx2dx4dx8dx16dx Fig. D.1.
Top:
Reynods α calculated at various bin sizes in r . The func-tion is strongly fluctuating at small bin sizes, and becomes smooth forbins sizes ≥ dx . Bottom:
Gravitational α calculated at various bin sizesin r . The function fluctuates strongly between positive and negative val-ues regardless of the increased bin size. It is impossible to measure thegravitational α in the simulations in this studies, probably due to carte-sian discretization. The resolution dx = ..