Improving the thin-disk models of circumstellar disk evolution. The 2+1-dimensional model
aa r X i v : . [ a s t r o - ph . GA ] J un Astronomy & Astrophysics manuscript no. disk2+1D˙5˙YP c (cid:13)
ESO 2017June 2, 2017
Improving the thin-disk models of circumstellar disk evolution.The 2+1-dimensional model
Eduard I. Vorobyov , , and Yaroslav N. Pavlyuchenkov Institute of Fluid Mechanics and Heat Transfer, TU Wien, 1060, Vienna, Austria; e-mail: [email protected] Research Institute of Physics, Southern Federal University, Stachki Ave. 194, Rostov-on-Don, 344090 Russia; University of Vienna, Department of Astrophysics, Vienna, 1180, Austria Institute of Astronomy of the Russian Academy of Sciences, Pyatnitskaya str. 48, Moscow, 119017 RussiaPreprint online version: June 2, 2017
ABSTRACT
Context.
Circumstellar disks of gas and dust are naturally formed from contracting pre-stellar molecular cores duringthe star formation process. To study various dynamical and chemical processes that take place in circumstellar disksprior to their dissipation and transition to debris disks, the appropriate numerical models capable of studying thelong-term disk chemodynamical evolution are required.
Aims.
We aim at improving the frequently used two-dimensional hydrodynamical model for disk evolution in the thin-disk limit by employing a better calculation of the disk thermal balance and adding a reconstruction of the disk verticalstructure. Together with the hydrodynamical processes, the thermal evolution is of great importance since it influencesthe strength of gravitational instability and the chemical evolution of the disk.
Methods.
We present a new 2+1-dimensional numerical hydrodynamics model of circumstellar disk evolution, in whichthe thin-disk model is complemented with the procedure for calculating the vertical distributions of gas volume densityand temperature in the disk. The reconstruction of the disk vertical structure is performed at every time step via thesolution of the time-dependent radiative transfer equations coupled to the equation of the vertical hydrostatic equilib-rium.
Results.
We perform a detailed comparison between circumstellar disks produced with our previous 2D model and withthe improved 2+1D approach. The structure and evolution of resulting disks, including the differences in temperatures,densities, disk masses and protostellar accretion rates, are discussed in detail.
Conclusions.
The new 2+1D model yields systematically colder disks, while the in-falling parental clouds are warmer.Both effects act to increase the strength of disk gravitational instability and, as a result, the number of gravitationallybound fragments that form in the disk via gravitational fragmentation as compared to the purely 2D thin-disk simu-lations with a simplified thermal balance calculation. The presented method has a low time overhead as compared tothe purely 2D models and is particularly suited for the long-term simulations of circumstellar disks including compactchemical reaction networks.
Key words. protoplanetary disks – stars: formation – stars: protostars – hydrodynamics
1. Introduction
Circumstellar disks hold the key to our understanding ofstellar mass accumulation and planet formation. They formthanks to the conservation of angular momentum of a col-lapsing cloud core when the in-spiralling material hits thecentrifugal barrier near the stellar surface before landingonto the star. As the core collapse continues, new lay-ers of infalling material are added to the newborn diskat progressively larger radial distances, causing the diskto grow in both mass and size. Numerical simulationsand observations indicate that this process starts in theClass 0 phase of stellar evolution and continue until thethe parental core depletes or dissipates (Machida et al.,2010; Tobin et al., 2015). In this so-called protostellar orembedded phase, disks are usually most massive owingto continuing mass loading from the parental core andare often prone to the development of gravitational in-stability (Vorobyov & Basu, 2005; Tsukamoto et al., 2013;Lomax et al., 2014; Dong et al., 2016; Mayer et al., 2016). This is also the phase when dust growth and perhaps thefirst phases of planet assembly are likely to take place(ALMA Partnership, 2015).In the subsequent T Tauri phase of stellar evolution, thedisk slowly loses its mass owing to accretion onto the hoststar and expands in reaction to angular momentum redis-tribution within the disk. In this phase, dust growth fromsmall grains to planetesimals and planetary cores takesplace, finally leading to the emergence of protoplanets cov-ered with primordial atmospheres (e.g., Benz et al., 2014)and basic chemical ingredients are converted into complex(organic) species (Henning & Semenov, 2013). Finally, thecombined action of stellar accretion, planet formation, diskwinds and photoevaporation leads to the dispersal of thedisk gaseous component, revealing a debris disk consistingof solids which are left from the planet formation process.Many aforementioned phenomena during the entire diskevolution are controlled by the radiative input from thehost star and external environment; both largely deter-
Vorobyov & Pavlyuchenkov: The 2+1D model for disk evolution mine the temperature in the disk upper layers and set theminimum temperature in its midplane. In the steady-statedisk models, the radiative input can be taken into accountrather accurately using sophisticated ray-tracing or Monte-Carlo techniques (e.g. Woitke et al., 2009; Dullemond,2012; Akimkin et al., 2013). However, due to high compu-tational costs, these models lack a dynamical aspect, whichcan be important when considering the young disks proneto gravitational instability or more evolved disks whereinplanetesimals/planetary cores are subject to growth andmigration. To properly consider these effects, fully dynam-ical models are needed which are usually based the equa-tions of hydrodynamics complemented with a module thattakes the radiative input from the stars and/or externalenvironment into account.In the full 3D hydrodynamics simulations, the radiativeinput can be taken into account by solving the equationsof radiation transfer, usually with simplifying assumptionssuch as frequency integrated opacities and diffusion approx-imation (e.g., Klahr & Kley, 2005; Tsukamoto et al., 2015),because the more accurate frequency-dependent models(e.g., Kuiper et al., 2010) are computationally expensive.Even with these simplifying assumptions, 3D models arean inconvenient tool to make a parameter space study ofthe disk evolution for many model realizations and manyorbital periods. For this purpose, 2D thin-disk models ofdisk evolution have been routinely employed (e.g. Vorobyov,2011; Zhu et al., 2012; Reg´aly, 2013; Gyergyovits et al.,2014) thanks to their low computational costs. In thesemodels, the equation for thermal balance is usually com-plemented with cooling and heating functions which takethe radiative input and local disk cooling into account. Theform of these cooling and heating terms may differ slightlyfrom study to study, but they are essentially hinged on thecalculated midplane and surface temperatures of the diskand the optical depth from the disk surface to the mid-plane. The midplane temperature is usually set equal to thehydrodynamic temperature, the surface temperature is cal-culated assuming the black-body character of the incidentstellar and background radiation, and the optical depth iscalculated assuming the vertically isothermal temperaturedistribution.While having advantages in simplicity and low computa-tional costs, this approach has obvious weaknesses. First, itis not clear if the 2D hydrodynamic temperature is indeedrepresentative of the disk midplane temperature. Second,this method provides little information on the disk ver-tical structure, essentially assuming that the disk is ver-tically isothermal, which may not be the case. For in-stance, passively irradiated disks are known to exhibit pos-itive temperature gradients in the vertical direction (e.g.,Dullemond et al., 2002), while dense gaseous clumps andspiral arms in gravitationally unstable disks may be char-acterized by a more complicated vertical temperature dis-tribution (Vorobyov et al., 2014). To circumvent this dif-ficulty, various forms of the vertical density distributioncan be adopted (Gaussian, exponential), but this proce-dure cannot be easily applied to the vertical temperaturedistribution. As a result, the 2D models have limited ap-plicability to studying, e.g., chemical reactions which areknown to sensitively depend on the gas temperature (whichis likely to be varying in the vertical direction as well).This paper presents a method that addresses the afore-mentioned weaknesses by means of coupling the gas dy- namics computations in the thin-disk limit and calcula-tions of the disk vertical structure. The gas dynamics inthe disk plane is calculated using the hydrodynamics equa-tions, while the disk vertical structure is calculated usingthe equations of radiation transfer and hydrostatic balancein the vertical direction. As a result, we retrieve the fullthree-dimensional density and temperature structure of thedisk at every time step, which is lacking in purely two-dimensional models. The presented method has low timeoverhead as compared to the purely 2D thin-disk modelsand it is faster than fully 3D models since the disk gravita-tional potential is found using a convolution theorem (seeEquation 4), which is not applicable to the fully 3D mod-els formulated in the curvilinear coordinate systems. Theadopted 1D radiative transfer in the vertical direction isalso much faster than the fully 3D method. Our method istherefore well suited for the long-term simulations of cir-cumstellar disks. We compare the disk evolution calculatedusing this 2+1D approach with purely 2D simulations andbriefly discuss the applicability of our new method to calcu-lating the chemical evolution in circumstellar disks. The pa-per is organized as follows. In Section 2, the hydrodynamicsequations in the thin-disk limit are reviewed. In Section 3.1,we formulate the modifications made to the thin-disk modelto improve the thermal balance calculations in the thin disk.In Section 4, we compare the disk evolution in the 2+1Dand 2D approaches. The model caveats and future improve-ments are discussed in Section 5. The main conclusions aresummarized in Section 6. The Appendix presents details ofthe solution procedure used to calculate the disk verticalstructure.
2. Model equations in the thin-disk limit
The equations of mass, momentum, and energy transportdescribing the dynamics of circumstellar disks in the thin-disk limit can be formulated as follows: ∂ Σ ∂t = −∇ · (Σ v ) , (1) ∂∂t (Σ v ) + ∇ · (Σ v ⊗ v ) = −∇P + Σ g + ∇ · Π , (2) ∂e∂t + ∇ · ( e v ) = −P ( ∇ · v ) − Λ + Γ + ( ∇ v ) : Π , (3)where Σ is the mass surface density, e is the internal energyper surface area, P is the vertically integrated gas pressurecalculated via the ideal equation of state as P = ( γ − e with γ = 7 / v = v r ˆ r + v φ ˆ φ is the velocity in the disk plane,and ∇ = ˆ r ∂/∂r + ˆ φ r − ∂/∂φ is the gradient along the planarcoordinates of the disk. The gravitational acceleration inthe disk plane, g = g r ˆ r + g φ ˆ φ , takes into account self-gravity of the disk and the gravity of the central protostarwhen formed. Disk self-gravity is found by solving for thePoisson integralΦ( r, φ ) = − G Z r out r sc r ′ dr ′ × Z π Σ( r ′ , φ ′ ) dφ ′ q r ′ + r − rr ′ cos( φ ′ − φ ) , (4) orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 3 where r sc and r out are the radial positions of the compu-tational inner and outer boundaries. This integral is calcu-lated using a FFT technique which applies the 2D Fourierconvolution theorem for polar coordinates and allows forthe non-periodic boundary conditions in the r -direction byeffectively doubling the computation domain in this co-ordinate direction and filling it with zero densities (seeBinney & Tremaine, 1987, Sect. 2.8). Turbulent viscosityis taken into account via the viscous stress tensor Π , theexpression for which can be found in Vorobyov & Basu(2010). The kinematic viscosity needed to calculate theviscous stress tensor is found adopting the Shakura andSunyaev parameterization (Shakura & Sunyaev, 1973), sothat ν = αc s h , where c s = p γ P / Σ is the sound speed and h is the disk scale height. The α -parameter is linked to the in-ferred strength of the magneto-rotational instability (MRI)following the method that takes the MRI active/inactivestates into account as described in (Bae et al., 2014).We use the following form for the cooling termΛ in equation (3) based on the analytical solution ofthe radiation transfer equations in the vertical direction(Dong et al., 2016):Λ = 4 τ P σT τ P + τ R τ P , (5)where τ R = κ R Σ / and τ P = κ P Σ / the Rosseland andPlanck optical depths to the disk midplane, κ P and κ R the Planck and Rosseland mean opacities, and Σ / = Σ / τ P σT τ P + τ R τ P , (6)where T irr is the irradiation temperature at the disk sur-face determined by the stellar and background black-bodyirradiation as T = T + F irr ( r ) σ , (7)where T bg is the uniform background temperature (in ourmodel set to the initial temperature of the natal cloud core)and F irr ( r ) is the radiation flux (energy per unit time perunit surface area) absorbed by the disk surface at radialdistance r from the central object. The latter quantity iscalculated as F irr ( r ) = L ∗ πr cos γ irr , (8)where γ irr is the incidence angle of radiation arriving atthe disk surface at radial distance r . The incidence an-gle is calculated using the disk surface curvature inferredfrom the radial profile of the disk vertical scale height(Vorobyov & Basu, 2010). The total stellar luminosity L ∗ includes contributions from the accretion and photosphericluminosities. Similar forms of the disk heating term wereemployed in other 1D axisymmetric and 2D thin-disk simulations (e.g. Rice & Armitage, 2009; Zhu et al., 2012;Bae et al., 2014).
3. Improving the thermal balance calculations: the2+1D approach
The numerical method described above is a fast and conve-nient means for computing the disk evolution with high nu-merical resolution and for many physical realizations (e.g.Vorobyov & Basu, 2010; Zhu et al., 2012). However, thismethod is essentially two-dimensional and, as such, it lacksthe information on the disk vertical structure. Some studies(e.g. Dong et al., 2016) assume a Gaussian or exponentialvertical density profile, but the same cannot be easily donefor the vertical temperature distribution. We therefore havedeveloped a straightforward modification to this method,which enables a calculation of the density and temperaturedistributions in the vertical direction concurrently with thecomputations of the gas dynamics in disk plane. In this section, we formulate the modifications made to thethin-disk model in order to improve the thermal balancecalculations in the disk. These modifications also enable areconstruction of the disk vertical structure, thus providinginformation on the disk volumetric density and temperaturedistributions. The new equations of mass, momentum, andenergy transport now read as follows: ∂ Σ ∂t = −∇ · (Σ v ) , (9) ∂∂t (Σ v ) + ∇ · (Σ v ⊗ v ) = −∇P + Σ g + ∇ · Π , (10) ∂e∂t + ∇ · ( e v ) = −P ( ∇ · v ) + ( ∇ v ) : Π . (11)While Equations (9) and (10) remain essentially sim-ilar to their thin-disk counterparts (apart from the effectof stellar motion), Equation (11) now updates the internalenergy only due to advection, viscous dissipation, and pres-sure work (adiabatic heating and cooling). To take the diskheating by the stellar and background irradiation and thedisk cooling due to its own infrared emission into account,we solve the moment equations describing the propagationof diffuse IR radiation in the vertical direction written inthe Eddington approximation: c V ∂T∂t = κ P c ( E − aT ) + S (12) ∂E∂t − ∂∂z (cid:18) c ρκ R ∂E∂z (cid:19) = − ρκ P c ( E − aT ) , (13)where E is the radiation energy density, T the gas tem-perature, ρ the gas volume density, c V the heat capacityof the gas, c the speed of light, a the radiation constant, z the vertical distance from the midplane, σ the column den-sity measured from the disk mid-plane, and S the heatingsource (per unit mass) by the stellar and interstellar UVradiation. though some modifications include a calculation of the diskvertical scale height and the incidence angle of stellar irradiation(Vorobyov & Basu, 2010). Vorobyov & Pavlyuchenkov: The 2+1D model for disk evolution Equations (12) and (13) are complemented with theequation describing the local vertical hydrostatic balancein the disk taking into account the gravity of the star aswell as the local self-gravity of the disk:
Rµρ d ( ρT ) dz = − GM ∗ r z − πGσ, (14)where M ∗ is the mass of the central star and µ = 2 . dσ = ρdz . Weassume that the disk vertical columns at various positionsin the disk do not influence each other, so that solving forEquations (12)-(14) reduces to a series of 1D problems foreach grid zone on the ( r, φ ) computational mesh. Our solution method for Equations (9)-(14) consists ofthree steps. In the first, so-called source step, we updatethe gas velocity and internal energy (per unit surface area)due to gravity, viscosity and pressure work by solving forthe following equations ∂∂t (Σ v ) = −∇P + Σ g + ∇ · Π , (15) ∂e∂t = −P ( ∇ · v ) + ( ∇ v ) : Π . (16)In the second, so-called ”thermal” step, we compute thechange in the disk gas temperature due to radiative cool-ing/heating and reconstruct the disk vertical structure. Todo that, we solve for the moment equations (12) and (13)describing the propagation of diffuse IR radiation in thevertical direction complemented with equation (14) for thevertical hydrostatic balance. We note that in Step 2 we usethe gas temperatures that are partly updated in Step 1.The detailed solution procedure for Step 2 is provided inthe Appendix.In the third, so-called transport step, we take advectionof hydrodynamical quantities into account by solving forthe following equations: ∂ Σ ∂t + ∇ · (Σ v ) = 0 , (17) ∂∂t (Σ v ) + ∇ · (Σ v ⊗ v ) = 0 (18) ∂e∂t + ∇ · ( e v ) = 0 . (19)In this final step, we use the gas velocities and internalenergies which are consequently updated during Steps 1and 2.To accomplish steps 1 and 3, we employ a combina-tion of the finite difference and finite volume methods witha time-explicit solution procedure similar in methodologyto the ZEUS code (Stone & Norman, 1992). The advectionin step 3 is treated using the third-order-accurate piece-wise parabolic scheme of Colella & Woodward (1984). Asmall amount of artificial viscosity is added to the code tosmooth out shocks over two grid zones in both coordinatedirections. The associated artificial viscosity torques inte-grated over the disk area are negligible in comparison withgravitational or turbulent viscosity torques. Equations (15), (16), and (17)-(19) are discretized in po-lar coordinates ( r, φ ) on a numerical grid with 512 ×
512 gridzones. The radial points are logarithmically spaced, whilethe azimuthal points are equidistant. The innermost gridpoint is located at the position of the sink cell r sc = 10 AU,and the size of the first adjacent cell is 0.14 AU which cor-responds to a radial resolution of ∆ r = 1 . R J = h v i /πG Σ,where h v i is the velocity dispersion in the disk plane(Vorobyov, 2013), is resolved by roughly 10–20 grid zones ineach coordinate direction to radial distance up to 500 AU,thus fulfilling the Truelove criterion (Truelove et al., 1998).Equations (12)–(14) are solved on an adaptive, non-equidistant grid with 32 grid points. The finest grid spacingis usually in the disk atmosphere where the largest densitygradients are found. The inner and outer boundaries on thepolar grid ( r, φ ) allow for material to freely flow out fromthe active computational domain, but prevent any materialto flow in. In the vertical direction, we assume a reflectingboundary in the disk midplane and a constant gas volumedensity of 10 cm − at the disk upper edge. Our numerical simulations start from a pre-stellar core withthe radial profiles of column density Σ and angular velocityΩ described as follows:Σ( r ) = r Σ p r + r , (20)Ω( r ) = 2Ω (cid:16) r r (cid:17) s (cid:18) rr (cid:19) − , (21)where Σ and Ω are the gas surface density and angularvelocity at the center of the core. These profiles have a smallnear-uniform central region of size r and then transitionto an r − profile; they are representative of a wide classof observations and theoretical models (Andr´e et al., 1993;Dapp & Basu, 2009).Our iterative solution procedure for Equations (12)–(14)(see the Appendix), requires us to make an initial guess forthe vertical structure of the core. We assume a constant gastemperature of 15 K and a Gaussian distribution of the gasvolume density of the form ρ ( z ) = ρ e − ( z/h ) , (22)where ρ = Σ / ( h √ π ). We note, however, that these ini-tial conditions is a mere initialization requirement and thevertical structure of the core quickly attains the form de-termined by the combined action of the external heating,radiative cooling, pressure gradients, and self-gravity of thecore.We have considered several model cores, but in this pa-per, for the sake of conciseness, we present only two. Theinitial parameters of these models are shown in Table 1,where M c is the initial core mass, β the ratio of rotationalenergy to the magnitude of gravitational potential energy,and r out the initial radius of the core. The two models aredifferent in the amount of initial rotation, as manifested bydistinct β and Ω . The initial parameters are chosen so thatthe cores are gravitationally unstable and begin to collapseat the onset of numerical simulations. We monitor the gas orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 5 Table 1.
Model core parameters
Model M core β Ω r Σ R out ( M ⊙ ) (%) (km s − pc − ) (AU) (g cm − ) (pc)A 1.38 0.5 1.53 2057 9 . × − . × − surface density in the sink cell, and when its value exceedsa critical one for the transition from isothermal to adia-batic evolution, we introduce a central point-mass star. Inthe subsequent evolution, 90% of the gas that crosses theinner boundary is assumed to land onto the growing star.The other 10% of the accreted gas is assumed to be carriedaway with protostellar jets The simulations continue intothe embedded phase of star formation, during which a pro-tostellar disk is formed. The simulations are terminated inthe T Tauri phase when nearly all material of the parentalcore has accreted onto the resulting star-plus-disk system.
4. Comparison of 2+1D and 2D approaches
In this section, we compare the properties of circumstellardisks obtained using the original 2D thin-disk model andthe improved 2+1D model with the purpose to understandhow an improved calculation of the thermal balance in thedisk can effect its dynamical evolution. First, we computethe disk evolution using the 2+1D models starting from thecollapse of pre-stellar cores and ending in the T Tauri stageof disk evolution when the disk age exceeds 0.5 Myr. Then,we compute the disk evolution in the 2D model startingfrom a certain time instance, usually, in the Class I phase,using the current values of Σ, e , and v taken from the 2+1Dsimulations. This allows us to directly determine the effectof different thermodynamical schemes on the disk dynami-cal evolution and avoid the possible dynamical influence ofthe early collapse phase. In the following text, the modelscomputed using the 2+1D approach are distinguished withthe ”plus” sign. Figure 1 presents the disk evolution in model A+ computedusing the 2+1D approach. Shown are the gas surface den-sity snapshots in the inner 2000 × box taken atseveral times elapsed since the onset of numerical simula-tions. The entire numerical domain is about 10 times largerand includes the infalling envelope. Evidently, the disk isstrongly gravitationally unstable and multiple fragmentsare seen forming in the disk’s middle and outer regions. TheToomre Q-parameter, Q = c s Ω /πG Σ, is smaller than unityin and around the fragments, as is demonstrated by theyellow contour lines in the inserts of Figure 1. The massesof these fragments range from about that of a Jupiter tothe upper-mass brown dwarfs.This behaviour is typical for massive, non-magnetizedor weakly-magnetized disks in the embedded phase ofevolution (especially, in the Class I phase), where grav-itational instability is fueled by continuous mass loadingfrom the infalling parental cloud (Vorobyov & Basu, 2005,2015; Kratter et al., 2008; Tsukamoto et al., 2013; Meyer,2017). In model A+, the embedded phase ends around t = 0 .
22 Myr, but fragmentation continues into a later phase because the disk is massive with the disk-to-star massratio ∼ . × box. The 2D simulations startfrom t = 0 .
25 Myr, explaining why the disk appearance at t ≤ .
25 Myr is identical in models A+ and A. At latertimes, however, notable differences in the disk evolutionbetween model A+ and model A arise. While the disk inmodel A+ is strongly unstable and often harbours multiplefragments, the disk in model A experiences only occasionalfragmentation after t = 0 .
25 Myr.This difference is further illustrated in Figure 3 showingthe number of fragments in the disk N frag at a certain timeinstance as a function of time elapsed since the onset of nu-merical simulations. Two conditions were used to identifythe fragments in the disk (see Vorobyov, 2013, for detail):1) they must be pressure supported, with a negative pres-sure gradient with respect to the center of the fragment;and 2) they must be kept together by gravity, with the po-tential well being deepest at the center of the fragment.Evidently, the number of fragments in the disk at a giventime instance is greater in model A+ than in model A. Inaddition, the duration of the disk fragmentation stage (de-fined as the time span during which fragments are presentin the disk) is longer in model A+. Both evidence suggestthat the disk in model A+ is more gravitationally unstable.This difference in the strength of gravitational insta-bility can in principle be caused by distinct disk and stel-lar masses in models A+ and A. It is known that systemswith a higher disk-to-star mass ratio are characterized bystronger gravitational instability because they have on av-erage lower angular velocities and/or higher surface den-sities, and are thus characterized by lower Q -parameters.However, we found that model A+ has a slightly lower diskmass and a slightly higher stellar mass than model A, mean-ing that the mass transport rate through the disk in thismodel is systematically higher than in model A. The highermass transport rates in model A+ are caused by systemat- Vorobyov & Pavlyuchenkov: The 2+1D model for disk evolution
Fig. 1.
Gas surface density distributions in model A+ shown at several time instances after the onset of numericalsimulations (the disk forms at t = 0 . Q -parameter smaller than unity. The red dashed line is a radial cut through thedisk used later to show the vertical volume density and temperature distributions. The scale bar is in log g cm − . Fig. 2.
Similar to Figure 1, but for model A.ically higher gravitational torques, as can be expected fordisks with stronger gravitational instability.The stronger gravitational instability in model A+ ascompared to that in model A is not related to the disk orstellar masses. It may then be related by differences in thedisk thermal structure arising from distinct approaches to calculating the disk cooling and heating in 2+1D and 2Dapproaches. To check if this is indeed the case, we plot inFigure 4 the radial gas temperatures profiles in model A+(red lines) and model A (black lines), calculated as T hydro = e ( γ − µ/ Σ R for every cell on the polar grid ( r, φ ) and thenarithmetically averaged over the azimuth. We note that in orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 7 N u m be r o f f r ag m en t s model A+Time (Myr) model A Fig. 3.
Number of fragments in the disk at a certain timeinstance as a function of time elapsed since the onset of nu-merical simulations in model A+ (top) and model A (bot-tom). The disk forms at t = 0 . T hydro represents the gas temperaturemass-weighted over the vertical column of gas in each cell( r, φ ) (see Equation A.17). For model A+, we also plot thetemperature in the disk midplane T mp (blue lines). Fourtime instances elapsed since the onset of both 2+1D and2D simulations are shown.Evidently, notable differences in the radial gas temper-ature distributions develop with time in models A+ andA. The disk in model A+ is systematically colder than inmodel A, no matter what temperature in model A+ is con-sidered: the hydrodynamic ( T hydro ) or the midplane one( T mp ). On the contrary, the inner envelope is warmer inmodel A+ than in model A. These differences can be un-derstood from the following analysis. Let us first considerthe 2D case and for a moment neglect the heating sourcesdue to viscosity and compressional heating due to P dV work. In the steady-state case, the temperature in the diskand envelope will be controlled by a balance between radia-tive cooling Λ (equation 5) and irradiation and backgroundheating Γ (equation 6), so that the midplane temperaturecan be written as T = T + F irr σ . (23)This steady-state temperature is plotted in Figure 4 withthe dashed black lines. Evidently, it is very close to theactual temperature in model A everywhere except the veryinner parts of the disk where viscous and compressionalheating become substantial and the midplane temperaturerises above the analytically predicted values.In the2+1D case, the midplane temperature can be ex-pressed in the following form if the medium is opticallythick to ultraviolet and infrared radiation T = F irr σ . (24)This expression follows from the assumption that the inci-dent UV flux absorbed by dust in the disk is radiated awayisotropically to the upper and lower hemisphere (see, e.g.,equation 12a in Chiang & Goldreich, 1997).A comparison of Equations (23) and (24) in the regionswhere the input from the background irradiation is much Radial distance (AU)
10 100 1000 10000204060
Radial distance (AU)
10 100 1000 G a s t e m pe r a t u r e ( K ) Radial distance (AU)
10 100 1000 G a s t e m pe r a t u r e ( K ) Radial distance (AU)
10 100 1000204060 t=0.27 Myr t=0.3 Myrt=0.35 Myr t=0.5 Myrmodel Amodel A (analytic)model A+ ( hydrodynamic )model A+ (midplane)
Fig. 4.
Azimuthally averaged radial profiles of the gas tem-perature in model A+ (the red and blue lines) and model A(the black solid lines) taken at several times since onset ofnumerical simulations. In particular, the red lines repre-sent the gas temperature mass-weighted over the verticalcolumn of gas, while the blue lines are the gas temperaturein the disk midplane. The dashed solid lines are the gastemperatures in model A derived for the steady-state caseneglecting compressional and viscous heating. The verticaldotted lines mark the radial positions where the disk joinsthe infalling envelope.smaller than from the stellar one (i.e., in the disk), demon-strates that the midplane temperature in model A+ is ex-pected to be a factor of 2 / ≈ . T = κ F κ R F irr σ , (25)where κ F is the dust opacity in the UV band. Evidently, themidplane temperature in the optically thin limit is higherthan in the optically thick one by a factor of ( κ F /κ R ) / .The optically thin limit is expected to take place in the en-velope and this explains an increase in the gas temperaturein model A+ at distances > T bg = 15 K. At the same time,Equation (25) for F irr = DT (only the interstellar com-ponent, see Appendix A) yields the midplane temperature ≈
14 K for κ F /κ R = 10 , demonstrating that the midplanetemperatures in both models converge to a similar value inthe outer envelope.The consequences of this distinct radial temperaturedistribution is twofold. First, a lower gas temperature Vorobyov & Pavlyuchenkov: The 2+1D model for disk evolution
Log of number density [cm -3 ] 0 50 100 150 200 250 300 350 400 z [ A U ] z [ A U ]
10 20 30 40 50 60 70
Fig. 5.
The distributions of gas volume density (top) andtemperature (bottom) in the ( r, z ) plane taken in model A+at t = 0 .
25 Myr through the radial cut with a position angleof 351 ◦ (shown in Figure 1 with the red dashed line).makes the disk in model A+ more gravitationally unsta-ble, as was already noted above. Second, higher infall ratesfrom the collapsing envelope onto the disk ˙ M infall ∼ c /G ,as implied by a higher gas temperature in the inner enve-lope, drive the disk in model A+ faster to the critical pointat which the disk becomes unstable to fragmentation, thuscreating more fragments in the disk. In other words, thecharacteristic time of mass infall onto the disk (or the diskmass replenishment) M disk / ˙ M infall is shorter in model A+.In Figure 5 we present the two-dimensional distribu-tions of the gas volume density and temperature in the( r, z ) plane obtained in model A+ at t = 0 .
25 Myr by tak-ing a vertical cut with a position angle of 351 ◦ (shown inFigure 1 with the red dashed line). The cut is chosen sothat it passes through an inner spiral arm at r ≈
48 AUand a dense fragment with mass ≈ . M Jup at r ≈
207 AU(shown in the insert of the upper-middle panel in Figure 1).Clearly, the gas volume density is highest in the midplaneand drops down toward the disk atmosphere. At the sametime, the temperature is higher in the upper layers and isminimal in the disk midplane everywhere in the disk, ex-cept for the fragments. This means that the temperature inall parts of the disk, except for the fragments, is controlledby the external radiation (stellar and background). In thefragments, however, the gas is intensively heated by com-pressional and viscous heating, while the radiative diffusionis not fast enough to cool the inner, optically thick partsof the fragment down to low temperatures observed aroundthe fragment. We also note that the fragment at r = 207 AU Fig. 6.
Temperature and gas surface density distributionsin model A+. Various definitions of temperature as de-scribed in the text are shown.is gravitationally bound, which results in a very low diskvertical height at and around this point. The signature ofdisk self-gravity is also seen at r = 230 −
250 AU cover-ing the spiral arc in which the fragment resides and wherethe disk vertical height is also significantly reduced as com-pared with the surroundings. It is therefore possible thatfragments may be deeply hidden in the disk, which wouldmake their detection via observations of the near-infraredlight scattered off the disk surface quite difficult, as was al-ready noted in Dong et al. (2016). Fully three-dimensionalsimulations are needed to determine how far from the mid-plane the fragments can be scattered via gravitational in-teraction with other fragments and spiral arms.The 2D temperature distributions in the disk midplaneand at the surface of the disk are shown for model A+in the upper-left and upper-right panels of Figure 6. Thebottom-left panel presents the gas hydrodynamic temper-ature defined as T hydro = e ( γ − µ/ (Σ R ), in a similarfashion as the midplane temperature in the 2D approach(see Equation 5). The gas surface density distribution isalso shown in the bottom-right panel for convenience. Alldistributions are plotted at t = 0 .
25 Myr.The hydrodynamic and midplane temperatures reflectthe structure of the disk – both are rather low in the outerdisk regions and dense spiral arms, but grow in the innerdisk thanks mainly to increased stellar irradiation. Densefragments heated by compressional and viscous heatingstand out as bright spots in the disk. The hydrodynamictemperature T hydro appears to be higher than the midplanetemperature T mp in the disk outer regions and in the innerenvelope. This trend can also be seen in Figure 4 and is ex-plained by the passively heated (by stellar and interstellarirradiation) nature of the outer disk and envelope – the up-per gas layers are always warmer than the midplane. As aresult, T hydro – a vertical average over the gas column – be-comes higher than T mp . In the inner disk regions, where effi-cient viscous and compressional heating operate in the disk orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 9 Vertical distance (AU) G a s v o l u m e den s i t y ( l og g c m - ) Vertical distance (AU) T e m pe r a t u r e ( K ) clump (r=206.7 AU)inter-arm region (r=100 AU)spiral arm (r=47.6 Au) inner disk (r=19 AU) Fig. 7.
One-dimensional vertical cuts taken at differentpositions of the disk (as indicated in the legend) and show-ing the vertical gas volume density (top) and temperature(bottom) distributions.midplane, the situation may reverse and T mp may becomehigher than T hydro . This trend is also evident in Figure 4.The surface temperature smoothly decreases with radialdistance from the star. The lack of azimuthal variations isa mere consequence of the adopted scheme for calculatingthe incidence angle of radiation onto the disk surface, whichuses an azimuthally averaged disk scale height to calculate γ irr . We note that the surface temperature is higher thanthe midplane temperature almost everywhere in the disk,except for the fragments as discussed in more detail below.Finally, in Figure 7 we show the one-dimensional ver-tical gas volume density and temperature profiles taken atseveral positions in the radial cut (the red dashed line inFigure 1). These positions were chosen so as to show thevertical distributions in various sub-structures that maybe present in the disk, such as spiral arms and fragments.More specifically, the blue lines represent the vertical pro-files taken through the fragment shown in the insert of theupper-middle panel in Figure 1, while the cyan lines aretaken through the spiral arm. The black lines represent thevertical profiles in a regular inter-arm disk region and thered lines are the vertical profiles in the inner disk. We notethat we apply an adaptive grid in the vertical direction thatenables an adequate numerical resolution in the disk regionswith strong volume density gradients and also in the diskatmosphere. Evidently, the fragment is characterized by the highestvolume density amounting to 10 cm − . At the same time,the fragment has the most compact structure with a verti-cal size of just 16 AU, notwithstanding the fact that it isactually located quite far away from the star. The midplanesize of this fragment is ≈
30 AU so that the clump has anellipsoidal shape with the semi-axis ratio of a z : a r = 1 : 2.This scaling is in agreement with the ratio of rotational-to-thermal energy of 44%, implying that the fragment hasa substantial rotational support against gravity (in caseof zero rotation, we would have expected a near-sphericalshape). The fragment is also characterized by a peculiar ver-tical temperature distribution: its midplane temperature ishigher than the temperature in its atmosphere. There is alsoa depression in the temperature profile in between the diskmidplane and atmosphere. A similar, non-regular verticaltemperature distribution can also found in the inner diskregions (the red line). In this case, however, the tempera-ture in the atmosphere is somewhat higher than in the diskmidplane. The non-regular temperature profiles are a directconsequence of the compressional and viscous heating op-erating in the optically thick interiors of the fragment andin the inner disk regions. These heating sources are moreefficient than heating due to stellar irradiation, the latterbeing effectively absorbed by the atmosphere of the frag-ment. This is an important phenomenon, which means thatradiation transfer codes that neglect hydrodynamic heatingsources (such as RADMC-3D) cannot accurately reproducethe temperature structure in the fragments formed via diskgravitational fragmentation and in the inner disk regions(see also Dong et al., 2016). The other two disk elements,the spiral arm and the inter-arm region have a vertical tem-perature distribution typical for passively irradiated disks,though the spiral arm demonstrates a mild increase of thegas temperature towards the midplane, caused probably bymild shock and viscous heating. Model B+ is distinct from model A+ in the amount ofrotation initially present in the parental collapsing core(see Table 1). More specifically, model B+ has a lower ra-tio of rotational to gravitational energy β , which impliesthat model B+ forms a less massive disk than model A+.Following the same procedure as was discussed in the be-ginning of Section 4, we first computed the disk evolutionin model B+ and then we computed model B starting froma time instance t = 0 .
215 Myr, which corresponds to theClass I phase of disk evolution.A comparison of model B+ and model B has revealedtrends that are essentially similar to those discussed in de-tails for model A+ and model A. The disk in model B+is more gravitationally unstable, more prone to fragmenta-tion, and, at the same time, is less massive than in model B.This apparent paradox is explained by a systematicallycolder disk temperatures in model B+ as compared tomodel B. For the sake of conciseness, we do not performan in-depth analysis of models B+ and B in this section,but simply present in Figure 8 the number of fragmentsin the disk N frag at a given time instance as a functionof time elapsed since the onset of numerical simulations.Clearly, N frag is larger in model B+ than in model B. Thenumber of fragments in model B+ can become as high as9 in the early evolutionary phase, whereas in model B the N u m be r o f f r ag m en t s Time (Myr) model B+model B
Fig. 8.
Number of fragments in the disk at a certain timeinstance as a function of time elapsed since the onset of nu-merical simulations in model B+ (top) and model B (bot-tom). The disk forms at t = 0 .
138 Myr.number of fragments never exceeds four. The duration ofthe disk fragmentation phase is also longer in Model B+,which implies a higher detection likelihood of disk fragmen-tation.
5. Model caveats and further improvements
In this section, we discuss the model caveats and furtherimprovements of our 2+1D model.
Disk self-gravity . When reconstructing the disk verti-cal structure, a local value for the gas column density σ was used in the last term of Equation (14). Under this ap-proximation, the vertical gas columns do not affect eachother when solving the equation of the vertical hydrostaticbalance. This may affect the disk vertical structure in theregions dominated by gravity. For instance, a sharp dropin the disk height seen around the fragment in the upperpanel of Figure 5 may be an artifact of the adopted approx-imation. In the future, the model needs to be improved bytaking the non-locality of disk gravity into account. Vertical motions.
Another assumption behind the pre-sented model is the neglection of vertical motions. Namely,we assume that the disk attains a vertical hydrostatic equi-librium on time scales much shorter than the dynamicaltimescale. This approximation is justified if the dynamicaltimescale is longer than any other pertinent timescale, suchas the timescales for diffusion of radiation, propagation ofsound waves, stellar heating, which was shown to be usuallythe case for circumstellar disks in Vorobyov et al. (2014).With the current approach, we cannot model interestingeffects, such as vertical oscillations and surface waves. Thisis, however, the price that we are willing pay for the abilityto follow the disk evolution on much longer timescales thanis currently possible with the full 3D numerical hydrody-namics simulations.
Realistic equation of state . In the current paper, we didnot take into account that the heat capacity c v and theadiabatic index γ of gas depend on the temperature sincethe rotation levels of molecular hydrogen can be excited inthe considered temperature ranges. We did not include thiseffect since it requires a more complicated procedure for the solution of the radiation transfer equations. We plan to dothis in follow-up papers.There are other modifications to our radiative transfermodel which can be implemented in the future. Currently,to calculate heating due to the stellar UV radiation, weevaluate the incidence angle of stellar radiation γ irr (seeEquation A.4) in each grid cell based on the local gradientsof the disk scale height and then average the resulting val-ues over the azimuth (Vorobyov & Basu, 2010). This pro-cedure provides the radially varying UV heating due to theglobal disk flaring, but does not allow us to model prop-erly the local effects, such as shadows in the disk causedby spiral arms. Our attempt to calculate the azimuthallyvarying γ irr , including the effect of the shadows, resultedin overheating of the disk surfaces directly exposed to stel-lar radiation. Allowing for the diffusion of the thermal IRemission in the equatorial plane may alleviate this problem.Finally, we now work on implementing the FARGO mech-anism for advection of hydrodynamical quantities (Masset,2000), which will allow us to ease the restrictive time steplimitations and move the sink cell boundary closer to thestar, hopefully, to the sub-AU region.
6. Conclusions
In this paper, we have improved the frequently used thin-disk models of circumstellar disk evolution and presentedthe method that includes a better calculation of the diskthermal balance and a reconstruction of the disk verticalstructure. Our method is based on the solution of the hy-drodynamics equations describing the gas dynamics in theplane of the disk, complemented with solution of the radia-tion transfer and hydrostatic balance equations describingthe disk vertical structure. We performed a detailed com-parison of this so-called 2+1D method with the purely 2Dthin-disk models of disk evolution. Our findings can be sum-marized as follows. – Improved 2+1D models yield systematically colderdisks, while the infalling parental clouds in the embed-ded evolutionary phase are warmer. Both effects act toincrease the strength of disk gravitational instability in2+1D models as compared to purely 2D simulations. – Disk gravitational fragmentation is more efficient andthe duration of the disk fragmentation phase is longerin 2+1D models, which implies an increased likelihoodfor detecting disk fragmentation in protostellar disks. – The outer disk regions are mostly characterized by apositive vertical temperature gradient, typical for so-called passive circumstellar disks heated mainly by stel-lar and background irradiation. The inner disk regionsusually have a more complex, non-regular vertical tem-perature distribution having local peaks in the midplaneand in the atmosphere separated by a mild depression.The temperature increase in the midplane is caused byefficient viscous and compressional heating operating inthe inner disk. – Fragments forming in the disk via gravitational frag-mentation are also characterized by a double-peakedvertical temperature profile, but unlike the inner diskregions, the center of the fragment is warmer than itsatmosphere. This means that the hydrodynamical heat-ing sources (compression, viscosity) are more efficientthan stellar irradiation heating, implying that radiation orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 11 transfer codes that neglect the hydrodynamical heatingsources cannot accurately compute the temperature inthe fragment interiors (see also Dong et al., 2016). – Fragments are significantly more compact than the sur-rounding disk, which could make their detection withthe scattered light techniques difficult unless significantexcursions away from the disk midplane are feasible.A detailed procedure for solving the radiation transferand hydrostatic balance equations in the vertical directionis presented in the Appendix. The improved 2+1D methodyield a full 3D structure of the disk, namely its volumetricdensity and temperature distributions, with a modest in-crease in the net computational time with respect to purely2D models. It also applies an adaptive mesh in the verticaldirection, enabling good resolution in the disk regions withstrong density gradients and in the disk atmosphere. Thismakes it possible to couple our 2+1D model with compactchemical reaction networks to follow the long-term chemo-dynamical evolution of young protostellar disks startingfrom the gravitational collapse of parental cloud cores. Thedetails of this method will be presented in the upcomingpaper.
Acknowledgments.
We are thankful to the anonymousreferee for constructive comments that helped us to im-prove the manuscript and the vertical reconstruction pro-cedure. This work was supported by the RFBR grant 17-02-00644. The simulations were performed on the ViennaScientific Cluster (VSC-2), on the Shared HierarchicalAcademic Research Computing Network (SHARCNET), onthe Atlantic Computational Excellence Network (ACEnet).This publication is supported by the Austrian Science Fund(FWF).
References
Akimkin, V. V., Pavlyuchenkov, Y. N., Launhardt, R., &Bourke, T. 2012, Astronomy Reports, 56, 915Akimkin, V., Zhukovska, S., Wiebe, D. et al. 2013, ApJ,766, 8ALMA Partnership; Brogan, C. L., Perez, L. M., Hunter,T. R., et al. 2015, ApJ, 808, L3Andr´e, Philippe; Ward-Thompson, Derek; Barsony, Mary1993, ApJ, 406, 122Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ,795, 61Baraffe, I., Elbakyan, V. G., Vorobyov, E. I., Chabrier, G.2017, A&A, 597, 19Basu, S., & Vorobyov E. I. 2012, ApJ, 750, 30Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C.2014, in Protostars and Planets VI, H. Beuther, R.S. Klessen, C. P. Dullemond, and T. Henning (eds.),University of Arizona Press, Tucson, 691Binney, J., & Tremaine, S. 1987, Galactic Dynamics(Princeton, NJ: Princeton Univ. Press)Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H.2010, Icarus, 207, 509Boss A. P. 1998, ApJ, 503, 923Chiang, E., & Goldreich, P. 1997, ApJ, 490, 368Colella, P., & Woodward, P. R. 1984, J. Comput. Phys.,54, 174Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092 Dong, R., Vorobyov, E. I., Pavlyuchenkov, Y., Chiang, E.,& Liu, H. B. 2016, ApJ, 823, 141Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002,A&A, 389, 464Dullemond, C. 2012, Astrophysics Source Code Library,record ascl:1202.015Elbakyan, V. G., Vorobyov, E. I., & Glebova, G. 2016,Astron. Rep., 60, 879Galvagni, M., & Mayer, L. 2014, MNRAS, 437, 2909Gyergyovits, M., Eggl, S., Pilat-Lohinger, E., & Theis, Ch.2014, A&A, 566, 114Henning, Th., & Semenov, D. 2013, Chem. Reviews, 113,9016Hubeny, I. 1990, ApJ, 351, 632Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131Klahr, H., & Kley, W. 2005, A&A, 445, 747Kley, W. & Crida, A. 2008, A&A, 487, L9Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008,ApJ, 681, 375Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning,T. 2010, A&A, 511, 81Lomax, O., Whitworth, A. P., Hubber, D. A., Stamatellos,D., Walch, S. 2014, MNRAS, 439, 3039Machida, M. N., Inutsuka, S., & Matsumoto, T. 2010, ApJ,724, 1006Masunaga, H., Miyama, S. M., & Inutsuka, S. 1998, ApJ,495, 346Mayer, L., Peters, T., Pineda, J. E., Wadsley, J. Rogers, P.2016, ApJ, 823, 36Masset, A. 2000, AAS, 141, 165Meyer, D. M.-A., Vorobyov, E. I., Kuiper, R., Kley, W.2017, MNRAS, 464, 90Nayakshin, S. 2010, MNRAS, 408, 2381Nayakshin, S. 2011, MNRAS, 413, 1462Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65Pavlyuchenkov, Y. N., Zhilkin, A. G., Vorobyov, E. I., &Fateeva, A. M. 2015, Astronomy Reports, 59, 133Reg´aly, Zs., S´andor, Zs., Csom´os, P., Ataiee, S. 2013,MNRAS, 433, 2626Reg´aly, Zs., & Vorobyov, E. I. 2017, A&A, in revisionRice, W. K. M., Armitage, P. J., 2009, MNRAS, 396, 2228Semenov, D., Henning, T., Helling, C., Ilgner, M., &Sedlmayr, E. 2003, A&A, 410, 611Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392,413Stamatellos, D., 2015, ApJL, 810, 11Stamatellos, D., & Herczeg, G. J. 2015, MNRAS, 449, 3432Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753Tobin, J. J., Looney, L.W., Wilner, D.J., Kwon, W.,Chandler, C.J., et al. 2015, ApJ, 805, 125Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1998, ApJ,495, 821Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2013,MNRAS, 436, 1667Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., &Inutsuka, S. 2015, 446, 1175Vorobyov, E. I., & Basu, S. 2005, ApJL, 633, 137Vorobyov, E. I., & Basu, S. 2009, MNRAS, 393, 822Vorobyov, E. I. 2010, ApJ, 723, 1294Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896Vorobyov, E. I. 2011, ApJL, 728, 45Vorobyov, E. I., 2013, A&A, 552, 129
Vorobyov, E. I., Zakhozhay, O. V., & Dunham, M. M. 2013,MNRAS, 433, 3256Vorobyov, E. I., Pavlyuchenkov, Y. N., & Trinkl, P. 2014,Astronomy Reports, 58, 522Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115Vorobyov, E. I. 2016, A&A, 590, 115Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F.2012, ApJ, 746, 110 orobyov & Pavlyuchenkov: The 2+1D model for disk evolution 13
Appendix A: The thermal step and thereconstruction of the disk vertical structure
The thermal step is placed in our algorithm between thesource step (Equations 15 and 16) and the transport step(Equations 17-19). While the source and transport stepsdeal with two-dimensional quantities integrated over thedisk vertical column (such as the gas surface density Σ,the vertically integrated pressure P and the internal energyper surface area e ) and are therefore defined on the ( r, φ )polar mesh, the thermal step deals with three-dimensionalquantities (such as the gas temperature T , the gas volumedensity ρ , the radiative energy E and the vertical columndensity from the disk midplane σ ) defined on the ( z, r, φ )cylindrical mesh. The initial values for all quantities (two-and three-dimensional) are provided by the initial setup asdescribed in Section 3.3.Here, we describe the algorithm for solvingEquations (12) and (13), which take the disk radia-tive cooling/heating into account. This algorithm alsoincludes the reconstruction of the disk vertical structureassuming the vertical hydrostatic equilibrium described byEquation (14). Our method is a time-dependent modifica-tion to the steady-state models presented in Akimkin et al.(2012) and Vorobyov et al. (2014). The method for solvingthe source and transport steps are described in Section 3.2and in more detail in Vorobyov & Basu (2010).Firstly, we note that the internal energy per surface area e is updated in the source step due to the P ( ∇ · v ) workand viscous heating. This may affect the gas temperaturedistribution in the disk and should be taken into accountbefore the thermal step is commenced. We assume that thegenerated (or consumed) heat in the source step is evenlyredistributed over the disk vertical column so that the gastemperature in each element of the vertical column can beupdated as T ∗ = T n e ∗ e n , (A.1)where index n refers to the quantities at the beginning ofthe source step and the asterix – to the quantities after thesource step. The updated gas temperature T ∗ is furtherused in the thermal step.In the second step, we compute the heating function S (Equation 12) due to sources other than the thermal radi-ation of the medium. We note that the dimension of S bydefinition is [erg s − g − ]. In our model, these sources in-clude the UV radiation from the central star S star and theinterstellar UV radiation S bg : S = S star + S bg , (A.2)where S star is defined as S star = 4 πκ starP J star , (A.3)where κ starP is the mean Planck opacity averaged over thestellar spectrum.To compute S star , we need to know the distribution ofthe mean intensity of UV radiation J star . We calculate J star taking into account the absorption of the UV radiation bythe disk as follows: J star = J exp ( − τ star /µ ) , (A.4)where J is the UV intensity at the surface of the disk, τ star is optical depth calculated from the surface of the disk, and µ = cos γ irr the cosine of the incidence angle of stellarradiation arriving at the disk surface with respect to thenormal (see Vorobyov & Basu, 2010, for details). The opti-cal depth as a function of the current column density σ isdefined as τ star = κ starP (cid:18) Σ2 − σ (cid:19) , (A.5)where σ is measured from the disk mid-plane. The intensity J at the disk surface can be found as J = 14 π L ∗ πr . (A.6)Here, L ∗ is the total stellar luminosity which includes con-tributions from photoshperic luminosity L ph and accretionluminosity L accr and r the radial distance to the star. Theaccretion luminosity is calculated as L acc = GM ∗ ˙ M / (2 R ∗ )using information on the current stellar mass M ∗ , stel-lar radius R ∗ , and accretion rate onto the star ˙ M . ThePhotospheric luminosity and stellar radius are found fromthe Lyon stellar evolution code coupled to the main hy-drodynamics code as described in Vorobyov & Basu (2015)and Baraffe et al. (2017).The background heating function S bg is calculated byanalogy to Equation (A.4), but with a fixed value of µ =0 .
5. As a boundary condition for the interstellar radiation,we use the following relation: J = D caT π , (A.7)where c is the speed of light, a the radiative constant, T bg and D the temperature and dilution of the interstellar radi-ation. In our simulations, we adopt 10 000 K and 8 × − ,correspondingly.In the third step, we calculate the change of the gastemperature T and radiative energy E in a given verti-cal column of the disk due to heating by the stellar UVand background radiation and cooling by the disk infraredemission. The thermal evolution of the vertical column isdescribed by the system of radiative transfer Equations (12)and (13).This system is a set of moment equations for diffuse IRradiation derived under the Eddington approximation. Wesolve equations (12) and (13) numerically to find T and E after one hydrodynamical time step ∆ t using the followingimplicit finite-difference scheme: c V T − T ∗ ∆ t = κ P c ( E − aT ) + S (A.8) E − E n ∆ t − ˆΛ E = − ρ n κ P c ( E − aT ) , (A.9)where T ∗ is the gas temperature updated after the sourcestep (see Equation A.1), and E n and ρ n the radiation en-ergy and gas volume density taken from the previous timestep. In the discretized equations, we use the following con-vention: index i refers to the left-hand-side grid interfaceand index i + 1 / i +1 / E denotes the finite-differenceapproximation to the diffusion term:ˆΛ E = 1∆ z i +1 / (cid:18) c ρ ni +1 κ R ,i +1 E i +3 / − E i +1 / ∆ z i +1 − c ρ ni κ R ,i E i +1 / − E i − / ∆ z i (cid:19) , (A.10) where∆ z i = z i +1 / − z i − / ∆ z i +1 / = z i +1 − z i ρ ni κ R ,i = 12 (cid:16) ρ ni − / κ R .i − / + ρ ni +1 / κ R ,i +1 / (cid:17) . The above non-linear system is solved with the iter-ative Newton method. Namely, we approximate T as( T k ) (cid:18) TT k − (cid:19) using the first two terms of the Taylorexpansion series, where superscript k refers to the previousiterative step and produce a system of linear equations for E / , ..., E M − / with a tridiagonal matrix, where M is thenumber of grid cells in the vertical direction. This tridiag-onal system is solved with the forward and back substitu-tion method (the Thomas algorithm). The obtained valuesof energies are back substituted in Equation (A.8) to de-rive the temperatures and both quantities are then used forthe next Newton iteration as T k +1 and E k +1 until conver-gence is achieved. The solution is then set to T = T k +1 and E = E k +1 .The boundary conditions for Equations (A.8) and (A.9)are the zero energy gradient near the mid-plane and follow-ing relation between the diffusion flux and radiation energyat the disk surface: − c ρ n κ R ∂E∂z (cid:12)(cid:12)(cid:12)(cid:12) z M = c (cid:0) E − aT (cid:1) , (A.11)where T cmb =2.73 K. The latter condition is derived underthe assumption that the incoming and outgoing radiationis isotropic over the upper and lower hemispheres. We notethat the presented method is conceptually similar to thatwhich is developed for the calculation of the pre-stellar corethermal evolution in Pavlyuchenkov et al. (2015).The next stage of the algorithm is to recover the hy-drostatic equilibrium along the vertical direction using theupdated temperatures. Equation (14) describing the verti-cal hydrostatic equilibrium can be written as follows: dzdσ = 1 ρ (A.12) Rµ d ( ρT ) dσ = − GM ∗ r z − πGσ, (A.13)where r is the radial distance to a given vertical column ofgas, G the gravitational constant, σ the gas column densitymeasured from the disk mid-plane. We note that T is knownfrom the previous step. Equation (A.13) accounts for thegravity of the central star and disk self-gravity in the plane-parallel limit. The boundary conditions for this system havethe following form: z (0) = 0 (A.14) ρ ( σ M ) = ρ ext , (A.15)where ρ ext is the gas volume density at the disk surface. Inour models, ρ ext corresponds to a hydrogen number densityof 10 cm − .The solution of the system (A.12) — (A.13) with thecorresponding boundary conditions can be found using animplicit scheme similar to that used when computing the internal structure of the stars. We linearize the right-handside of (A.12) as follows:1 ρ ≈ ρ k − ρ k ) ( ρ − ρ k ) , (A.16)which transforms the initial system into a linear sys-tem of ordinary differential equations. The implicit finite-difference approximation of the latter one generates a sys-tem of algebraic linear equations with a tridiagonal matrixwhose solution can be easily found using the forward andbackward substitution method. The resulting solution isused to form a new approximation, and iterations over k are carried out until convergence is achieved (usually, aftera few iterations).At the last step of the algorithm, we calculate the up-dated thermal energy e per unit area by summing up thethermal energies over the vertical grid: e = 2 c V M X k =0 T i +1 / ( σ i +1 − σ i ) , (A.17)where a factor of 2 accounts for the fact that we adoptedan equatorial symmetry with resect to the disk midplane.This value will be further used in the transport step.Finally, we note that the transport step updates the val-ues of Σ and e (see Equations 17 and 19), which in turnaffects the volumetric distributions of the gas temperature T and volume density ρ . To take these changes into ac-count, we assume that the mass and thermal energy thatare carried with the flow are evenly redistributed over thevertical cells so that the updated distributions of T and ρ can be calculated as: ρ n +1 = ρ Σ n +1 Σ , (A.18) T n +1 = T ΣΣ n +1 e n +1 e , (A.19)where T , ρ , and Σ are the values of the gas temperature, itsvolume and surface densities before the transport step. Thiscompletes one cycle of integration and the updated values T n +1 , ρ n +1 , and Σ n +1 are used at the next time step.Since the method is fully implicit, it is in general verystable. However, in very few cases the Newton proceduremay not converge for some values of the time step. In thiscase, we divide the hydrodynamic time step into severalsub-cycles, which usually solves the problem. In the limitof large time steps, the method yields a steady-state so-lution similar to that obtained by the steady-state modeldescribed in Vorobyov et al. (2014).The presented algorithm was carefully tested and com-pared with other methods. In particular, we benchmark oursteady-state solutions (in the limit of large time steps) withthe results of 1+1D code “diskstruct” developed by C.Dullemond and discussed in Dullemond et al. (2002). Theresults of this comparison are shown in Figure A.1. The ver-tical temperature distributions calculated with our methodturn out to be very similar to the results obtained withthe approximate method “vertrt”, which also adopted greyopacities for the dust thermal radiation included in the“diskstruct” package. The deviation of our solution from ∼ dullemond/software/diskstruct/index.shtmlorobyov & Pavlyuchenkov: The 2+1D model for disk evolution 15
10 100 1000 0 0.04 0.08 0.12 0.16 0.2 0.24 T , K z, AU R = 1 AU 10 100 1000 0 1 2 3 4 5 T , K z, AU R = 10 AU 10 100 1000 0 10 20 30 40 50 60 70 80 90 T , K z, AU R = 100 AUdiskstruct-fullvrtdiskstruct-vertrtour model Fig. A.1.
Vertical temperature distributions of a proto-planetary disk calculated with our method and with the“diskstruct” package written by C. Dullemond. More specif-ically, the results from the “diskstruct-fullrt” method withthe full treatment of the frequency and angular dependentradiative transfer are shown with the thick blue lines. Thethick red curves represent the results from the approximate“diskstruct-vertrt” method, which adopts grey opacities forthe dust thermal radiation. The temperature distributionscalculated with our method are shown with the thin blacklines. The presented vertical distributions are calculated forthree radial positions in the disk, namely for R = 1 AU(top), R = 10 AU (middle), and R = 100 AU (bottom) foran optically thick protoplanetary disk illuminated by a starwith T eff = 6000 K.the results produced by the more accurate method “fullrt”,which takes into account the frequency and angular depen-dence of radiative transfer, is notable near the equatorialplane. However, we consider these differences to be accept-able for our simulations.To illustrate the time-dependent aspect of the radia-tive transfer model used in the thermal step, we show inFigure A.2 the evolution of the temperature distributionat r = 10 AU with the density distribution manually fixed.For the parameters of the considered vertical column of gas,
0 0.5 1 1.5 2 2.5 3 T , K z, AU T star =4000KR=10AU Σ g =10 g/cm t=1e9t=1e8t=1e7t=1e6t=1e5t=0 Fig. A.2.
Evolution of the gas temperature distribution at r = 10 AU in a toy model of the protoplanetary disk. Thedensity structure is fixed. The initial temperature distribu-tion is uniform with T = 4 K (first case) and T = 400 K(second case). The initial distributions are shown with theshades of blue, the steady-state solution is shown with thered color. The arrows illustrate the evolution of the gastemperature profiles with time. The corresponding time inseconds is shown in the color legend. -16 -15 -14 -13 -12 -11 ρ , g / c m z, AU T star =4000KR=10AU Σ g =10 g/cm t=1e9t=1e8t=1e7t=1e6t=1e5t=0 Fig. A.3.
Evolution of the gas volume density distributionat r = 10 AU in a model of protoplanetary disk. The initialtemperature distribution is uniform with T = 4 K (firstcase) and T = 400 K (second case). The initial density dis-tributions are shown with the shades of blue, the steady-state solution is shown with the red color. The arrows il-lustrate the evolution of the gas density profiles with time.The corresponding time in seconds is shown in the colorlegend.the relaxation to a steady-state solution is achieved withinabout a few years (about an order of magnitude faster thanthe dynamical time at r = 10 AU) starting from either awarmed-up disk (with an initial uniform temperature of T = 400 K) or from a cooled-down disk (with a uniformtemperature of T = 4 K). This relaxation time turns outto be in agreement with the characteristic thermal timeestimated in Vorobyov et al. (2014). As expected, the re-laxation to the steady-state solution is faster in the upperlayers of the disk where the gas volume density is lower. While producing the temperature distributions inFigure A.2, the gas volume density distribution was man-ually fixed. When the density is allowed to evolve togetherwith the temperature, then the density distribution devel-ops as shown in Figure A.3. We see that at low and highinitial temperatures the disk vertical height is low and high,respectively. In the steady-state profile near z = 1 ..