1+1D implicit disk computations
Florian Ragossnig, Ernst A. Dorfi, Bernhard Ratschiner, Lukas Gehrig, Daniel Steiner, Alexander Stökl, Colin P. Johnstone
11 +
1D implicit disk computations
Florian Ragossnig a, ∗ , Ernst Dorfi a , Bernhard Ratschiner a , Lukas Gehrig a , Daniel Steiner a ,Alexander St¨okl a , Colin P. Johnstone a a University of Vienna, Department of Astrophysics, T¨urkenschanzstr. 17, 1180 Vienna, Austria
Abstract
We present an implicit numerical method to solve the time-dependent equations of radiationhydrodynamics (RHD) in axial symmetry assuming hydrostatic equilibrium perpendicular tothe equatorial plane (1 + Keywords:
Numerical methods, partial di ff erential equations, radiation hydrodynamics,adaptive grid
1. Introduction
During the gravitational collapse of tenuous interstellar gas towards a condensed object likea protostar, the initial angular momentum has to be reduced by several orders of magnitude toovercome the centrifugal barrier [1, 2]. A solution to this angular momentum problem is theformation of an accretion disk around the central star where turbulence-induced viscosity and / ormagnetic fields provide a mechanism to transport angular momentum outwards accompanied bya mass transport towards the interior regions.We point out that accretion disks do not only occour as protoplanetary disks (PPDs) but alsoaround other astrophysical objects like black holes and white dwarfs. Such accretion disks areusually characterized by large disk masses, high surface irradiation and / or the disks are cut o ff due to the Roche lobe of a companion star feeding the disk [e.g. 3]. Since non-PPDs are eithermassive (disks around black holes) or experience a mass flux over the outer boundary (dwarf-novae), they are unlikely to be dynamically stable. Although, we are confident that our modelcan be adopted to such disks, we want to emphasize that our current approach of an isothermal,vertical disk-structure is approximatevely valid for PPDs. However, the vertical structure ofdisks around x-ray binaries and white dwarfs requires a more careful treatment of the verticalradiative transport, in particular to get the correct, opacity-dependent photospheric temperature ∗ Corresponding author.
E-mail address: fl[email protected]
Preprint submitted to Computer Physics Communications June 24, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J un tratification. We accentuate that this work focuses on protoplanetary accretion disks with diskmasses of less than 10% of the stellar mass and an outer boundary located in the interstellarmedium ISM (cf. Sec. 2.1).In the simplest situation, we can consider a disk with negligible mass around a central objectwith mass M ∗ in centrifugal balance and obtain a Keplerian angular velocity at a radius r Ω = (cid:114) G M ∗ r , (1)with an orbital velocity of u ϕ = Ω r ∝ r − / . Hence, assuming an orbital ratio of 100, the cor-responding orbital periods vary according to Kepler’s 3 rd law by a factor of 1000. This basicphysical requirement immediately hinder all numerical simulations of disk-like structures. Tak-ing ∆ x as an orbital distance, any explicit numerical computation is restricted by the Courant,Friedrichs and Lewy (CFL) condition [4] δ t ≤ min all cells ∆ x | u ϕ | + c s (cid:39) min IB ∆ x | u ϕ | . (2)For stability, the smallest timesteps δ t at the inner boundary (IB) must be taken. Although it ispossible to increase the timestep of an explicit scheme beyond the CFL limit by special methods(e.g. FARGO, local timestepping, [5, 6]), long term evolution studies of PPDs with an innerboundary within a few stellar radii away from the star, are only possible by utilizing implicitmethods. All disks are dominated by centrifugal forces, and thus the sound velocity of the gas, c s , is small compared to the orbital velocity, u ϕ (i.e. c s (cid:28) u ϕ ).Almost all astrophysical disk simulations are carried out by explicit numerical schemes andas a result of the CFL condition, the size of the computational domain towards the inner regionsof the disk is limited [e.g. 7]. For matter orbiting the host star, Keplerian rotation leads to verylarge velocity gradients which results in very small timesteps, typically of the order of one orbitalperiod of the inner region divided by the number of azimuthal grid points. Taking into accountinfrared-observations [e.g. 8] of protoplanetary disks, the estimated lifespan of a disk can beestimated to be up to 10 years. Therefore, the application of explicit schemes is limited, if onewants to investigate the long term evolution of disks. To overcome this restriction, we utilize animplicit time integration scheme which allows us to extend the inner boundary to the surface ofthe host star and is not restricted by timestep length [e.g. 9]. For physical reasons we usuallyset the inner boundary to the co-rotation radius, which is where the rotational velocity of the starmatches the angular velocity of the disk gas (see Sec. 2.1) and is typically at few stellar radii.We emphasize the importance of the inner boundary condition as the solution of disk evolution isstrongly coupled to changes of flux and temperature within the innermost disk. Recent studies ofviscous disk models support this argument that an appropriate description of the innermost regionis a significant factor for global disk evolution [10] and cannot be approximated by speciallydesigned inner sink cells.Usually, boundary conditions for partial di ff erential equations (PDE) determine the solutionout of a class of functions and are therefore more important and more complicated to imple-ment compared to ordinary di ff erential equations. In particular, computing disk-like structurescan strongly depend on the location of the boundaries because neither the inner nor the outerboundary is clearly determined by geometrical conditions. Consequently, the disk mass dependsentirely on the choice of the inner and outer boundary condition.2urthermore, the time step in implicit schemes does not depend on cell size but is entirelydetermined by the accuracy and convergence of the utilized method. If short timescale phenom-ena are present, implicit schemes bear no advantage over explicit schemes. Utilizing an implicitmethod and implementing axial symmetry can lead to time steps of several orbital periods andallows us to study the evolution of interplanetary disks throughout the entire disk lifetime.Although our method adopts axial symmetry and hence no azimuthal resolution can be ob-tained, recent ALMA (Atacama Large Millimeter Array) observations show that 75% of the disksaround young T Tauri stars are axially symmetric protoplanetary disks, at least within the reso-lution limit of the instrument (approximately 5 AU) [11]. Disks showing a spiral-like structuresor small anti-symmetric features either involve a companion star or are most likely disturbed byan undetected planet [12, 13]. Furthermore, recent simulations show that especially less mas-sive disks tend towards axial symmetry shortly after the initial collapse and remain about 80%of their life in an axisymmetric state [see e.g. 10]. Encouraged by such observations and non-axisymmetric simulations, we adopt a cylindrical geometry to study the long term evolution ofprotoplanetary disks.This paper describes recent developments of our initially spherical symmetric
TAPIR code(The Adaptive Implicit RHD Code) [14], now modified for cylindrical configurations and inparticular the geometrical changes necessary to ensure also a conservative formulation of theequations of radiation hydrodynamics on an adaptive grid. Furthermore, self-gravity of the diskmaterial has also been included by solving the corresponding Poisson equation for an axisym-metric density distribution. In Section 2, we describe the physical equations implemented in ourmodel. In Section 3, we present the numerical methods. Section 4 is devoted to numerical ex-amples and accuracy tests to validate our computations. Finally, in Section 5, we discuss somesimulations concerning their astrophysical perspectives.
2. Physical equations
We adopt a cylindrical coordinate system with ( r , ϕ, z ) and assume hydrostatic equilibriumin the z -direction and therefore no vertical velocities. The velocity vector is therefore u = ( u r ( r , t ) , u ϕ ( r , t ) , c s perpendicular to the equatorial plane. To eliminate fur-ther vertical dependences (e.g. in the equation of motion) we adopt the thin disk approximation( z (cid:28) r ), which in particular is given by ( r + z ) − / ∼ r − and results in a local constant verticalgravitational acceleration. Since all quantities are independent on the azimuthal angle ϕ , we canuse vertically integrated quantities. For example, the surface density Σ is related to the massdensity ρ by Σ ( r , t ) = (cid:90) ∞−∞ ρ ( r , z , t ) dz . (3)Correspondingly, the equation of continuity is given by [e.g. 15] ∂ Σ ∂ t + r ∂∂ r ( ru r Σ ) = . (4)As a result of utilizing cylindrical geometry, the equation of motion appears in a radial and3n angular part. Applying the above assumptions to the radial equation of motion gives ∂∂ t ( Σ u r ) + r ∂∂ r (cid:16) r Σ u r (cid:17) − ρ u ϕ r + ∂ P ∂ r + Σ ∂ Ψ tot ∂ r − π c κ R Σ H r + r ∂∂ r ( r Q rr ) − Q ϕϕ r = , (5)where κ R is the Rosseland-mean opacity [e.g. 16], the thermodynamic quantity P represents thevertically integrated gas pressure and is retrieved from tables, H r represents the radial radiationflux component, i.e. 1 st moment of the specific intensity and Q rr and Q ϕϕ are the radial andangular components of the viscous pressure tensor, respectively [17]. The generalised form of Q writes (written as unity matrix) Q = µ Q (cid:34) ( ∇ u ) − ∇ · u (cid:35) (6)The quantity ∂ Ψ tot /∂ r = ∂ Ψ s /∂ r + ∂ Ψ d /∂ r represents the total gravitational potential, where ∂ Ψ s /∂ r = GM ∗ / r is the star’s potential acceleration (cf Sect. 2.3).Since thermal, gravitational and radiation pressures do not contribute to the ϕ -component ofthe equation of motion, the appropriate terms can be neglected. After some simplifications, theangular component of the equation of motion becomes ∂∂ t ( r Σ u ϕ ) + r ∂∂ r (cid:16) r Σ u r u ϕ (cid:17) + ∂∂ r (cid:16) r Q r ϕ (cid:17) + Q ϕ r = . (7)Hydrostatic equilibirum in the z -direction leads to a simple vertial density structure given by[e.g. 18] ρ ( r , z ) = ρ ( r ) exp − z H , (8)where ρ ( r ) is the equatorial density. The last expression uses a defintion of scale height, H p ,given by H p = c g z , (9)and note, that c s and the vertical gravitational acceleration g z depend on the radial distance r .Since we adopt the thin disk approximation [e.g. 18, page 39] (valid for z (cid:28) r ) g z = GM ∗ zr . (10)A more elaborated treatment involves the solution of the Poisson equation (see Sect. 2.3).The generalized form of the equation of internal energy in cylindrical geometry is ∂∂ t ( Σ e r ) + r ∂∂ r ( ru r Σ e r ) + P r ∂∂ r ( ru r Σ ) + (cid:15) Q − πκ P ( J − S ) + ∆ E rad + r ∂∂ r (cid:16) ru φ (cid:17) = , (11)where (cid:15) Q = Q : ∇ u (where : is the tensor product) denotes the viscous energy dissipation [17], κ P is the Planck opacity [e.g. 16] and ∆ E rad is the radiative heating / cooling term which will be4iscussed in Sect. 3 in more detail. Since the radial velocity u r is much smaller than the angularvelocity u ϕ , we can truncate the radiation equations to their dominant terms κ P Σ ( J − S ) = −∇ · H (12)and κ R Σ H = −∇ · K ∼ −∇ ( f edd J ) , (13)where f edd represents the Eddington factor [e.g. 16]. We emphasize that these dominant termsare purely spatial and will not obey any time scale. As the time-dependent terms in the radiationflux equation are usually not important compared to the pure radiation di ff usion terms, they canbe neglected. Implementing Eq. 12 and 13 into Eq. 11, we get ∂∂ t ( Σ e r ) + r ∂∂ r ( ru r Σ e r ) + P r ∂∂ r ( ru r Σ ) + (cid:15) Q + π r ∂∂ r ( rH r ) + ∆ E rad + r ∂∂ r (cid:16) ru φ (cid:17) = . (14)For a consistent description of the RHD equations, we would have to formulate the equationof the radiation flux and the equation for the radiation energy as individual equations. Since theinternal energy budget is dominated by irradiation and radiative cooling, we can simplify thedescription for the radiation temperature J and the radiative flux H further, and implement bothdirectly into the internal energy equation. A more detailed description of how this problem hasbeen tackled is given in Sect. 3.1. As known from the theory of PDEs, boundary conditions will determine the solution to aparticular problem. In the case of astrophysical disks, the locations and the physical propertiesof the boundaries are di ffi cult to specify. Since our computations are restricted to axisymmetricproblems, we have to avoid massive disks where local gravitational instabilities can lead to acollapse of individual gas blobs [7]. In our case, we only assume non-massive disks where M disk (cid:28) M ∗ and where the gravitational potential is dominated by the central star with mass M ∗ .Typically, non-axial structures will grow if M disk exceeds approximately 0 . M ∗ [e.g. 18]. Thedisk’s gravitational potential is smaller by several orders of magnitude (cf. Sect. 2.3 for details).The exact position of the inner and outer boundary of PPDs is hard to define since the innerboundary depends on the physical properties of the star and the outer boundary is basicallydefined at the location of the interstellar medium (ISM). At the inner boundary, a star with afinite radius rotates at a di ff erent angular speed Ω ∗ than the material of the circumstellar disk Ω ( r ). The corotation radius, r co , is defined as the radius where the stellar rotation rate matchesthe orbital speed of the disk, i.e. Ω ∗ = Ω ( r co ), and depends on the topology and the strength of thestellar magnetic field [19]. Due to the centrifugal balance of the disk (cf. Eq. 1), r co determinesthe gap between the star and the disk. Details about the mass flow of disk material onto thestar can only be analyzed by detailed 3D-MHD simulations ([e.g 20]). The role of the stellarmagnetic field emerges also in the so-called magnetic truncation radius r mt [19], where studiesplace the inner boundary of a disk at the location where the stellar magnetic field matches thedisk’s magnetic field [X-point, see 21]. In some cases, depending on the disk mass and in caseof a very weak magnetic field, it is possible that the disk reaches to the surface of the star [22].Regardless of the definite process, the position of the inner boundary for the majority of PPDs is5lose to r co and within a few stellar radii away from the host star. Although our model is capableof setting the inner boundary at any radius (also at the star’s surface), the intention of this workis to present a model to investigate the long term evolution of PPDs. Thus we do not focus ondetailed astrophysical applications but on the importance of the ability of a model to set the innerdisk boundary close to the star. In Sect. 4 we show that the solutions are critically a ff ected ifa simplistic central hole is cut out of the computational domain. Hence we emphasize that thetreatment of the inner boundary is crucial to obtain reasonable physical models and neglectingthe innermost regions can critically distort the disk structure.Since most PPDs do not have a companion and therefore no mass inflow over the outer diskboundary, the outer boundary is less important at distances far from the central source. Notethat the exact position of the outer boundary does become more important when dealing withmassive disks (e.g. around black holes) and / or the disk is getting fed by a companion [3]. Again,examples are given in Sect. 4. In principle, our model allows to define the outer boundary at verylarge radii but to avoid numerical issues, we set the outer disk radius to a value where the denistydoes not fall below ρ out = − g / cm according to the standard disk model [18]. For example, ifthe density drops below the numerical accuracy, we can set it to values typical for the ISM nearthe outer disk edge. Since implicit numerical schemes are not restricted by the CFL timestepcondition, the implementation of appropriate boundary conditions near to the centre as well asfar outside enables large space regions to be covered and therefore opens another advantage forthe global disk simulations described in the next sections. Astrophysical disks are characterized by viscous forces due to turbulent motions which en-able the transport of circulating material towards the central source. This accretion process re-duces the mass of the disk but the physical nature of the viscosity is still under debate [23, 24].The timescale of disk dispersal by viscosity is much longer than the Keplerian orbital time whichis demonstrated in Sect 4.3 and Fig. 7.The usual viscosity description relies on the turbulence based approach for the kinematicviscosity ν according to Shakura & Sunyayev [25] ν = α c s H P , (15)where α is a free parameter, c s is the sound velocity and H P the vertical scale height (cf. 9). Sucha formulation is based on the idea that turbulent gas blobs are moving at speed similar to thesound speed and that the largest turbulent eddies are naturally limited by the scale height of thedisk. Numerical simulations of the magnetic rotational instability MRI [e.g. 26] can be used todetermine a more reliable estimate of α and typical values are assumed to be α ∼ .
01. Notethat this value for alpha represents the ideal MHD limit. In a realistic astrophysical scenarioalpha certainly can have a more complicated radial dependency [e.g. 23]. In this work we do notintent to investigate such realistic astrophysical applications but describe in detail a novel methodto study the evolution of PPDs. The dependency of the alpha parameter and the evolution of asimple viscous disk is discussed in more detail in Sect. 4.
Although the disk mass is small compared to the stellar mass, changes in the disk’s gravita-tional potential can play an important role when mass is redistributed within the disk. Althoughwe do not focus on detailed astrophysical scenarios where the disk potential becomes relevant6e.g. so-called FU Ori outbursts) we include in our numerical method the disk’s gravitationalpotential to obtain a more consistent modelling. Furthermore, we mention that the majority ofPPDs have masses around 1% to 5% of the stellar mass. Such disks are in any case gravitationallystable with respect to the Toomre parameter of Q > Q =
1, i.e. the disk mass reaches around 10% of the stellar mass[e.g. 27, 28, 29].The gravitational potential and the gravitational acceleration towards the central star of a thinaxisymmetric disk can be derived according to [30, page 73, Eq. 2-14b; Eq. 2-146] and further[31, 32] Ψ d ( r ) = − G (cid:90) r out r in Σ ( r (cid:48) ) k K ( k ) (cid:114) r (cid:48) r (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) (cid:13) dr (cid:48) (16)and ∂ Ψ d ∂ r = Gr (cid:90) r out r in Σ ( r (cid:48) ) (cid:34) K ( k ) − (cid:32) k − k (cid:33)(cid:32) r (cid:48) r − rr (cid:48) (cid:33) E ( k ) (cid:35) (cid:114) r (cid:48) r k (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) (cid:13) dr (cid:48) , (17)where r in and r out represent the inner and outer boundary of the disk and K ( k ) and E ( k ) is thecomplete elliptical integral of the first and second kind respectively. The variable k is defined as k = (cid:112) rr (cid:48) / ( r + r (cid:48) ) . For better readability we further substitute the terms 1 (cid:13) and 2 (cid:13) with M ij and N ij respectively. Since Σ is constant within each cell per timestep, the discrete form of Eq. 16and 17 can be written as Ψ d ( r i ) = n p − (cid:88) j = Σ ( r j ) (cid:90) r j + r j M i j dr j (18)and ∂ Ψ d ( r i ) ∂ r = n p − (cid:88) j = Σ ( r j ) (cid:90) r j + r j N i j dr j . (19)The indices i and j represent the i -th and j -th grid points and n p the total number of grid pointswithin the computational domain. The individual integrals (cid:82) r j + r j are numerically solved via atrapezoidal rule. The remaining elliptical integrals E ( k ) and K ( k ) are approximated with thearithmetic-geometric mean (AGM) [33]. A validation for the implemetation of Ψ d and a moredetailed discussion is presented in Sect. 4.5.
3. Numerical method
Since the method is based on the same principles as those presented in [9] we restrict thediscussion to specific issues typical for cylindrical geometry. The advection scheme is identicaland based on a second-order advection method according to [34]. Since we adopt a staggeredmesh, we have to distinguish scalar and vector quantities. The scalar quantities are defined insidea finite volume and vector quantities are located at the cell boundaries. The exact locations are7 igure 1: The overall picture of the discretization of an axial-symmetric configuration. The upper part of the figure showsthe projection of this logical discretization scheme into the physical domain. The boundaries of the physical domain arelocated at r and r npt − , respectively. The cells between the grid points r and r , or r npt − and r npt have no physicalvolume. These cells are necessary to provide the boundary conditions for the model variables and are therefore referredto as ‘ghost cells’. depicted in Fig. 1. Considering the shape of computational cells in cylindrical geometry, we haveto take care of the corresponding divergence terms, use a volume-weighted density, and use theappropriate definition of the grid velocity through adaptive Reynolds theorem. Hence, the scalarvolume V S is defined by V S , i ( t ) = π (cid:104) r + ( t ) − r i ( t ) (cid:105) . (20)For brevity, we omit from our notation the temporal dependence of the radial grid points r i in theremaining test. The flux term across a cell boundary is composed of two contributions resultingfrom the fluid flow with velocity u i and from a motion of the boundary itself, controlled by theadaptive grid, i.e. F S = π r i u i δ t − πδ (cid:16) r i (cid:17) , (21)where the symbol δ denotes the temporal di ff erence between of quantities between two timesteps.For vector-like variables, the position within a cell is given by r + = (cid:16) r + + r i (cid:17) , (22)which demands a definition of the corresponding volume V V for vector quantities through V V , i = π (cid:18) r + − r i − (cid:19) = π (cid:16) r + − r − (cid:17) (23)to ensure mass conservation.In order to allow a correct treatment of the divergence terms within the adopted staggeredcylindrical geometry, we calculate and conserve the flux by u i + r i + =
12 ( u i r i + u i + r i + ) , (24)when transforming for the cell boundaries to the cell center and obtain the advection over theupper cell boundary at r i + F v = π r i u i + r i + u i + ) δ t − π δ (cid:16) r i + r i + (cid:17) . (25)8or a coordinate system moving with the velocity u grid , we introduce the adaptive grid trans-port theorem of [35] which is a generalization of the well-known Reynolds transport theorem inhydrodynamics, given by ddt (cid:32)(cid:90) V f dV (cid:33) = (cid:90) V (cid:34) ∂ f ∂ t + ∇· (cid:16) f u grid (cid:17)(cid:35) dV . (26)The motion of the grid points has to be defined in such a way that a constant function f re-mains constant if the grid moves. Assuming only radial motions, we end up with an appropriatedefintion of the grid velocity through2 π ru r grid δ t = πδ (cid:16) r (cid:17) . (27)Note that a Lagrangeian motion, i.e. u grid = u , is then consistent with no fluxes across the cellboundaries. The discretization of the hydrodynamic equations has been discussed in detail by [9] and thuswe refrain from defining the advected density ˜ Σ which is obtained using a second order advectionmethod [see 34]. Accordingly, the derivation of the relative velocity u rel between the grid cellsand the gas is not part of this work. Furthermore, in order to present the hydrodynamic eqationsin a more compact form, we are applying an index free notation and we denote the temporaldi ff erences with δ and the spatial di ff erences with ∆ .Applying Eq. 26 to Eq. 4 we obtain ∂∂ t (cid:90) V ( t ) Σ dV + (cid:90) ∂ V ( t ) Σ u r dS = . (28)Utilizing Eq. 20 the discrete form of the equation of contiuity is given by δ [ ρ V s ] + (cid:88) (cid:101) Σ ∆
Vol s = , (29)where (cid:101) Σ denotes the surface density to be advected with the advective flux during the time step δ t over the cell boundary and therefore (cid:80) (cid:101) Σ ∆
Vol s represents the advection term.Integrating Eq. 5 over the volume and since u grid only has a radial component, the discretisa-tion of the radial equation of motion takes the form δ [ Σ u r V v ] + (cid:88) (cid:103) Σ u r ∆ Vol v − Σ u ϕ u ϕ r V v δ t + π r √ π ∆ (cid:104) H p P gas (cid:105) δ t − ∂ r Ψ ( r ) tot Σ ∆
Vol v δ t − π r ∆ (cid:26) r µ Qr (cid:20) ∆ u r ∆ r − u r r (cid:21)(cid:27) δ t = , (30)where ∂ r Ψ ( r ) tot tags the total radial gravitational acceleration (see Sect. 4.5) and µ Qr is the ra-dial viscosity coe ffi cient. Quantities denoted with F are the averaged values transformed to thecorresponding grid point in the staggered grid. The terms from left to right can be identified asthe temporal di ff erence, advection, the centrifugal force, the gas pressure gradient, gravitationalacceleration and viscosity. To simplify the notation, we have abbreviated the radial di ff erencebetween r i and r i + by the spatial operator ∆ . 9s a consequence of our 1D-method the velocity component u ϕ is projected on the purelyradial grid. Hence, Eq. 7 is discretized as a scalar-like quantity δ [ r Σ u ϕ V s ] + (cid:88) (cid:103) r Σ u ϕ ∆ Vol s − π ∆ (cid:26) r µ Q ϕ (cid:20) ∆ u ϕ ∆ r − u ϕ r (cid:21)(cid:27) δ t = , (31)where µ Q ϕ is the angular scaling factor of the viscosity. In principle, µ Qr and µ Q ϕ can havedi ff erent values, but presenting and validating the numerical method, we set them to identicalvalues allowing a direct comparison with analytical viscous disk solutions (see Sect. 4.4).Equation 14 requires additional considerations. Since the internal energy balance of the diskis dominated by the irradiation term and the radiative cooling term, both implemented directlyinto the equation of internal energy, the temperature structure does not follow the radiation tem-perature (J), which in itself is computed correctly. The obvious solution for this problem is toimplement the radiative energy flux also directly into the internal energy equation versus intothe radiation energy equation, i.e. assuming J = S (since the disk is optacally thick in radialdirection). Further assuming an almost radial isotropic radiation field, i.e f edd = /
3, we candrop the 3 K − J term and inserting Eq. 3 leads to √ π H p π ∆ ( rH ) δ t + κ P Σ ( J − S ) V s δ t = . (32)Moreover, the time-dependent terms in the radiation flux equation are usually not importantcompared to the pure radiation di ff usion terms. Hence, it makes sense to skip the time-dependentterms and the 3 K rr − J and H ∇ · u terms, which gives23 r σ ∆ ( T ) + κ R ρ H V v = , (33)leading to H = − r σ ∆ ( T ) κ R Σ V v . (34)Utilizing the last equation, the discrete form of the equation of internal energy simplifies to δ [ Σ e V s ] + (cid:88) (cid:103) ( Σ e ) ∆ Vol s + P (cid:108) π ∆ ( ru r ) δ t + π √ π H p σ ∆ − r ∆ ( T ) κ R ρ V v δ t + ∆ E rad − µ Q (cid:26)(cid:20) ∆ u r ∆ r − u r r (cid:21) + (cid:20) ∆ u ϕ ∆ r − u ϕ r (cid:21) (cid:27) V s δ t = , (35)where the last term represents viscous energy generation due to friction. The updown arrow (cid:108) denotes z-integrated values. The ∆ E rad stands for the radiation heating / cooling term in verticaldirection ∆ E rad (see Eq. 14), which is discussed in detail in Sect. 3.2. Radiation transport is a crucial component when simulating viscous disks because radiationheating and cooling leads to changes in the thermal disk profile and consequently to changes in10he mass flux through the disk. Since astrophysical disks are predominantly optically thick in theradial direction, the radial radiation transport is only defined locally and thus has a minor impacton the long term evolution of disks. However, the radiative energy from the central star (andconsequently the vertical radiation transport), becomes the dominant heating / cooling process.Our method adopts an axisymmetric geometry with vertically integrated values and thereforethe vertical radiation transport needs a deeper discussion. The star’s radiation hits the disk ata shallow angle that depends on the inclination of the disk and the stellar radius. The verticallocation at which the star’s radiation penetrates the disk is defined as H phot and represents thevertical distance from the disk midplane to the photosphere of the gas disk. Therefore, thevertical scale height of the disk H p and H phot are constrained via a constant value α phot H phot = α phot H p . (36)The area which is exposed to solar radiation is given by 2 π r ∆ ( H phot / r ), where ∆ ( H phot / r ) repre-sents the inclination of the surface area to the line-of-sight of the central star. Thus, the irradiationheating rate from the star is ∂ t E irr = − a ) L (cid:63) π r π r ∆ (cid:32) H phot r (cid:33) , (37)where a repersents the disk’s albedo and L (cid:63) the stellar luminosity. The factor 2 appears whentaking into account that the disk is radiated from both the top and bottom. Utilizing the relationin Eq. 36, combining all constants in a factor f irr = (1 − a ) α phot and, in order to prevent negativevalues for E irr , the lower limit of E irr per timestep δ t gives E irr = L (cid:63) f irr max (cid:34) ∆ (cid:32) H p r (cid:33) , (cid:35) . (38)Furthermore, the physical size of the central star has to be taken into account as the stellar ra-diation is assumed to originate from a certain distance H (cid:63) from the midplane of the disk. Thedistance H (cid:63) can be computed from a half circle with stellar radius R (cid:63) using H (cid:63) = π R (cid:63) . (39)Implementing Eq. 39 into Eq. 38 leads to the radiation enrgy from a star on the disk surface at acertain orbital radius, given by E irr = L (cid:63) f irr max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) . (40)To describe the radiative cooling of the disk, we assume that the radiation of a gid-cell volume V s can be specified by the black body radiation and hence the cooling energy E cool per timestepis E cool = V s σ T δ t , (41)where σ is the Stefan-Boltzmann constant and T surf the local surface temperature of the disk.Taking into account that the heating energy from the star and the cooling energy from the diskhave to be transported in the z -direction, the radiative energy balance ∆ E rad (see. Eq. 14) is ∆ E rad = π H z V s δ t . (42)The final step is to calculate the radiation flux in the vertical direction H z which is illustrated inSect. 3.2.1. 11 .2.1. Radiation flux in z-direction In order to compute the radiation flux in the vertical direction, we utilize the stationary radi-ation transport equation n · ∇ I ν = η ν − χ ν I ν . (43)Integration over all directions and frequencies gives the radiation energy equation (zero th mo-ment) 14 π (cid:90) d Ω : ∇ · H ν = η ν − χ ν J ν (cid:90) d ν : ∇ · H = η − χ P J η = χ S = χ P ( S − J ) ∂ H z ∂ z = − κ P ρ ( J − S ) (44)while integration with n yields the first moment equation of radiation14 π (cid:90) n d Ω : ∇ · K ν = − χ ν n I ν = − χ ν H ν (cid:90) d ν : ∇ · K = − χ R H z ∂ K zz ∂ z = − χ R H z = − κ R ρ H z . (45)In proceeding further we list and apply our idealizing assumptions: Eddington factor: K = f edd J , or in z-components K zz = f edd J = J . Local thermal equilibrium (LTE): J = S = σπ T Integration of the radiation energy equation thus gives ∂ H z ∂ z = − κ P ρ ( J − S ) = −→ H = const. (46)Integration of the radiation flux equation yields ∂ K zz ∂ z = ∂ J ∂ z = − κ R ρ H z = σπ ∂ T ∂ z = − κ R ρ H z = σπ T − T H p = − κ R ρ H z (47)The position of the disk’s photosphere H phot would be the proper z-localization for the surfacetemperature T surf , however, the position of midplane temperature T is less clear as the hydro-dynamic discretization still assumes an isothermal structure. Therefore, using H p as the spatialseparation of T surf and T in the last equation above seems to be a reasonable compromise. Herewe also specified the relevant density to be the midplane density ρ , since we only have one12pacity which corresponds to the midplane, and the optical depth would likely be dominated bythe dense layers close to the midplane.Rewriting the last line, we obtain an expression for H z , given by H z = − σ πκ R ρ H p (cid:16) T − T (cid:17) . (48)Identifying the optical depth in z direction as τ = κ R ρ H p (i.e. τ = κ R ρ (cid:108) √ π ), we obtain H z = − σ πτ (cid:16) T − T (cid:17) (49)The contribution to the energy balance of the scalar cell i in the time interval δ t , computed fromthe H z flux above, is given by ∆ E = π c ∆ J = π H z V s δ t = − π σ πτ (cid:16) T − T (cid:17) V s δ t , (50)where a factor 2 enters to allow for radiation over both the upper and lower surfaces, and thefactor 4 π relates the first moment of radiation H to the energy flux F rad . In order to treat the fluxdirection correctly, we assume the radiation from the disk to the surrounding space as positivefor H z , a positive ∆ E and therefore the cell loses energy (i.e. δ [ ρ e ] + ∆ E + . . . =
0, or δ [ J ] + c π ∆ E + . . . =
0, respectively).
To calculate the vertical energy flux and thus the contribution to the internal energy equation,we have to compute an estimate for the surface temperature, T surf , of the disk. In order to do so,we assume local balance between irradiation and radiative cooling on the disk surface using thedescription from Sect. 3.2. Additionally, we consider the vertical energy transport from the diskinterior to the disk surface which corresponds to the irradiation / cooling term ∆ E rad in the internalenergy equation (Eq. 14) and in the following, we denote this quantity F vert . Hence, the energybalance for the disk surface is F vert + E irr + E amb − E cool = , (51)where E amb is the irradiation from the ambient radiation field for which we assume, as for theradiative cooling E cool , black body radiation. The contributions then is F vert = σ τ (cid:16) T − T (cid:17) V s δ t , (52) E irr = L (cid:63) f irr max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) δ t , (53) E amb = V s σ T δ t , (54) E cool = V s σ T δ t . (55)Summing these terms up, skipping δ t , and dividing by 2 σ V s gives43 τ (cid:16) T − T (cid:17) + L (cid:63) f irr σ V s max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) + T − T = τ T + L (cid:63) f irr σ V s max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) + T = T (cid:32) + τ (cid:33) . This leads to an expression for T , given by T = τ + τ T + + τ L (cid:63) f irr σ V s max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) + + τ T , (56)Since we have adopted the Eddington-approximation for calculating the vertical temperaturestructure, the optical thin limit usually overestimates of the radiative losses. In order to correctthis feature a flux-limited di ff usion approximation can be utilized [e.g. 16]. The radiative energylosses can be computed for both, the optically thin and optically thick case, by rewriting theterms with τ as τ + τ = + τ , + τ = τ + τ . Insertion of the result for T into the equation for the vertical radiation flux gives ∆ E rad = σ τ (cid:16) T − T (cid:17) V s δ t . (57)This leads to ∆ E rad = σ τ T − + τ T − τ + τ L (cid:63) f irr σ V s max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) − τ + τ T (cid:19) V s δ t . Using 1 − + τ = τ + τ , and rewriting the leading factor, gives ∆ E rad = σ V s τ ( τ + τ T − τ + τ L (cid:63) f irr σ V s max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) − τ + τ T (cid:19) δ t . (58)By canceling τ and introducing a modified optical depth τ (cid:48) with the definition τ (cid:48) = τ (cid:32) = κ R ρ H p (cid:33) (59)14he result can be written in the form of ∆ E rad = σ V s + τ (cid:48) (cid:16) T − T (cid:17) δ t − + τ (cid:48) L (cid:63) f irr max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35) δ t . (60)Note that Eq. 60 is valid if the disk is optically thick. As the optically thick part of theirradiation is always dominant, the optically thin part can be neglected (if we reach the opticallythin limit, we simply have an isothermal structure). Nevertheless, we have to account for theweakening of the irradiation energy when approaching the optically thin limit. This can beachieved by adding a weighted factor to Eq. 60. For the optically thick limit, we can define τ n + τ n → (cid:40) τ (cid:28) , τ (cid:29) . In the optically thin case, we have τ → + τ n → , + τ (cid:48) → . Therefore, to account for an optical thin disk, the exponent must satisfy n ≥
2. Utilizing τ (cid:48) = / τ leads to the final form of the heating / cooling input into the disk ∆ E rad = τ + τ δ t (cid:18) σ V s + τ (cid:48) (cid:16) T − T (cid:17) − L (cid:63) f irr max (cid:34) ∆ (cid:32) H p − H (cid:63) r (cid:33) , (cid:35)(cid:19) . (61) To provide a su ffi cient radial resolution, we use an adaptive grid which is augmented to thephysical equations and redistributes the gridpoints at every timestep. Following the strategy of[14], where all details of the adaptive grid are presented, the basic parameters are controllingthe spatial variation and the temporal smoothing of the grid motion. In particular, the temporalgrid-scale τ g ( r ) is typically set to the global di ff usive timescale ¯ τ ν (see Sect. 4.4) and we use τ g = ¯ τ ν . Since the spatial variation of the physical variables is expected to be rather smooth, onlythe surface density Σ enters the desired resolution.In Fig. 2 we illustrate the redistribution of gridpoints towards an artificial density feature.Starting from an initial model with a homogeneous density profile on which we impose a dis-turbance, Σ , at a certain radial distance R (top panel). After initiating the simulation and onlysolving the grid equation, one can see that the initially logarithmic equidistant gridpoints movetowards the density peak leading to a higher gridpoint concentration around the disturbance (bot-tom panel). This behaviour demonstrates the dynamic improvement of accuracy at regions withphysical features, e.g. wave-fronts. 15 igure 2: Propagation of gridpoints towards a density feature imposed on a homogeneous background. The upper panelshows the surface density Σ in values of a reference density Σ against the radius R in values of a reference radius R . Thedensity feature has a maximum of Σ at the location of R . Starting from an initially logarithmic equidistant gridpointdistribution, the gridpoints accumulate around the density feature to increase the local resolution (lower panel). .4. Numerical solution procedure Summing up all radial indices of the discrete equations (see last sections), we find that thedi ff erence equations contain variables at 5 gridpoints, i.e. the discrete equations connect [ i − , i − , i , i + , i + m × m -matrix, where m is thenumber of unknowns per gridpoint. In the simplest case, we have m = Σ , u r , u ϕ , e , H , Ψ d and r and for 500 radial gridpoints the whole system of algebraic equationscontains 7 × = δ x of δ x > − . More details onthese numerical solution issues can be found in [9]. Boundary conditions according to Sect.4.1are implemented through ghost cells to keep the overall matrix pattern simple.
4. Numerical examples, accuracy and comparison to analytical solutions
To assess the validity of our model when it is applied to astrophysical viscous disks and tomeet their specific requirements mentioned in the previous sections, we tested our code againstseveral theoretical models. In Sect. 1, we emphasized the importance of the inner boundarycondition and hence in Sect. 4.1, we show that leaving out the very inner regions of the diskleads to an entirely di ff erent disk structure. Thus results from models with large disk-to-star gapshave to be treated carefully. In general, defining an exact outer boundary for the disk is di ffi cultas there is a pressure gradient towards the surrounding ISM. Our model has no restrictions ofwhere to set the outer boundary and hence we can cover several orders of magnitude of diskradius. As viscous disks are defined via a thermal and gravitational profile, permuting the disk ineither the thermal or the gravitational variables should lead to a wave propagation in either of thevariable spaces (temperature-radius and density-radius space). In Sect. 4.2 and 4.3, we show thatour model reproduces the expected behaviour which validates the model setup. Viscosity playsan essential role in astrophysical disks as it determines the mass accretion onto the star. Sincewe have implemented a viscosity model (see Sect. 2), we test our model against an analyticalsolution (cf. Subsect. 4.4). As mentioned in Sect. 1, we emphasize that the boundary conditions play an important rolein the long term evolution of astrophysical, viscous disks and hence we test our model accord-ingly. We have stated that changing the inner boundary of the disk consequently leads to entirelydi ff erent disk structures whereas varying the outer boundary, the disk structure remains similar.To verify this, we start from a reference disk with a fixed inner and outer radius ( R in = . R out = / R ref (min values) and 10 × R ref (maximum values). Note that the exact radii values are notimportant for the qualitative simulation result but should be mentioned for completeness.To make disks comparable, the mass flux of the di ff erent models must be equal as the fluxis entirely determined by viscosity. In Fig. 3, we plot the disk’s surface density profile againstthe radius in dimensionless variables. Panel a) shows simulations in which the inner boundaryis varied with a constant outer boundary and panel b) shows simulations in which the outer17 nner radius outer radius reference disk 0.05 AU 25 AU R in;min R in;max R out;min R out;max Table 1: List of boundary conditions tested with our model. boundary is varied with a constant inner boundary. We assume pressureless outflow, ∂ P /∂ r = P = / ff erent surfacedensity profile. As the surface density is tightly coupled to the viscosity, such a variation mustresult in a di ff erent evolution track for the disk (e.g FU Ori bursts, Ex Lup objects, etc.).If varying the outer boundary condition, represented by Panel b), we see that the disk structurestays una ff ected. This outcome also supports the assumption that the density in the outer regionsof the disk is of the order of the ISM density ( ρ ISM = − gcm − ) [36] and therefore viscousdrag is negligible. We emphasize that our model can cover large spatial dimensions, meaningthat we do not have any limits in where to set the outer boundary. For numerical accuracy andsimulation time optimisation, we set the outer boundary to that of the outer base radius accordingto the standard disk model [18]. When increasing the host stars luminosity, the thermal pressure in the disk will increase dueto radiation. The higher pressure breaks hydrostatic equilibrium in the vertical disk structureand thus the disk will expand in the vertical direction. Changes in the disk dynamics can onlyappear on the viscous timescale of the disk which increases with τ ν = r /ν further from thestar. Additionally, the surface temperature of the disk decreases with increasing distance fromthe centre and as a result, the innermost disk regions are expected to adjust faster to the star’sthermal input compared to regions further away. On the contrary, due to the inwards-orientedaccretion flow, regions closer to the star are cooled by cold gas from the outer regions of thedisk. Hence, the local τ ν decreases when getting closer to the centre which leads to a delay inthe thermal adjustment of these disk regions. Outer disk regions are less a ff ected and thereforeadjust faster to the thermal input of the star and as a consequence, the thermal profile of the diskis expected to approach the final state from both inwards and outwards, with less e ff ective zonesin the middle of the disk.To test our model accordingly, we started from a stationary viscous disk around a Sun-likestar. We then spontaneously increased the stars luminosity by an arbitrary factor ∆ L . We foundthat regardless of the exact value of ∆ L , the qualitative outcome of our simulation remains similarand therefore all our simulations have been carried out with ∆ L =
10 for illustration purposes.In Fig. 4, we show the thermal evolution of the disk from the initial model ( τ , solid) to the final18 )1.00.50.0 d i s k / a)a)a)a)a)a) R in , min reference diskR in , max )b)b)b)b)b)b) R out , min reference diskR out , max Figure 3: Surface density Σ disk in values of an arbitrary reference surface density Σ against disk radius R in values ofa reference radius R . Panel a) shows the structure of a stationary viscous disk when only varying the inner boundaryand panel b) shows the structure of a disk when only varying the outer boundary and leave the inner boundary constant.In both cases, we modify the inner / outer radius by 1 / R ref (dashed) and 10 R ref (dashed-dot), where R ref represents theinner / outer radius of the reference disk. In both cases, all models have equal mass flux. One can see that varying theinner boundary leads to entirely di ff erent disk structures whereas when changing the outer disk radius, the disk remainsunaltered. stationary solution ( τ , dashed). The plot shows the gas temperature as a function of the diskradius. As expected, the inner region responds immediately to the updated thermal input fromthe host star ( τ , dashed-dotted) as the very inner disk-cells contain less material according tothe pressureless outflow condition at the inner boundary (compare Sect. 4.1). Consequently, thelocal viscosity is less e ff ective and additionally the thermal input of the star is near the maximum.Therefore, the disk gas temperature can adjust to the updated thermal input on a short timescale.According to the radially increasing τ ν , the temperature only gradually increases with distancefrom the disk centre ( τ , dash-dot dotted) and due to mass accretion from the outer still colddisk regions, the local temperature adjustment is slowed down significantly at around 30 R .Simultaneously, the very outer grid cell is not a ff ected by accretion of cold gas and thus thetemperature can freely adjust to the radiation from the star ( τ , dotted). The thermal profileof the disk then adapts gradually from the outer and the inner regions until the final stationarysolution is reached. We emphasize that the timescales on which the thermal adjustments occur( τ to τ ) are about three times the order of the mean Keplerian rotation period and hence utilizingan implicit scheme is necessary in order to get reasonably large timesteps to investigate the fulllifespan of astrophysical viscous disks. As mass is transported inwards due to the gravitational pull of the host star and viscousfriction within the disk (compare Sect. 2), mass flow through the inner boundary of the disk isgenerated. Part of this mass is accreted onto the central star and therefore increasing the star’smass results in a change of the gravitational star-disk potential. A change in the gravitationalpotential directly results in an alteration of the Keplerian velocities in the disk. Areas close to19 .0 0.5 1.0 1.5 2.0 2.5 3.0log(R / R )2.01.51.00.50.0 l o g ( T g a s / T ) Figure 4: Temperature evolution of a viscous disk when spontaneously increasing the host star’s luminosity. We showthe gas temperature against the radius in dimensionless values. The initial model is illustrated by τ (solid) and the finalmodel after changing the luminosity of the central star is shown by τ (dashed). One can see that the innermost grid cellsrespond quickly to the updated thermal input ( τ , dashed-dotted). Since the viscous timescale τ ν decreases with distancefrom the centre and because of accretion flow from the outer cold regions of the disk, the temperature gradually adaptsto τ from the outer and inner parts of the disk ( τ , τ ). l o g ( d i s k / ) a)a)a)a) )420 l o g ( d i s k / ) b)b)b)b) Figure 5: Propagation of the expected density wave in the density-radius space. All values are in dimensionless values.The black solid lines represent the initial (Panel a) and the final models (Panel b). Panel a shows the onset and propagationof the wave at several di ff erent times ( τ to τ , dashed). Once the wave reaches the outer boundary, the settling begins( τ to τ , dot-dashed) which is shown in Panel b. The solution converges towards the final, stationary model ( τ ) byslowly balancing the mass-flux at the inner and outer boundaries of the disk. the star will be a ff ected much quicker than regions further from the centre and thus mass willbe redistributed in the disk. This redistribution of mass is expected to be visible as an outwardspropagating ‘density wave’. To test that our model reproduces this behaviour, we start from astationary viscous disk and multiply the stellar mass by an arbitrary factor ∆ m . All other stellarparameters remain constant. Note that the exact value of ∆ m only alters the intensity of thewave and therefore all our calculations have been carried out with ∆ m = .
5. The change in thestar’s mass does not just increase the gravitational potential, but also alters the e ff ective viscosity(cf. Sect. 2.2). We emphasize that in our model a change in stellar mass instantly changesthe Keplerian velocity and gravitational pull in each grid cell. Thus, the disk experiences anincreased mass flow in each cell which leads to mass-flux instabilities if the time scale of themass accretion onto the star ∆ t exceeds the local time scale of the viscous disk τ ν = r ν − .Therefore, the mass accretion onto the star must be distributed over a timescale that is shorterthan the viscous timescale of at a reference radius ∆ t < τ ν ( r ref ) (cf. Fig. 6). The shortest viscoustimescale appears to be at the position with the highest e ff ective viscosity. We found that τ ν atthe inner boundary is always close to τ ν ;min and thus we set our reference radius to r ref = r in .Figure 5 shows the temporal evolution of the surface density after gradually increasing the21 l o g ( d i s k / ) t = 0.1 t visc t = 0.1 t visc t = 0.1 t visc t = 0.1 t visc ) t = 1.0 t visc t = 1.0 t visc t = 1.0 t visc t = 1.0 t visc t = 10.0 t visc t = 10.0 t visc t = 10.0 t visc t = 10.0 t visc Figure 6: Simulation of the protoplanetary disk subjected to a density wave for di ff erent values ∆ t = . , . , . t visc as the timescale for the gradual increase of the stellar mass. The plots show five arbitrary times starting from the initialstationary model ( τ ; solid) across the transmission of the density wave ( τ τ ; dash-dotted, dotted) until stationarity isattained again ( τ ; dashed). The graphs in the panels only show about 30% of the outer disk as here the density wavesappear with the most intensity. stellar mass of a stationary viscous disk ( τ ; solid, Panel a). As the central mass increases,mass accumulates at the inner boundary resulting in a higher surface density ( τ ; dashed, Panela). At disk regions further away from the centre, viscosity is not e ffi cient enough to transportmass accordingly and a pile-up of mass occurs. This ‘mass-bumps’ occur because density wavespropagate towards the outer boundary ( τ and τ ; dashed-dotted and dotted, Panel a). Once thedensity waves reach the outer boundary ( τ ; dashed, Panel b), the mass flow within the disksettles but leaves the disk with an imbalanced mass accretion over the inner and outer boundaries(mass accretion is higher over the inner boundary due to more e ffi cient viscosity). Hence, thedisk loses mass ( τ and τ ; dot-dashed and dotted, Panel b) until a stationary solution is reached.As expected, a higher central mass causes a larger gravitational pull and a higher rotational speedof the disk and thus, compared to the initial model, slightly less (in this case ∼ τ ; solid, Panel b). These results show that our modelcan reproduce the behaviour expected in a scenario where the mass of the host star is increasedand is therefore capable of producing comprehensive results when studying viscous disks.We have mentioned that the mass accretion onto the central star has to be spread over atimescale that is shorter than the viscous timescale at a reference radius r ref . Since r ref withthe shortest viscous timescale is di ffi cult to determine but the timescale on that mass is accretedonto the star might influence the result, we have studied the impact on the results when usingdi ff erent mass acrretion timescales. Figure 6 shows the behaviour of a stationary viscous diskwhen the central mass is increased by a factor of ∆ m = . ff erent timescales of ∆ t = . τ ν, in (left panel), ∆ t = τ ν, in (middlepanel) and ∆ t = τ ν, in (right panel), where τ ν, in is the viscous time at the inner boundary. Onecan see that if mass is applied to the central star too slowly, the disk can adjust to the changinggravitational potential and the density waves get absorbed. This leads to the conclusion that ifone wants to study mass infall correctly (e.g. when investigating FU-Ori bursts, infall of planets,etc.), it is even more important to compute the very inner disk regions consistently in order notto lose any physical feedback onto the disk. Moreover, density adjustments occur on timescaleswhich are up to seven orders of magnitude larger than the mean Keplerian rotation period. Thus,for studying the long term evolution of astrophysical, viscous disks, only an implicit scheme can22rovide reasonably large timesteps. From the previous section, we have seen that viscosity plays a central role in the evolutionof disks and thus the reliability of our viscosity model is an important factor. Accretion diskscan be considered as axisymmetric thin disks [18]. The disk can be separated into annuli withdi ff erent average densities (the density distribution of one single annulus can be considered tobe constant). According to v r = r Ω ( r , t ), the radial velocity is di ff erent for every disk radiusand as the disk gas is a viscous fluid, friction between the annuli causes an energy loss from theinnermost to the outermost disk regions. As a result, the inwards oriented mass-flux ˙ m occurs. If˙ m is constant over all radial rings, the disk can be called a viscous stationary disk. If an externalmass is applied anywhere onto the disk, the mass has to be dispersed over the disk (inwardsand outwards) until the disk reaches a stationary solution again. To test our model, we consideran axisymmetric disk with a given surface density Σ . According to [37] or [18], we want toreproduce the time-dependent analytical solutions for the evolution of the surface density Σ ( r , t )of a geometrically thin disk under the action of internal angular momentum transport ∂ Σ ∂ t = r ∂∂ r (cid:34) √ r ∂∂ r (cid:16) ν Σ √ r (cid:17)(cid:35) . (62)This equation represents a di ff usive partial di ff erential equation for the surface density, which canbe derived by applying a variable substitution, assuming constant viscosity ν = const . . Defining X = √ r and f = Σ X , (63)Equation 62 takes the form of a typical di ff usion equation ∂ f ∂ t = D ∂ f ∂ X , (64)where D is the di ff usion coe ffi cient given by D = ν X . (65)Although a constant viscosity is not necessarily realistic for a protoplanetary disk, a Green’sfunction solution and the qualitative illustration of the behaviour of Eq. 62 is possible. Initialconditions at t = m at radius r Σ ( r , t = = m π r δ ( r − r ) , (66)where δ ( r − r ) is the Dirac delta function. Boundary conditions that enforce zero-torque at r = r → ∞ yield (see [37]) Σ ( x , τ ) = m π r τ x − exp (cid:34) − (1 + x ) τ (cid:35) I / (cid:32) x τ (cid:33) (67)for the time-dependent solution for Eq. 62, where x and τ represent unitless variables x = rr τ = ν r − t , I / is the modified Bessel function of the first kind.Note that contrary to the pure analytical solution of [18], our model solves a coupled systemof di ff erential equations and thus we have to adapt our numerical model to fit the analytical case.To do so, we start with a disk with an outer radius of r out =
10 AU surrounding a central starwith M (cid:63) = (cid:12) . Our initial model has a density profile similar to Eq. 66 and since a Diracdelta function is numerically not possible, we adopt a Gaussian distribution and normalized thesurface density Σ with an initially constant surface density Σ Σ ( r ) = Σ r √ π r | a | exp (cid:34) − ( r − r ) a (cid:35) + Σ , (69)where Σ represents a background density, r = ( r out − r in ) is the initial position of the grid celland a is comparable to the standard deviation. In our simulation, we solve the the equation ofcontinuity, the equation of motion and the equation of energy for a constant kinematic viscosity ν . Since our initial model needs time until it has adapted to the global viscosity and the densityprofile at t = ∆ t to align with theanalytical model. We found that if we substitute the averaged viscous time ¯ t visc for τ in Eq. 67 sothat τ visc = ν r − ¯ t visc , we can compute the ∆ t that best fits numerical model best to the analyticalmodel. The constant time correction factor is given by ∆ t = ¯ t ν t = const . (70)The remaining value m in Eq. 67 is given by m = (cid:88) π r Σ ( r ) ∆ r . (71)In Figure 7 we show the viscous evolution of a disk fragment with a background surfacedensity profile Σ on which we imposed a disturbance Σ at a reference radius r , compared tothe analytical solution at di ff erent times τ n . The numerical results are represented by solid linesand the analytical results by dashed lines. One can see that with time the density feature isdissolving to the left and the right of r . Simultaneously, the maximum drifts inwards causedby the inwards-oriented mass-flux ˙ m . Our model shows the same behaviour as predicted by theanalytical solution and thus verifies the correctness of our viscosity model. In Section 4.3 we show that altering the gravitational potential of the star can lead to seri-ous changes in the disk structure. Taking into account that our main goal is to investigate thebehaviour of viscous disks, the transport of mass through the disk will a ff ect the local viscositydue to changes in the disk’s gravitational potential. Even though the gravitational potential ofthe disk is negligible compared to the star’s potential, even small local changes might a ff ect theglobal outcome when investigating the long term evolution of viscous disks. To test the correctimplementation and behaviour of Eq. 16 and Eq. 17, we test whether the total gravitational po-tential is conserved, meaning (cid:82) ∞ Ψ tot ( t ) dr = const. ∀ t . As our computational domain is radiallyrestricted from a few stellar-radii to several hundred AU, there is a di ff erence between the fullyanalytical solution and the numerical result. However, the conservation-error in our simulationsis less than 0 .
3% which is acceptable in the scope of astrophysical disk simulations. Further-more, we emphasize that the disk’s potential has to adapt to a change in the mass distribution of24 .4 0.6 0.8 1.0 1.2 1.4 1.6r / r / analyticalnumericalanalyticalnumericalanalyticalnumericalanalyticalnumericalanalyticalnumerical Figure 7: Comparison of the evolution of a surface density feature in a viscous disk with the analytical model presentedby [18]. Numerical (solid lines) and analytical models (dashed lines) are plotted at di ff erent times ( τ to τ ) in the surfacedensity-radius space in unitless variables. We see that our viscosity model well reproduces the analytical solution andthus supports the correct implementation into our code. the disk. Hence, we tested this behaviour by inserting a mass M test , which is distributed over anannulus within the axisymmetric description at an arbitrary position R M into the disk. Figure 8(lower panel) shows the evolution of the disk potential at di ff erent times τ i . The upper panelrepresents the simultaneous adaptation of the surface density. Focusing on the lower panel, wesee the formation of a potential cavity ( τ , solid) after inserting M test into the initial model ( τ ,solid-triangle). Due to viscous forces, M test is dissolved radially ( τ , dotted), which results in adepletion of the cavity. The slightly higher disk-to-star mass ratio causes a change of accretionrate at the outer and inner boundary ( ˙ m in and ˙ m out ). As this change in ˙ m is not equal for the innerand outer boundary (the inner region is a ff ected more due to the gravitational pull from the hoststar), more mass is accreted over the inner boundary than the outer boundary, resulting in a massloss of the disk. Once the disk mass is equivalent to the mass of the initial model, ˙ m is againconstant over the entire disk, representing a stationary solution ( τ ; dashed) which is equivalentto the initial model. We emphasize that there is a maximum disk-to-star mass ratio if when ex-ceeded, no stationary solution is possible [38]. A detailed investigation of this result should beconsidered in further disk evolution studies.
5. Conclusion
In this paper, we present an implicit numerical method to solve axially symmetric time-dependent equations of radiation hydrodynamics. We assume hydrostatic equilibrium perpen-dicular to the equatorial plane, which can be described as a 1 +
1D approach. Although the nu-merical method is currently restricted to 1-dimensional problems, i.e. all variables X = X ( r , t )depend on radius and time, such implicit computations allow an accurate description of axialsymmetric configurations and their evolution over time scales much larger than several orbitalperiods. Global phenomena like mass accretion and angular momentum transport can be studied25 igure 8: Reaction of the disk by perturbing the initial model with a mass M test at an arbitrary radius r ref . The upperpanel shows the evolution of the disk’s relative surface density, ∆Σ ( r ), and the lower panel the adaptation of the relativegravitational potential, ∆Ψ ( r ), where ∆Σ ( r ) = Σ ( r ) − Σ ( r ref ) and ∆Ψ ( r ) = Ψ ( r ) − Ψ ( r ref ) respectively. The dash-dottedvertical line in both panels represents the position of the mass insertion r ref . On the horizontal axis, we have the radiusin values of the reference radius R = R out , where R out represents the outer disk radius. Focusing on the lower panel,we see the formation of a potential cavity ( τ ; solid) after applying the M test to the initial model ( τ ; solid-triangle).Due to viscous forces, mass will be transported inwards (to the left of r ref ) and outwards (to the right of r ref ), leadingto a widening of the cavity ( τ ; dotted). The drift of the cavity minimum to the left is caused by the inwards-orientedmass-flux. The mass distribution continues until the perturbation M test is dispersed over the entire disk. This resultsin di ff erent mass accretion rates ˙ m over the inner and outer boundaries and consequently, due to the inwards-orientedmass-flux, to a mass loss over the inner boundary. The mass loss continues until the disk mass is equal to that of theinitial model ( τ ; dashed).
26n long evolutionary timescales. In order to test the validity of our method, we have tested ourmodel for various configurations.Since the paper is focused on the numerical method, the presented examples are calculatedwith the simplest material functions, e.g. an ideal equation of state, constant opacity, and con-stant viscosity. In our astrophysical computations, additional features are introduced by the morerealistic material functions, i.e. opacity changes due to ionization or dissociation or variations inthe viscosity caused by the onset of magnetic instabilities [26, e.g.]. From spherical computa-tions of pulsating stars [39, e.g], we can conclude that such features can be nicely resolved bythe adaptive grid and do not lead to additional complications as long as the derivatives on thedependent variables remain smooth.We emphasize that radiation transport is an essential process for the thermal structure of thedisk and thus has to be considered in astrophysical disk simulations. Changes in the thermalprofile lead to changes in viscosity and consequently to a di ff erent mass flux through the disk.Since the radial radiation flux is only defined locally, because the disk is optically thick in radialdirection, the vertical radiation transport becomes the dominant heating / cooling term. Since weadopt cylindrical geometry, a proper description for the radiation transport in vertical direction(see Sect. 3.2) has to be defined.Furthermore, our results from Sect. 2.1 show that a consistent description of the entire diskis necessary in order to produce physically realistic results. Especially when leaving out the veryinner parts of the disk, the structure of the disk changes drastically. Because of large orbitalvelocities in the central disk regions, explicit schemes are limited by the CFL condition and thusutilizing an implicit scheme can overcome this timestep restriction and produce more realisticlong term evolution results. References [1] L. Mestel, Problems of Star Formation - I, Quarterly Journal of the Royal Astronomical Society 6 (1965) 161.[2] E. Schatzman, A theory of the role of magnetic activity during star formation, Annales d’Astrophysique 25 (1962)18.[3] J.-M. Hameury, K. Menou, G. Dubus, J.-P. Lasota, J.-M. Hure, Accretion disc outbursts: a new version of an oldmodel, MNRAS 298 (4) (1998) 1048–1060. arXiv:astro-ph/9803242 , doi:10.1046/j.1365-8711.1998.01773.x .[4] R. Courant, K. Friedrichs, H. Lewy, ¨Uber die partiellen Di ff erenzengleichungen der mathematischen Physik, Math-ematische Annalen 100 (1928) 32–74. doi:10.1007/BF01448839 .[5] F. Masset, FARGO: A fast eulerian transport algorithm for di ff erentially rotating disks, Astronomy and Astro-physics, Supplement 141 (2000) 165–173. arXiv:astro-ph/9910390 , doi:10.1051/aas:2000116 .[6] M. Liska, C. Hesp, A. Tchekhovskoy, A. Ingram, M. van der Klis, S. Marko ff , Formation of precessing jets bytilted black hole discs in 3D general relativistic MHD simulations, Monthly Notices of the Royal AstronomicalSociety 474 (1) (2018) L81–L85. arXiv:1707.06619 , doi:10.1093/mnrasl/slx174 .[7] E. Vorobyov, V. Akimkin, O. Stoyanovskaya, Y. Pavlyuchenkov, H. B. Liu, The early evolution of viscousand self-gravitating circumstellar disks with a dust component arXiv:1801.06898 , doi:10.1051/0004-6361/201731690 .[8] J. Hern´andez, L. Hartmann, N. Calvet, R. D. Je ff ries, R. Gutermuth, J. Muzerolle, J. Stau ff er, A Spitzer View ofProtoplanetary Disks in the γ Velorum Cluster, ApJ 686 (2008) 1195–1208. arXiv:0806.2639 , doi:10.1086/591224 .[9] E. Dorfi, Computational Methods for Astrophysical Fluid Flow, Vol. 27 of Saas-Fee Advanced Courses, Springer-Verlag, Berlin / Heidelberg, 1998. doi:10.1007/3-540-31632-9 .[10] E. I. Vorobyov, A. M. Skliarevskii, V. G. Elbakyan, Y. Pavlyuchenkov, V. Akimkin, M. Guedel, Global evolu-tion of a gravitoviscous protoplanetary disk. I. The importance of the inner sub-au region, arXiv e-prints (2019)arXiv:1905.11335 arXiv:1905.11335 .[11] S. M. Andrews, J. Huang, L. M. P´erez, A. Isella, C. P. Dullemond, N. T. Kurtovic, V. V. Guzm´an, J. M. Car-penter, D. J. Wilner, S. Zhang, Z. Zhu, T. Birnstiel, X.-N. Bai, M. Benisty, A. M. Hughes, K. I. ¨Oberg, L. Ricci, he Disk Substructures at High Angular Resolution Project (DSHARP): I. Motivation, Sample, Calibration, andOverview arXiv:1812.04040 , doi:10.3847/2041-8213/aaf741 .[12] L. M. P´erez, M. Benisty, S. M. Andrews, A. Isella, C. P. Dullemond, J. Huang, N. T. Kurtovic, V. V. Guzm´an,Z. Zhu, T. Birnstiel, S. Zhang, J. M. Carpenter, D. J. Wilner, L. Ricci, X.-N. Bai, E. Weaver, K. I. ¨Oberg, TheDisk Substructures at High Angular Resolution Project (DSHARP): X. Multiple rings, a misaligned inner disk,and a bright arc in the disk around the T Tauri star HD 143006 arXiv:1812.04049 , doi:10.3847/2041-8213/aaf745 .[13] A. Isella, J. Huang, S. M. Andrews, C. P. Dullemond, T. Birnstiel, S. Zhang, Z. Zhu, V. V. Guzm´an, L. M. P´erez,X.-N. Bai, M. Benisty, J. M. Carpenter, L. Ricci, D. J. Wilner, The Disk Substructures at High Angular ResolutionProject (DSHARP) - IX. A high definition study of the HD 163296 planet forming disk arXiv:1812.04047 , doi:10.3847/2041-8213/aaf747 .[14] E. Dorfi, L. Drury, Simple adaptive grids for 1 - D initial value problems, Journal of Computational Physics 69 (1)(1987) 175–195. doi:10.1016/0021-9991(87)90161-6 .[15] L. D. Landau, E. M. Lifshitz, Fluid Mechanics, 1987.[16] D. Mihalas, B. Mihalas, Foundations of radiation hydrodynamics, 1984.[17] W. M. Tscharnuter, K. H. A. Winkler, A Method for Computing Self-gravitating Gas Flows with Radiation, Com-puter Physics Communications 18 (1979) 171–199. doi:10.1016/0010-4655(79)90111-5 .[18] P. J. Armitage, Astrophysics of Planet Formation, 2013.[19] L. Hartmann, G. Herczeg, N. Calvet, Accretion onto Pre-Main-Sequence Stars, Annual Review of Astronomy andAstrophysics 54 (2016) 135–180. doi:10.1146/annurev-astro-081915-023347 .[20] M. M. Romanova, G. V. Ustyugova, A. V. Koldoba, R. V. E. Lovelace, Launching of conical winds and axial jetsfrom the disc-magnetosphere boundary: axisymmetric and 3D simulations, MNRAS 399 (4) (2009) 1802–1828. arXiv:0907.3394 , doi:10.1111/j.1365-2966.2009.15413.x .[21] F. H. Shu, J. R. Najita, H. Shang, Z. Y. Li, X-Winds Theory and Observations, in: V. Mannings, A. P. Boss, S. S.Russell (Eds.), Protostars and Planets IV, 2000, pp. 789–814.[22] M. A. Belyaev, R. R. Rafikov, J. M. Stone, Angular Momentum Transport by Acoustic Modes Generated in theBoundary Layer. I. Hydrodynamical Theory and Simulations, Astrophysical Journal 770 (1) (2013) 67. arXiv:1212.0580 , doi:10.1088/0004-637X/770/1/67 .[23] S. A. Balbus, J. F. Hawley, A powerful local shear instability in weakly magnetized disks. I - Linear analysis. II -Nonlinear evolution, ApJ 376 (1991) 214–233. doi:10.1086/170270 .[24] D. N. C. Lin, J. E. Pringle, A viscosity prescription for a self-gravitating accretion disc, MNRAS 225 (1987)607–613. doi:10.1093/mnras/225.3.607 .[25] N. I. Shakura, R. A. Sunyaev, Black Holes in Binary Systems: Observational Appearances, in: H. Bradt, R. Giac-coni (Eds.), X- and Gamma-Ray Astronomy, Vol. 55 of IAU Symposium, 1973, p. 155.[26] X.-N. Bai, J. M. Stone, Magnetic Flux Concentration and Zonal Flows in Magnetorotational Instability Turbu-lence arXiv:1409.2512 , doi:10.1088/0004-637X/796/1/31 .[27] G. Lodato, W. K. M. Rice, Testing the locality of transport in self-gravitating accretion discs, Monthly Noticesof the Royal Astronomical Society 351 (2) (2004) 630–642. arXiv:astro-ph/0403185 , doi:10.1111/j.1365-2966.2004.07811.x .[28] A. C. Boley, A. C. Mej´ıa, R. H. Durisen, K. Cai, M. K. Pickett, P. D’Alessio, The Thermal Regulation of Grav-itational Instabilities in Protoplanetary Disks. III. Simulations with Radiative Cooling and Realistic Opacities,Astrophysical Journal 651 (1) (2006) 517–534. arXiv:astro-ph/0607112 , doi:10.1086/507478 .[29] P. Cossins, G. Lodato, C. J. Clarke, Characterizing the gravitational instability in cooling accretion discs, MonthlyNotices of the Royal Astronomical Society 393 (4) (2009) 1157–1173. arXiv:0811.3629 , doi:10.1111/j.1365-2966.2008.14275.x .[30] J. Binney, S. Tremaine, Galactic dynamics, 1987.[31] J. T. Conway, Analytical solutions for the Newtonian gravitational field induced by matter within axisymmet-ric boundaries, Monthly Notices of the Royal Astronomical Society 316 (2000) 540–554. doi:10.1046/j.1365-8711.2000.03523.x .[32] I. S. Gradshteyn, I. M. Ryzhik, Table of integrals, series and products, 1980.[33] I. A. Gerasimov, Calculation of complete elliptical integrals by the AGM-method., Trudy Gosudarstvennogo As-tronomicheskogo Instituta 60 (1988) 3–6.[34] B. van Leer, Towards the Ultimate Conservation Di ff erence Scheme. II. Monotonicity and Conservation Com-bined in a Second-Order Scheme, Journal of Computational Physics 14 (1974) 361–370. doi:10.1016/0021-9991(74)90019-9 .[35] K. H. Winkler, M. Norman, D. Mihalas, Adaptive-mesh radiation hydrodynamics - I. The radiation transport equa-tion in a completely adaptive coordinate system., Journal of Quantitative Spectroscopy and Radiative Transfer 31(1984) 473–489. doi:10.1016/0022-4073(84)90054-2 .[36] C. F. McKee, J. P. Ostriker, A theory of the interstellar medium: three components regulated by supernova explo- ions in an inhomogeneous substrate., ApJ 218 (1977) 148–169. doi:10.1086/155667 .[37] D. Lynden-Bell, J. E. Pringle, The evolution of viscous discs and the origin of the nebular variables., MNRAS 168(1974) 603–637. doi:10.1093/mnras/168.3.603 .[38] J. E. Pringle, Accretion discs in astrophysics, Annual review of astronomy and astrophysics 19 (1981) 137–162. doi:10.1146/annurev.aa.19.090181.001033 .[39] E. A. Dorfi, A. Gautschy, Where Are the Regularly Pulsating Massive Stars?, ApJ 545 (2000) 982–991. doi:10.1086/317861 ..