EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 606 (October 2017) A&A, 606 (2017) A5 Full HTML
Free Access
Issue
A&A
Volume 606, October 2017
Article Number A5
Number of page(s) 14
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201730516
Published online 22 September 2017

© ESO, 2017

1. Introduction

Circumstellar disks hold the key to our understanding of stellar mass accumulation and planet formation. They form thanks to the conservation of angular momentum of a collapsing cloud core when the inspiralling material hits the centrifugal barrier near the stellar surface before landing onto the star. As the core collapse continues, new layers of infalling material are added to the newborn disk at progressively larger radial distances, causing the disk to grow in both mass and size. Numerical simulations and observations indicate that this process starts in the Class 0 phase of stellar evolution and continue until the parental core depletes or dissipates (Machida et al. 2010; Tobin et al. 2015). In this protostellar or embedded phase, disks are usually most massive owing to continuing mass loading from the parental core and are often prone to the development of gravitational instability (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 the first phases of planet assembly are likely to take place (ALMA Partnership 2015).

In the subsequent T Tauri phase of stellar evolution, the disk slowly loses its mass owing to accretion onto the host star and expands in reaction to angular momentum redistribution within the disk. In this phase, dust growth from small grains to planetesimals and planetary cores takes place, finally leading to the emergence of protoplanets covered with primordial atmospheres (e.g., Benz et al. 2014) and basic chemical ingredients are converted into complex (organic) species (Henning & Semenov 2013). Finally, the combined action of stellar accretion, planet formation, disk winds, and photoevaporation leads to the dispersal of the disk gaseous component, revealing a debris disk consisting of solids that are left from the planet formation process.

Many of the above-mentioned phenomena during the entire disk evolution are controlled by the radiative input from the host star and external environment; both largely determine the temperature in the disk upper layers and set the minimum temperature in its midplane. In the steady-state disk models, the radiative input can be taken into account rather accurately using sophisticated ray-tracing or Monte Carlo techniques (e.g., Woitke et al. 2009; Dullemond 2012; Akimkin et al. 2013). However, because of the high computational costs, these models lack a dynamical aspect that can be important when considering young disks prone to gravitational instability or more evolved disks where planetesimals/planetary cores are subject to growth and migration. To properly consider these effects, fully dynamical models are needed; these models are usually based the equations of hydrodynamics complemented with a module that takes the radiative input from the stars and/or external environment into account.

In the full 3D hydrodynamics simulations, the radiative input can be taken into account by solving the equations of radiation transfer, usually with simplifying assumptions such as frequency integrated opacities and diffusion approximation (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, a 3D model is an inconvenient tool to make a parameter space study of the disk evolution for many model realizations and many orbital periods. For this purpose, 2D thin-disk models of disk evolution have been routinely employed (e.g., Vorobyov 2011; Zhu et al. 2012; Regály 2013; Gyergyovits et al. 2014) thanks to their low computational costs. In these models, the equation for thermal balance is usually complemented with cooling and heating functions that take the radiative input and local disk cooling into account. The form of these cooling and heating terms may differ slightly from study to study, but they are essentially based on the calculated midplane and surface temperatures of the disk and the optical depth from the disk surface to the midplane. The midplane temperature is usually set equal to the hydrodynamic temperature, the surface temperature is calculated assuming the blackbody character of the incident stellar and background radiation, and the optical depth is calculated assuming the vertically isothermal temperature distribution.

Although it has advantages in simplicity and low computational costs, this approach has obvious weaknesses. First, it is not clear whether the 2D hydrodynamic temperature is indeed representative of the disk midplane temperature. Second, this method provides little information on the disk vertical structure, essentially assuming that the disk is vertically isothermal, which may not be the case. For instance, passively irradiated disks are known to exhibit positive temperature gradients in the vertical direction (e.g., Dullemond et al. 2002), while dense gaseous clumps and spiral arms in gravitationally unstable disks may be characterized by a more complicated vertical temperature distribution (Vorobyov et al. 2014). To circumvent this difficulty, various forms of the vertical density distribution can be adopted (Gaussian, exponential), but this procedure cannot be easily applied to the vertical temperature distribution. As a result, the 2D models have limited applicability to studying chemical reactions that are known to sensitively depend on the gas temperature (which likely varies in the vertical direction as well).

This paper presents a method that addresses the above-mentioned weaknesses by means of coupling the gas dynamics computations in the thin-disk limit and calculations of the disk vertical structure. The gas dynamics in the disk plane is calculated using hydrodynamics equations, while the disk vertical structure is calculated using the equations of radiation transfer and hydrostatic balance in the vertical direction. As a result, we retrieve the full 3D density and temperature structure of the disk at every time step, which is lacking in purely 2D models. The presented method has lower time overhead than the purely 2D thin-disk models and it is faster than fully 3D models since the disk gravitational potential is found using a convolution theorem (see Eq. (4)), which is not applicable to the fully 3D models formulated in the curvilinear coordinate systems. The adopted 1D radiative transfer in the vertical direction is also much faster than the fully 3D method. Our method is therefore well suited for the long-term simulations of circumstellar disks. We compare the disk evolution calculated using this 2+1D approach with purely 2D simulations and briefly discuss the applicability of our new method to calculating the chemical evolution in circumstellar disks.

The paper is organized as follows. In Sect. 2, the hydrodynamics equations in the thin-disk limit are reviewed. In Sect. 3, we formulate the modifications made to the thin-disk model to improve the thermal balance calculations in the thin disk. In Sect. 4, we compare the disk evolution in the 2+1D and 2D approaches. The model caveats and future improvements are discussed in Sect. 5. The main conclusions are summarized in Sect. 6. The Appendix presents details of the solution procedure used to calculate the disk vertical structure.

2. Model equations in the thin-disk limit

The equations of mass, momentum, and energy transport describing the dynamics of circumstellar disks in the thin-disk limit can be formulated as follows: where Σ is the mass surface density, e is the internal energy per surface area, is the vertically integrated gas pressure calculated via the ideal equation of state as with γ = 7 / 5, is the velocity in the disk plane, and is the gradient along the planar coordinates of the disk. The gravitational acceleration in the disk plane, , takes into account the self-gravity of the disk and the gravity of the central protostar when formed. Disk self-gravity is found by solving for the Poisson integral (4)where rsc and rout are the radial positions of the computational inner and outer boundaries. This integral is calculated using a FFT technique which applies the 2D Fourier convolution theorem for polar coordinates and allows for the nonperiodic boundary conditions in the r-direction by effectively doubling the computation domain in this coordinate direction and filling it with zero densities (see Binney & Tremaine 1987, Sect. 2.8). Turbulent viscosity is taken into account via the viscous stress tensor Π, the expression for which can be found in Vorobyov & Basu (2010). The kinematic viscosity needed to calculate the viscous stress tensor is found adopting the Shakura and Sunyaev parameterization (Shakura & Sunyaev 1973), so that ν = αcsh, where is the sound speed and h is the disk scale height. The α-parameter is linked to the inferred strength of the magneto-rotational instability (MRI) following the method that takes the MRI active/inactive states into account as described in Bae et al. (2014).

We use the following form for the cooling term Λ in Eq. (3) based on the analytical solution of the radiation transfer equations in the vertical direction (Dong et al. 2016) (5)where τR = κRΣ1 / 2 and τP = κPΣ1 / 2 are the Rosseland and Planck optical depths to the disk midplane, κP and κR are the Planck and Rosseland mean opacities, and Σ1 / 2 = Σ / 2 is the gas surface density from the disk surface to the midplane. Similar forms of the disk cooling term were employed in other 1D axisymmetric and 2D thin-disk simulations (e.g., Kley & Crida 2008; Rice & Armitage 2009; Vorobyov & Basu 2010; Zhu et al. 2012; Bae et al. 2014) adopted from analytic studies of the disk vertical structure by Hubeny (1990) followed by slight modifications to the needs of numerical modeling by Johnson & Gammie (2003).

The heating function Γ is expressed by analogy to the cooling function as (Dong et al. 2016) (6)where Tirr is the irradiation temperature at the disk surface determined by the stellar and background blackbody irradiation as (7)where Tbg is the uniform background temperature (in our model set to the initial temperature of the natal cloud core) and Firr(r) is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central object. The latter quantity is calculated as (8)where γirr is the incidence angle of radiation arriving at the disk surface at radial distance r. The incidence angle is calculated using the disk surface curvature inferred from the radial profile of the disk vertical scale height (Vorobyov & Basu 2010). The total stellar luminosity L∗ includes contributions from the accretion and photospheric luminosities. Similar forms of the disk heating term were employed 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: the 2+1D approach

The numerical method described above is a fast and convenient tool for computing the disk evolution with high numerical resolution and for many physical realizations (e.g., Vorobyov & Basu 2010; Zhu et al. 2012). However, this method is essentially two-dimensional1 and, as such, it lacks the information on the disk vertical structure. Some studies (e.g., Dong et al. 2016) assume a Gaussian or exponential vertical density profile, but the same assumption cannot be easily made for the vertical temperature distribution. We therefore have developed a straightforward modification to this method, which enables a calculation of the density and temperature distributions in the vertical direction concurrently with the computations of the gas dynamics in disk plane.

3.1. The 2+1D approach

In this section, we formulate the modifications made to the thin-disk model in order to improve the thermal balance calculations in the disk. These modifications also enable a reconstruction of the disk vertical structure, thus providing information on the disk volumetric density and temperature distributions. The new equations of mass, momentum, and energy transport now read as follows: While Eqs. (9) and (10) remain essentially similar to their thin-disk counterparts (apart from the effect of stellar motion), Eqs. (11) now updates the internal energy that is only due to advection, viscous dissipation, and pressure work (adiabatic heating and cooling). To take the disk heating by the stellar and background irradiation and the disk cooling due to its own infrared emission into account, we solve the moment equations describing the propagation of diffuse IR radiation in the vertical direction written in the Eddington approximation where E is the radiation energy density, T the gas temperature, ρ the gas volume density, cV the heat capacity of the gas, c the speed of light, a the radiation constant, z the vertical distance from the midplane, σ the column density measured from the disk mid-plane, and S the heating source (per unit mass) by the stellar and interstellar UV radiation.

Equations (1212) and (13) are complemented with the equation describing the local vertical hydrostatic balance in the disk taking into account the gravity of the star as well as the local self-gravity of the disk (14)where M∗ is the mass of the central star and μ = 2.33 the mean molecular weight. We note that dσ = ρdz. We assume that the disk vertical columns at various positions in the disk do not influence each other, so that solving for Eqs. (1212)−(14) reduces to a series of 1D problems for each grid zone on the (r,φ) computational mesh.

Table 1

Model core parameters.

3.2. Method

Our method for solving Eqs. (9)−(14) consists of three steps. In the first step, called the source step, we update the gas velocity and internal energy (per unit surface area) due to gravity, viscosity, and pressure work by solving for the following equations:

In the second step, called the thermal step, we compute the change in the disk gas temperature due to radiative cooling/heating and reconstruct the disk vertical structure. To do this, we solve for the moment Eqs. (1212) and (13) describing the propagation of diffuse IR radiation in the vertical direction complemented with Eq. (14) for the vertical hydrostatic balance. We note that in step 2 we use the gas temperatures that are partly updated in step 1. The detailed solution procedure for step 2 is provided in the appendix.

In the third step, called the transport step, we take advection of hydrodynamical quantities into account by solving for the following equations: In this final step, we use the gas velocities and internal energies which are consequently updated during steps 1 and 2.

To accomplish steps 1 and 3, we employ a combination of the finite difference and finite volume methods with a time-explicit solution procedure similar in methodology to the ZEUS code (Stone & Norman 1992). The advection in step 3 is treated using the third-order-accurate piecewise parabolic scheme of Colella & Woodward (1984). A small amount of artificial viscosity is added to the code to smooth out shocks over two grid zones in both coordinate directions. The associated artificial viscosity torques integrated over the disk area are negligible in comparison with gravitational or turbulent viscosity torques.

Equations (15), (16), and (17)−(19) are discretized in polar coordinates (r,φ) on a numerical grid with 512 × 512 grid zones. The radial points are logarithmically spaced, while the azimuthal points are equidistant. The innermost grid point is located at the position of the sink cell rsc = 10 AU, and the size of the first adjacent cell is 0.14 AU, which corresponds to a radial resolution of Δr = 1.4 AU at 100 AU. With this grid spacing, the Jeans length RJ = ⟨ v2 ⟩ /πGΣ, where ⟨ v2 ⟩ is the velocity dispersion in the disk plane (Vorobyov 2013), is resolved by roughly 10–20 grid zones in each coordinate direction to radial distance up to 500 AU, thus fulfilling the Truelove criterion (Truelove et al. 1998).

Equations (1212)–(14) are solved on an adaptive, non-equidistant grid with 32 grid points. The finest grid spacing is usually in the disk atmosphere where the largest density gradients are found. The inner and outer boundaries on the polar grid (r,φ) allow material to freely flow out from the active computational domain, but prevent any material from flowing in. In the vertical direction, we assume a reflecting boundary in the disk midplane and a constant gas volume density of 103 cm-3 at the disk upper edge.

3.3. Initial conditions

Our numerical simulations start from a pre-stellar core with the radial profiles of column density Σ and angular velocity Ω described as where Σ0 and Ω0 are the gas surface density and angular velocity at the center of the core. These profiles have a small near-uniform central region of size r0 and then transition to an r-1 profile; they are representative of a wide class of observations and theoretical models (André et al. 1993; Dapp & Basu 2009).

Our iterative solution procedure for Eqs. (1212)–(14) (see the appendix), requires us to make an initial guess for the vertical structure of the core. We assume a constant gas temperature of 15 K and a Gaussian distribution of the gas volume density of the form (22)where . We note, however, that these initial conditions are a mere initialization requirement and the vertical structure of the core quickly attains the form determined by the combined action of the external heating, radiative cooling, pressure gradients, and self-gravity of the core.

We considered several model cores, but in this paper, for the sake of conciseness, we present only two. The initial parameters of these models are shown in Table 1, where Mc is the initial core mass, β the ratio of rotational energy to the magnitude of gravitational potential energy, and rout the initial radius of the core. The two models are different in the amount of initial rotation, as manifested by distinct β and Ω0. The initial parameters are chosen so that the cores are gravitationally unstable and begin to collapse at the onset of numerical simulations. We monitor the gas surface density in the sink cell, and when its value exceeds a critical one for the transition from isothermal to adiabatic evolution, we introduce a central point-mass star. In the subsequent evolution, 90% of the gas that crosses the inner boundary is assumed to land on the growing star. The other 10% of the accreted gas is assumed to be carried away with protostellar jets The simulations continue into the embedded phase of star formation, during which a protostellar disk is formed. The simulations are terminated in the T Tauri phase when nearly all material of the parental core has accreted onto the resulting star-plus-disk system.

thumbnail Fig. 1

Gas surface density distributions in model A+ shown at several time instances after the onset of numerical simulations (the disk forms at t = 0.1 Myr). The insets zoom in on several fragments in the disk. The yellow contour lines delineate regions with the Toomre Q-parameter smaller than unity. The red dashed line is a radial cut through the disk used later to show the vertical volume density and temperature distributions. The scale bar is in log g cm-2.

4. Comparison of 2+1D and 2D approaches

In this section, we compare the properties of circumstellar disks obtained using the original 2D thin-disk model and the improved 2+1D model with the purpose of understanding how an improved calculation of the thermal balance in the disk can effect its dynamical evolution. First, we compute the disk evolution using the 2+1D models starting from the collapse of pre-stellar cores and ending in the T Tauri stage of disk evolution when the disk age exceeds 0.5 Myr. Then, we compute the disk evolution in the 2D model starting from a certain time instance, usually in the Class I phase, using the current values of Σ, e, and v taken from the 2+1D simulations. This allows us to directly determine the effect of different thermodynamical schemes on the disk dynamical evolution and avoid the possible dynamical influence of the early collapse phase. In the following text, the models computed using the 2+1D approach are distinguished with the “plus” sign.

thumbnail Fig. 2

Similar to Fig. 1, but for model A.

4.1. Models A+ and A

Figure 1 presents the disk evolution in model A+ computed using the 2+1D approach. Shown are the gas surface density snapshots in the inner 2000 × 2000 AU2 box taken at several times since the onset of numerical simulations. The entire numerical domain is about 10 times larger and includes the infalling envelope. Clearly, the disk is strongly gravitationally unstable and multiple fragments are seen forming in the disk’s middle and outer regions. The Toomre Q-parameter, Q = csΩ /πGΣ, is smaller than unity in and around the fragments, as is demonstrated by the yellow contour lines in the insets of Fig. 1. The masses of these fragments range from about a Jupiter-mass to the upper limit on the brown dwarf mass.

This behavior is typical for massive nonmagnetized or weakly-magnetized disks in the embedded phase of evolution (especially, in the Class I phase), where gravitational instability is fueled by continuous mass loading from 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 mass ratio ~ 0.2 at the end of the embedded phase. Many fragments do not live long and migrate into the star, causing FU Orionis-type luminosity outbursts (Vorobyov & Basu 2005, 2015; Machida et al. 2011). Others are dispersed by tidal torques (Vorobyov 2011; Zhu et al. 2012) or are ejected from the disk via multi-body gravitational interactions (Stamatellos & Whitworth 2009; Basu & Vorobyov 2012; Vorobyov 2016). However, some fragments may survive, contract to planetary sizes, and form planets or brown dwarfs at various distances from the star (Boss 1998; Nayakshin 2010, 2011; Boley et al. 2010; Vorobyov et al. 2013; Stamatellos & Herczeg 2015a,b; Galvagni & Mayer 2014). Disk fragmentation, therefore, can be an important channel for the formation of planets and brown dwarfs, either as companions to the host star or as freely floating objects.

Figure 2 presents the disk evolution in model A computed using the 2D approach. The gas surface density maps at the same evolutionary times as in model A+ are shown in the inner 2000 × 2000 AU2 box. The 2D simulations start from t = 0.25 Myr, explaining why the disk appearance at t ≤ 0.25 Myr is identical in models A+ and A. At later times, however, notable differences in the disk evolution between model A+ and model A arise. While the disk in model A+ is strongly unstable and often harbors multiple fragments, the disk in model A experiences only occasional fragmentation after t = 0.25 Myr.

This difference is further illustrated in Fig. 3 showing the number of fragments in the disk Nfrag at a certain time instance as a function of time elapsed since the onset of numerical simulations. Two conditions were used to identify the fragments in the disk (see Vorobyov 2013, for detail): they must be pressure supported, with a negative pressure gradient with respect to the center of the fragment, and they must be kept together by gravity, with the potential well being deepest at the center of the fragment. Clearly, the number of fragments in the disk at a given time instance is greater in model A+ than in model A. In addition, the duration of the disk fragmentation stage (defined as the time span during which fragments are present in the disk) is longer in model A+. Both suggest that the disk in model A+ is more gravitationally unstable.

thumbnail Fig. 3

Number of fragments in the disk at a certain time instance as a function of time elapsed since the onset of numerical simulations in model A+ (top) and model A (bottom). The disk forms at t = 0.1 Myr.

This difference in the strength of gravitational instability can in principle be caused by distinct disk and stellar masses in models A+ and A. It is known that systems with a higher disk-to-star mass ratio are characterized by stronger gravitational instability because they have on average lower angular velocities and/or higher surface densities, and are thus characterized by lower Q-parameters. However, we find that model A+ has a slightly lower disk mass and a slightly higher stellar mass than model A, meaning that the mass transport rate through the disk in this model is systematically higher than in model A. The higher mass transport rates in model A+ are caused by systematically higher gravitational torques, as can be expected for disks with stronger gravitational instability.

The stronger gravitational instability in model A+ compared to that in model A is not related to the disk or stellar masses. It may then be related by differences in the disk thermal structure arising from distinct approaches to calculating the disk cooling and heating in 2+1D and 2D approaches. To check if this is indeed the case, we plot in Fig. 4 the radial gas temperature profiles in model A+ (red lines) and model A (black lines), calculated as Thydro = e(γ − 1)μ/ Σℛ for every cell on the polar grid (r,φ) and then arithmetically averaged over the azimuth. We note that in the 2+1D approach, Thydro represents the gas temperature mass-weighted over the vertical column of gas in each cell (r,φ) (see Eq. (A.17)). For model A+, we also plot the temperature in the disk midplane Tmp (blue lines). Four time instances since the onset of the 2+1D and 2D simulations are shown.

thumbnail Fig. 4

Azimuthally averaged radial profiles of the gas temperature in model A+ (red and blue lines) and model A (black solid lines) taken at several times since onset of numerical simulations. In particular, the red lines represent the gas temperature mass-weighted over the vertical column of gas, while the blue lines are the gas temperature in the disk midplane. The dashed solid lines are the gas temperatures in model A derived for the steady-state case neglecting compressional and viscous heating. The vertical dotted lines mark the radial positions where the disk joins the infalling envelope.

Clearly, notable differences in the radial gas temperature distributions develop with time in models A+ and A. The disk in model A+ is systematically colder than in model A,no matter what temperature in model A+ is considered, the hydrodynamic (Thydro) or the midplane (Tmp). On the contrary, the inner envelope is warmer in model A+ than in model A. These differences can be understood from the following analysis. Let us first consider the 2D case and for a moment neglect the heating sources due to viscosity and compressional heating due to PdV work. In the steady-state case, the temperature in the disk and envelope will be controlled by a balance between radiative cooling Λ (Eq. (5)) and irradiation and background heating Γ (Eq. (6)), so that the midplane temperature can be written as (23)This steady-state temperature is plotted in Fig. 4 with the dashed black lines. Clearly, it is very close to the actual temperature in model A everywhere except the very inner parts of the disk where viscous and compressional heating become substantial and the midplane temperature rises above the analytically predicted values.

In the 2+1D case, the midplane temperature can be expressed in the following form if the medium is optically thick to ultraviolet and infrared radiation (24)This expression follows from the assumption that the incident UV flux absorbed by dust in the disk is radiated away isotropically to the upper and lower hemisphere (see Eq. (12a) in Chiang & Goldreich 1997).

A comparison of Eqs. (23) and (24) in the regions where the input from the background irradiation is much smaller than from the stellar irradiation (i.e., in the disk), demonstrates that the midplane temperature in model A+ is expected to be a factor of 21 / 4 ≈ 1.2 smaller than that of model A. A similar difference between the gas temperatures in models A+ and A is also seen in Fig. 4. We note that the stellar luminosity in our models exhibits short-term variations caused by the time-variable protostellar accretion (e.g., Vorobyov & Basu 2015). As a result, the disk thermal state may take some time to adjust to the time-variable stellar flux, meaning that the midplane temperature may not exactly coincide with the analytical values derived in the steady-state limit of constant stellar luminosity.

In the optically thin (to UV radiation) limit, the midplane temperature in the 2+1D model can be written as (25)where κF is the dust opacity in the UV band. Clearly, the midplane temperature in the optically thin limit is higher than in the optically thick limit by a factor of (κF/κR)1 / 4. The optically thin limit is expected to take place in the envelope and this explains an increase in the gas temperature in model A+ at distances > 1000 AU. In the 2D case, the gas temperature in the envelope is controlled by the background irradiation with Tbg = 15 K. At the same time, Eq. (25) for (only the interstellar component, see Appendix A) yields the midplane temperature ≈ 14 K for κF/κR = 103, demonstrating that the midplane temperatures in the two models converge to a similar value in the outer envelope.

The consequences of this distinct radial temperature distribution is twofold. First, a lower gas temperature makes the disk in model A+ more gravitationally unstable, as was already noted above. Second, higher infall rates from the collapsing envelope onto the disk , as implied by a higher gas temperature in the inner envelope, drive the disk in model A+ faster to the critical point at which the disk becomes unstable to fragmentation, thus creating more fragments in the disk. In other words, the characteristic time of mass infall onto the disk (or the disk mass replenishment) Mdisk/Ṁinfall is shorter in model A+.

thumbnail Fig. 5

Distributions of gas volume density (top) and temperature (bottom) in the (r,z) plane taken in model A+ at t = 0.25 Myr through the radial cut with a position angle of 351° (red dashed line in Fig. 1).

In Fig. 5 we present the 2D distributions of the gas volume density and temperature in the (r,z) plane obtained in model A+ at t = 0.25 Myr by taking a vertical cut with a position angle of 351° (shown in Fig. 1 with the red dashed line). The cut is chosen so that it passes through an inner spiral arm at r ≈ 48 AU and a dense fragment with mass ≈ 7.8 MJup at r ≈ 207 AU (shown in the inset of the upper middle panel in Fig. 1). Clearly, the gas volume density is highest in the midplane and drops down toward the disk atmosphere. At the same time, the temperature is higher in the upper layers and is minimal in the disk midplane everywhere in the disk, except for the fragments. This means that the temperature in all parts of the disk, except for the fragments, is controlled by the external radiation (stellar and background). In the fragments, however, the gas is intensively heated by compressional and viscous heating, while the radiative diffusion is not fast enough to cool the inner optically thick parts of the fragment down to low temperatures observed around the fragment. We also note that the fragment at r = 207 AU is gravitationally bound, which results in a very low disk vertical height at and around this point. The signature of disk self-gravity is also seen at r = 230 − 250 AU covering the spiral arc in which the fragment resides and where the disk vertical height is also significantly reduced compared with the surroundings. It is therefore possible that fragments may be deeply hidden in the disk, which would make their detection via observations of the near-infrared light scattered off the disk surface quite difficult, as was already noted in Dong et al. (2016). Fully 3D simulations are needed to determine how far from the midplane the fragments can be scattered via gravitational interaction with other fragments and spiral arms.

thumbnail Fig. 6

Temperature and gas surface density distributions in model A+. Various definitions of temperature as described in the text are shown.

The 2D temperature distributions in the disk midplane and at the surface of the disk are shown for model A+ in the upper left and upper right panels of Fig. 6. The bottom left panel presents the gas hydrodynamic temperature defined as Thydro = e(γ − 1)μ/ (Σℛ), in a similar fashion to that of the midplane temperature in the 2D approach (see Eq. (5)). The gas surface density distribution is also shown in the bottom right panel for convenience. All distributions are plotted at t = 0.25 Myr.

The hydrodynamic and midplane temperatures reflect the structure of the disk: both are rather low in the outer disk regions and dense spiral arms, but grow in the inner disk thanks mainly to increased stellar irradiation. Dense fragments heated by compressional and viscous heating stand out as bright spots in the disk. The hydrodynamic temperature Thydro appears to be higher than the midplane temperature Tmp in the disk outer regions and in the inner envelope. This trend can also be seen in Fig. 4 and is explained by the passively heated (by stellar and interstellar irradiation) nature of the outer disk and envelope: the upper gas layers are always warmer than the midplane. As a result, Thydro – a vertical average over the gas column – becomes higher than Tmp. In the inner disk regions, where efficient viscous and compressional heating operate in the disk midplane, the situation may reverse and Tmp may become higher than Thydro. This trend is also evident in Fig. 4. The surface temperature smoothly decreases with radial distance from the star. The lack of azimuthal variations is a mere consequence of the adopted scheme for calculating the incidence angle of radiation onto the disk surface, which uses an azimuthally averaged disk scale height to calculate γirr. We note that the surface temperature is higher than the midplane temperature almost everywhere in the disk except for the fragments, which is discussed in more detail below.

Finally, in Fig. 7 we show the 1D vertical gas volume density and temperature profiles taken at several positions in the radial cut (the red dashed line in Fig. 1). These positions were chosen so as to show the vertical distributions in various sub-structures that may be present in the disk, such as spiral arms and fragments. More specifically, the blue lines represent the vertical profiles taken through the fragment shown in the inset of the upper middle panel in Fig. 1, while the cyan lines are taken through the spiral arm. The black lines represent the vertical profiles in a regular inter-arm disk region, and the red lines are the vertical profiles in the inner disk. We note that we apply an adaptive grid in the vertical direction that enables an adequate numerical resolution in the disk regions with strong volume density gradients and also in the disk atmosphere.

Clearly, the fragment is characterized by the highest volume density amounting to 1012 cm-3. At the same time, the fragment has the most compact structure with a vertical size of just 16 AU, even though it is actually located quite far away from the star. The midplane size of this fragment is ≈ 30 AU so that the clump has an ellipsoidal shape with the semi-axis ratio of az:ar = 1:2. This scaling is in agreement with the ratio of rotational-to-thermal energy of 44%, implying that the fragment has a substantial rotational support against gravity (in case of zero rotation, we would have expected a near-spherical shape). The fragment is also characterized by a peculiar vertical temperature distribution: its midplane temperature is higher than the temperature in its atmosphere. There is also a depression in the temperature profile in between the disk midplane and atmosphere. A similar irregular vertical temperature distribution can also found in the inner disk regions (the red line). In this case, however, the temperature in the atmosphere is somewhat higher than in the disk midplane. The irregular temperature profiles are a direct consequence of the compressional and viscous heating operating in the optically thick interiors of the fragment and in the inner disk regions. These heating sources are more efficient than heating due to stellar irradiation, the latter being effectively absorbed by the atmosphere of the fragment. This is an important phenomenon, which means that radiation transfer codes that neglect hydrodynamic heating sources (such as RADMC-3D) cannot accurately reproduce the temperature structure in the fragments formed via disk gravitational 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 temperature distribution that is typical for passively irradiated disks, though the spiral arm demonstrates a mild increase in the gas temperature towards the midplane, caused probably by mild shock and viscous heating.

thumbnail Fig. 7

One-dimensional vertical cuts taken at different positions of the disk (as indicated in the legend) and showing the vertical gas volume density (top) and temperature (bottom) distributions.

4.2. Models B+ and B

Model B+ is distinct from model A+ in the amount of rotation initially present in the parental collapsing core (see Table 1). More specifically, model B+ has a lower ratio of rotational to gravitational energy β, which implies that model B+ forms a less massive disk than model A+. Following the same procedure as was discussed in the beginning of Sect. 4, we first computed the disk evolution in model B+ and then we computed model B starting from a time instance t = 0.215 Myr, which corresponds to the Class I phase of disk evolution.

thumbnail Fig. 8

Number of fragments in the disk at a certain time instance as a function of time elapsed since the onset of numerical simulations in model B+ (top) and model B (bottom). The disk forms at t = 0.138 Myr.

A comparison of model B+ and model B has revealed trends that are very similar to those discussed in detail for model A+ and model A. The disk in model B+ is more gravitationally unstable, more prone to fragmentation, and, at the same time, is less massive than in model B. This apparent paradox is explained by a systematically colder disk temperatures in model B+ than in model B. For the sake of conciseness, we do not perform an in-depth analysis of models B+ and B in this section, but simply present in Fig. 8 the number of fragments in the disk Nfrag at a given time instance as a function of time elapsed since the onset of numerical simulations. Clearly, Nfrag is larger in model B+ than in model B. There can be as many as nine fragments in model B+ in the early evolutionary phase, whereas in model B the number of fragments never exceeds four. The duration of the disk fragmentation phase is also greater in Model B+, which implies a higher detection likelihood of disk fragmentation.

5. Model caveats and further improvements

In this section, we discuss the model caveats and further improvements of our 2+1D model.

Disk self-gravity. When reconstructing the disk vertical structure, a local value for the gas column density σ was used in the last term of Eq. (14). Under this approximation, the vertical gas columns do not affect each other when solving the equation of the vertical hydrostatic balance. This may affect the disk vertical structure in the regions dominated by gravity. For instance, a sharp drop in the disk height seen around the fragment in the upper panel of Fig. 5 may be an artifact of the adopted approximation. In the future, the model needs to be improved by taking the non-locality of disk gravity into account.

Vertical motions. Another assumption behind the presented model is that vertical motions are neglected. Specifically, we assume that the disk attains a vertical hydrostatic equilibrium on timescales much shorter than the dynamical timescale. This approximation is justified if the dynamical timescale is longer than any other pertinent timescale, such as the timescales for diffusion of radiation, propagation of sound waves, and stellar heating, which was shown to usually be the case for circumstellar disks at distances greater than 1 AU in Vorobyov et al. (2014). With the current approach, we cannot model interesting effects, such as vertical oscillations and surface waves. This is, however, the price that we are willing pay for the ability to follow the disk evolution on much longer timescales than is currently possible with the full 3D numerical hydrodynamics simulations.

Realistic equation of state. In the current paper, we did not take into account that the heat capacity cV and the adiabatic index γ of gas depend on the temperature since the rotation levels of molecular hydrogen can be excited in the considered temperature ranges. We did not include this effect since it requires a more complicated procedure for the solution of the radiation transfer equations. We plan to do this in follow-up papers.

There are other modifications to our radiative transfer model which can be implemented in the future. Currently, to calculate heating due to the stellar UV radiation, we evaluate the incidence angle of stellar radiation γirr (see Eq. (A.4)) in each grid cell based on the local gradients of the disk scale height, and then average the resulting values over the azimuth (Vorobyov & Basu 2010). This procedure provides the radially varying UV heating due to the global disk flaring, but does not allow us to properly model the local effects, such as shadows in the disk caused by spiral arms. Our attempt to calculate the azimuthally varying γirr, including the effect of the shadows, resulted in overheating of the disk surfaces directly exposed to stellar radiation. Allowing for the diffusion of the thermal IR emission in the equatorial plane may alleviate this problem. Finally, we now work on implementing the FARGO mechanism for advection of hydrodynamical quantities (Masset 2000), which will allow us to ease the restrictive time step limitations and move the sink cell boundary closer to the star, 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 presented the method that includes a better calculation of the disk thermal balance and a reconstruction of the disk vertical structure. Our method is based on the solution of the hydrodynamics equations describing the gas dynamics in the plane of the disk, complemented with solution of the radiation transfer and hydrostatic balance equations describing the disk vertical structure. We performed a detailed comparison of this 2+1D method with the purely 2D thin-disk models of disk evolution. Our findings can be summarized as follows:

  • Improved 2+1D models yield systematically colder disks, whilethe infalling parental clouds in the embedded evolutionary phaseare warmer. Both effects act to increase the strength of diskgravitational instability in 2+1D models as compared to purely2D simulations.

  • Disk gravitational fragmentation is more efficient and the duration of the disk fragmentation phase is longer in 2+1D models, which implies an increased likelihood for detecting disk fragmentation in protostellar disks.

  • The outer disk regions are mostly characterized by a positive vertical temperature gradient, typical for passive circumstellar disks heated mainly by stellar and background irradiation. The inner disk regions usually have a more complex, non-regular vertical temperature distribution having local peaks in the midplane and in the atmosphere separated by a mild depression. The temperature increase in the midplane is caused by efficient viscous and compressional heating operating in the inner disk.

  • Fragments forming in the disk via gravitational fragmentation are also characterized by a double-peaked vertical temperature profile, but unlike the inner disk regions, the center of the fragment is warmer than its atmosphere. This means that the hydrodynamical heating sources (compression, viscosity) are more efficient than stellar irradiation heating, implying that radiation transfer codes that neglect the hydrodynamical heating sources cannot accurately compute the temperature in the fragment interiors (see also Dong et al. 2016).

  • Fragments are significantly more compact than the surrounding disk, which could make their detection with the scattered light techniques difficult unless significant excursions away from the disk midplane are feasible.

A detailed procedure for solving the radiation transfer and hydrostatic balance equations in the vertical direction is presented in the appendix. The improved 2+1D method yield a full 3D structure of the disk, namely its volumetric density and temperature distributions, with a modest increase in the net computational time with respect to purely 2D models. It also applies an adaptive mesh in the vertical direction, enabling good resolution in the disk regions with strong density gradients and in the disk atmosphere. This makes it possible to couple our 2+1D model with compact chemical reaction networks to follow the long-term chemodynamical evolution of young protostellar disks starting from the gravitational collapse of parental cloud cores. The details of this method will be presented in the upcoming paper.

1

Some modifications include a calculation of the disk vertical scale height and the incidence angle of stellar irradiation (Vorobyov & Basu 2010).

Acknowledgments

We thank the anonymous referee for the constructive comments that helped us to improve the manuscript and the vertical reconstruction procedure. This work was supported by RFBR grant 17-02-00644. The simulations were performed on the Vienna Scientific Cluster (VSC-2), on the Shared Hierarchical Academic Research Computing Network (SHARCNET), on the Atlantic Computational Excellence Network (ACEnet). This publication is supported by the Austrian Science Fund (FWF).

References

  1. Akimkin, V. V., Pavlyuchenkov, Y. N., Launhardt, R., & Bourke, T. 2012, Astron. Rep., 56, 915 [NASA ADS] [CrossRef] [Google Scholar]
  2. Akimkin, V., Zhukovska, S., Wiebe, D., et al. 2013, ApJ, 766, 8 [NASA ADS] [CrossRef] [Google Scholar]
  3. ALMA Partnership,Brogan, C. L., Perez, L. M., Hunter, T. R., et al. 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
  4. André, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baraffe, I., Elbakyan, V. G., Vorobyov, E. I., & Chabrier, G. 2017, A&A, 597, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Basu, S., & Vorobyov E. I. 2012, ApJ, 750, 30 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  9. Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 691 [Google Scholar]
  10. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Univ. Press) [Google Scholar]
  11. Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H. 2010, Icarus, 207, 509 [NASA ADS] [CrossRef] [Google Scholar]
  12. Boss A. P. 1998, ApJ, 503, 923 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chiang, E., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
  14. Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  15. Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092 [NASA ADS] [CrossRef] [Google Scholar]
  16. Dong, R., Vorobyov, E. I., Pavlyuchenkov, Y., Chiang, E., & Liu, H. B. 2016, ApJ, 823, 141 [NASA ADS] [CrossRef] [Google Scholar]
  17. Dullemond, C. 2012, Astrophysics Source Code Library, record [tascl:1202.015] [Google Scholar]
  18. Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Elbakyan, V. G., Vorobyov, E. I., & Glebova, G. 2016, Astron. Rep., 60, 879 [NASA ADS] [CrossRef] [Google Scholar]
  20. Galvagni, M., & Mayer, L. 2014, MNRAS, 437, 2909 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gyergyovits, M., Eggl, S., Pilat-Lohinger, E., & Theis, Ch. 2014, A&A, 566, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Henning, Th., & Semenov, D. 2013, Chem. Rev., 113, 9016 [Google Scholar]
  23. Hubeny, I. 1990, ApJ, 351, 632 [NASA ADS] [CrossRef] [Google Scholar]
  24. Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131 [NASA ADS] [CrossRef] [Google Scholar]
  25. Klahr, H., & Kley, W. 2005, A&A, 445, 747 [Google Scholar]
  26. Kley, W., & Crida, A. 2008, A&A, 487, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 681, 375 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Lomax, O., Whitworth, A. P., Hubber, D. A., Stamatellos, D., & Walch, S. 2014, MNRAS, 439, 3039 [NASA ADS] [CrossRef] [Google Scholar]
  30. Machida, M. N., Inutsuka, S., & Matsumoto, T. 2010, ApJ, 724, 1006 [NASA ADS] [CrossRef] [Google Scholar]
  31. Masunaga, H., Miyama, S. M., & Inutsuka, S. 1998, ApJ, 495, 346 [Google Scholar]
  32. Masset, A. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Mayer, L., Peters, T., Pineda, J. E., Wadsley, J., & Rogers, P. 2016, ApJ, 823, 36 [Google Scholar]
  34. Meyer, D. M.-A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017, MNRAS, 464, 90 [Google Scholar]
  35. Nayakshin, S. 2010, MNRAS, 408, 2381 [NASA ADS] [CrossRef] [Google Scholar]
  36. Nayakshin, S. 2011, MNRAS, 413, 1462 [NASA ADS] [CrossRef] [Google Scholar]
  37. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 [NASA ADS] [CrossRef] [Google Scholar]
  38. Pavlyuchenkov, Y. N., Zhilkin, A. G., Vorobyov, E. I., & Fateeva, A. M. 2015, Astron. Rep., 59, 133 [NASA ADS] [CrossRef] [Google Scholar]
  39. Regály, Zs., & Vorobyov, E. I. 2017, A&A, submitted [Google Scholar]
  40. Regály, Zs., Sándor, Zs., Csomós, P., & Ataiee, S. 2013, MNRAS, 433, 2626 [NASA ADS] [CrossRef] [Google Scholar]
  41. Rice, W. K. M., & Armitage, P. J. 2009, MNRAS, 396, 2228 [NASA ADS] [CrossRef] [Google Scholar]
  42. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  44. Stamatellos, D. 2015, ApJ, 810, L11 [NASA ADS] [CrossRef] [Google Scholar]
  45. Stamatellos, D., & Herczeg, G. J. 2015, MNRAS, 449, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  46. Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413 [Google Scholar]
  47. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [Google Scholar]
  48. Tobin, J. J., Looney, L.W., Wilner, D. J., et al. 2015, ApJ, 805, 125 [NASA ADS] [CrossRef] [Google Scholar]
  49. Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1998, ApJ, 495, 821 [NASA ADS] [CrossRef] [Google Scholar]
  50. Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2013, MNRAS, 436, 1667 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  51. Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 446, 1175 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  52. Vorobyov, E. I. 2010, ApJ, 723, 1294 [NASA ADS] [CrossRef] [Google Scholar]
  53. Vorobyov, E. I. 2011, ApJ, 728, L45 [NASA ADS] [CrossRef] [Google Scholar]
  54. Vorobyov, E. I. 2013, A&A, 552, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Vorobyov, E. I. 2016, A&A, 590, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137 [NASA ADS] [CrossRef] [Google Scholar]
  57. Vorobyov, E. I., & Basu, S. 2009, MNRAS, 393, 822 [NASA ADS] [CrossRef] [Google Scholar]
  58. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  59. Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115 [NASA ADS] [CrossRef] [Google Scholar]
  60. Vorobyov, E. I., Zakhozhay, O. V., & Dunham, M. M. 2013, MNRAS, 433, 3256 [NASA ADS] [CrossRef] [Google Scholar]
  61. Vorobyov, E. I., Pavlyuchenkov, Y. N., & Trinkl, P. 2014, Astron. Rep., 58, 522 [NASA ADS] [CrossRef] [Google Scholar]
  62. Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The thermal and step the reconstruction of the disk vertical structure

The thermal step is placed in our algorithm between the source step (Eqs. (15) and (16)) and the transport step (Eqs. (17)−(19)). While the source and transport steps deal with 2D quantities integrated over the disk vertical column (such as the gas surface density Σ, the vertically integrated pressure , and the internal energy per surface area e) and are therefore defined on the (r,φ) polar mesh, the thermal step deals with 3D quantities (such as the gas temperature T, the gas volume density ρ, the radiative energy E, and the vertical column density from the disk midplane σ) defined on the (z,r,φ) cylindrical mesh. The initial values for all quantities (2D and 3D) are provided by the initial setup as described in Sect. 3.3.

Here, we describe the algorithm for solving Eqs. (1212) and (13), which take the disk radiative cooling/heating into account. This algorithm also includes the reconstruction of the disk vertical structure assuming the vertical hydrostatic equilibrium described by Eq. (14). Our method is a time-dependent modification to the steady-state models presented in Akimkin et al. (2012) and Vorobyov et al. (2014). The method for solving the source and transport steps are described in Sect. 3.2 and in more detail in Vorobyov & Basu (2010).

First, we note that the internal energy per surface area e is updated in the source step, due to the work and viscous heating. This may affect the gas temperature distribution in the disk and should be taken into account before the thermal step is commenced. We assume that the generated (or consumed) heat in the source step is evenly redistributed over the disk vertical column so that the gas temperature in each element of the vertical column can be updated as (A.1)where index n refers to the quantities at the beginning of the source step and the asterix refers to the quantities after the source step. The updated gas temperature T∗ is also used in the thermal step.

In the second step, we compute the heating function S (Eq. (1212)) due to sources other than the thermal radiation of the medium. We note that the dimension of S by definition is [ erg s-1 g-1 ]. In our model, these sources include the UV radiation from the central star Sstar and the interstellar UV radiation Sbg, (A.2)where Sstar is defined as (A.3)where is the mean Planck opacity averaged over the stellar spectrum.

To compute Sstar, we need to know the distribution of the mean intensity of UV radiation Jstar. We calculate Jstar taking into account the absorption of the UV radiation by the disk as (A.4)where 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 stellar radiation arriving at the disk surface with respect to the normal (see Vorobyov & Basu 2010, for details). The optical depth as a function of the current column density σ is defined as (A.5)where σ is measured from the disk mid-plane. The intensity at the disk surface can be found as (A.6)Here, L∗ is the total stellar luminosity, which includes contributions from photoshperic luminosity Lph and accretion luminosity Laccr, and r the radial distance to the star. The accretion luminosity is calculated as Lacc = GM∗Ṁ/ (2R∗) using information on the current stellar mass M∗, stellar radius R∗, and accretion rate onto the star Ṁ. The photospheric luminosity and stellar radius are found from the Lyon stellar evolution code coupled to the main hydrodynamics code as described in Vorobyov & Basu (2015) and Baraffe et al. (2017).

The background heating function Sbg is calculated by analogy to Eq. (A.4), but with a fixed value of μ = 0.5. As a boundary condition for the interstellar radiation, we use the relation (A.7)where c is the speed of light and a the radiative constant. The temperature and dilution of the interstellar radiation are Tbg and D; in our simulations we adopt 10 000 K and 8 × 10-15, respectively.

In the third step, we calculate the change in the gas temperature T and radiative energy E in a given vertical column of the disk due to heating by the stellar UV and background radiation and cooling by the disk infrared emission. The thermal evolution of the vertical column is described by the system of radiative transfer Eqs. (1212) and (13).

This system is a set of moment equations for diffuse IR radiation derived under the Eddington approximation. We solve Eqs. (1212) and (13) numerically to find T and E after one hydrodynamical time step Δt using the implicit finite-difference scheme where T∗ is the gas temperature updated after the source step (see Eq. (A.1)), and En and ρn the radiation energy and gas volume density taken from the previous time step. In the discretized equations, we use the following convention: index i refers to the left-hand side grid interface and index i + 1 / 2 to the grid center. In the above system, all quantities are defined at the grid center i + 1 / 2 (the index omitted for simplicity) and denotes the finite-difference approximation to the diffusion term: (A.10)where The above nonlinear system is solved with the iterative Newton method. Namely, we approximate T4 as using the first two terms of the Taylor expansion series, where superscript k refers to the previous iterative step and produce a system of linear equations for E1 / 2,...,EM − 1 / 2 with a tridiagonal matrix, where M is the number of grid cells in the vertical direction. This tridiagonal system is solved with the forward and back substitution method (the Thomas algorithm). The obtained values of energies are back substituted in Eq. (A.8) to derive the temperatures and both quantities are then used for the next Newton iteration as Tk + 1 and Ek + 1 until convergence is achieved. The solution is then set to T = Tk + 1 and E = Ek + 1.

The boundary conditions for Eqs. (A.8) and (A.9) are the zero energy gradient near the mid-plane and following relation between the diffusion flux and radiation energy at the disk surface (A.11)where Tcmb = 2.73 K. This condition is derived under the assumption that the incoming and outgoing radiation is isotropic over the upper and lower hemispheres. We note that the presented method is conceptually similar to that developed for the calculation of the pre-stellar core thermal evolution in Pavlyuchenkov et al. (2015).

The next stage of the algorithm is to recover the hydrostatic equilibrium along the vertical direction using the updated temperatures. Equation (14) describing the vertical hydrostatic equilibrium can be written as where r is the radial distance to a given vertical column of gas, G the gravitational constant, and σ the gas column density measured from the disk mid-plane. We note that T is known from the previous step. Equation (A.13)accounts for the gravity of the central star and disk self-gravity in the plane-parallel limit. The boundary conditions for this system have the form where ρext is the gas volume density at the disk surface. In our models, ρext corresponds to a hydrogen number density of 103 cm-3.

The solution of the system (A.12)−(A.13)with the corresponding boundary conditions can be found using an implicit scheme similar to that used when computing the internal structure of the stars. We linearize the right-hand side of (A.12)as (A.16)which transforms the initial system into a linear system of ordinary differential equations. The implicit finite-difference approximation of this equation generates a system of algebraic linear equations with a tridiagonal matrix whose solution can be easily found using the forward and backward substitution method. The resulting solution is used to form a new approximation, and iterations over k are carried out until convergence is achieved (usually after a few iterations).

At the last step of the algorithm, we calculate the updated thermal energy e per unit area by summing up the thermal energies over the vertical grid (A.17)where a factor of 2 accounts for the fact that we adopted an equatorial symmetry with resect to the disk midplane. This value will also be used in the transport step.

Finally, we note that the transport step updates the values of Σ and e (see Eqs. (17) and (19)), which in turn affects the volumetric distributions of the gas temperature T and volume density ρ. To take these changes into account, we assume that the mass and thermal energy that are carried with the flow are evenly redistributed over the vertical cells so that the updated distributions of T and ρ can be calculated as where T, ρ, and Σ are the values of the gas temperature, its volume, and surface densities before the transport step. This completes one cycle of integration and the updated values Tn + 1, ρn + 1, and Σn + 1 are used at the next time step.

Since the method is fully implicit, it is in general very stable. However, in a very few cases the Newton procedure may not converge for some values of the time step. In this case, we divide the hydrodynamic time step into several subcycles, which usually solves the problem. In the limit of large time steps, the method yields a steady-state solution similar to that obtained by the steady-state model described in Vorobyov et al. (2014).

thumbnail Fig. A.1

Vertical temperature distributions of a protoplanetary disk calculated with our method and with the “diskstruct” package written by C. Dullemond. More specifically, the results from the “distruct-fullrt” method with the full treatment of the frequency and angular dependent radiative transfer are shown with the thick blue lines. The thick red curves represent the results from the approximate “diskstruct-vertrt” method, which adopts gray opacities for the dust thermal radiation. The temperature distributions calculated with our method are shown with the thin black lines. The presented vertical distributions are calculated for three radial positions in the disk, namely for R = 1 AU (top), R = 10 AU (middle), and R = 100 AU (bottom) for an optically thick protoplanetary disk illuminated by a star with Teff = 6000 K.

The presented algorithm was carefully tested and compared with other methods. In particular, we benchmark our steady-state solutions (in the limit of large time steps) with the results of 1+1D code “diskstruct”2 developed by C. Dullemond and discussed in Dullemond et al. (2002). The results of this comparison are shown in Fig. A.1. The vertical temperature distributions calculated with our method are very similar to the results obtained with the approximate method “vertrt”, which also adopted gray opacities for the dust thermal radiation included in the “diskstruct” package. The deviation of our solution from the results produced by the more accurate method “fullrt”, which takes into account the frequency and angular dependence of radiative transfer, is notable near the equatorial plane. However, we consider these differences to be acceptable for our simulations.

thumbnail Fig. A.2

Evolution of the gas temperature distribution at r = 10 AU in a toy model of the protoplanetary disk. The density structure is fixed. The initial temperature distribution is uniform with T = 4 K (first case) and T = 400 K (second case). The initial distributions are shown in shades of blue, the steady-state solution is shown in red. The arrows illustrate the evolution of the gas temperature profiles with time. The corresponding time in seconds is shown in the color legend.

To illustrate the time-dependent aspect of the radiative transfer model used in the thermal step, we show in Fig. A.2 the evolution of the temperature distribution at r = 10 AU where the density distribution is manually fixed. For the parameters of the considered vertical column of gas, the relaxation to a steady-state solution is achieved within a few years (about an order of magnitude faster than the dynamical time at r = 10 AU) starting from either a warmed-up disk (with an initial uniform temperature of T = 400 K) or from a cooled-down disk (with a uniform temperature of T = 4 K). This relaxation time is in agreement with the characteristic thermal time estimated in Vorobyov et al. (2014). As expected, the relaxation to the steady-state solution is faster in the upper layers of the disk where the gas volume density is lower.

thumbnail Fig. A.3

Evolution of the gas volume density distribution at r = 10 AU in a model of protoplanetary disk. The initial temperature distribution is uniform with T = 4 K (first case) and T = 400 K (second case). The initial density distributions are shown in shades of blue, the steady-state solution is shown in red. The arrows illustrate the evolution of the gas density profiles with time. The corresponding time in seconds is shown in the color legend.

While producing the temperature distributions in Fig. A.2, the gas volume density distribution was manually fixed. When the density is allowed to evolve together with the temperature, then the density distribution develops as shown in Fig. A.3.

We see that at low and high initial temperatures the disk vertical height is low and high, respectively. In the steady-state profile near z = 1.2 AU a slight change in the slope can be seen that corresponds to the change in the temperature profile.

All Tables

Table 1

Model core parameters.

In the text

All Figures

thumbnail Fig. 1

Gas surface density distributions in model A+ shown at several time instances after the onset of numerical simulations (the disk forms at t = 0.1 Myr). The insets zoom in on several fragments in the disk. The yellow contour lines delineate regions with the Toomre Q-parameter smaller than unity. The red dashed line is a radial cut through the disk used later to show the vertical volume density and temperature distributions. The scale bar is in log g cm-2.
In the text
thumbnail Fig. 2

Similar to Fig. 1, but for model A.
In the text
thumbnail Fig. 3

Number of fragments in the disk at a certain time instance as a function of time elapsed since the onset of numerical simulations in model A+ (top) and model A (bottom). The disk forms at t = 0.1 Myr.

In the text
thumbnail Fig. 4

Azimuthally averaged radial profiles of the gas temperature in model A+ (red and blue lines) and model A (black solid lines) taken at several times since onset of numerical simulations. In particular, the red lines represent the gas temperature mass-weighted over the vertical column of gas, while the blue lines are the gas temperature in the disk midplane. The dashed solid lines are the gas temperatures in model A derived for the steady-state case neglecting compressional and viscous heating. The vertical dotted lines mark the radial positions where the disk joins the infalling envelope.

In the text
thumbnail Fig. 5

Distributions of gas volume density (top) and temperature (bottom) in the (r,z) plane taken in model A+ at t = 0.25 Myr through the radial cut with a position angle of 351° (red dashed line in Fig. 1).
In the text
thumbnail Fig. 6

Temperature and gas surface density distributions in model A+. Various definitions of temperature as described in the text are shown.

In the text
thumbnail Fig. 7

One-dimensional vertical cuts taken at different positions of the disk (as indicated in the legend) and showing the vertical gas volume density (top) and temperature (bottom) distributions.

In the text
thumbnail Fig. 8

Number of fragments in the disk at a certain time instance as a function of time elapsed since the onset of numerical simulations in model B+ (top) and model B (bottom). The disk forms at t = 0.138 Myr.

In the text
thumbnail Fig. A.1

Vertical temperature distributions of a protoplanetary disk calculated with our method and with the “diskstruct” package written by C. Dullemond. More specifically, the results from the “distruct-fullrt” method with the full treatment of the frequency and angular dependent radiative transfer are shown with the thick blue lines. The thick red curves represent the results from the approximate “diskstruct-vertrt” method, which adopts gray opacities for the dust thermal radiation. The temperature distributions calculated with our method are shown with the thin black lines. The presented vertical distributions are calculated for three radial positions in the disk, namely for R = 1 AU (top), R = 10 AU (middle), and R = 100 AU (bottom) for an optically thick protoplanetary disk illuminated by a star with Teff = 6000 K.
In the text
thumbnail Fig. A.2

Evolution of the gas temperature distribution at r = 10 AU in a toy model of the protoplanetary disk. The density structure is fixed. The initial temperature distribution is uniform with T = 4 K (first case) and T = 400 K (second case). The initial distributions are shown in shades of blue, the steady-state solution is shown in red. The arrows illustrate the evolution of the gas temperature profiles with time. The corresponding time in seconds is shown in the color legend.

In the text
thumbnail Fig. A.3

Evolution of the gas volume density distribution at r = 10 AU in a model of protoplanetary disk. The initial temperature distribution is uniform with T = 4 K (first case) and T = 400 K (second case). The initial density distributions are shown in shades of blue, the steady-state solution is shown in red. The arrows illustrate the evolution of the gas density profiles with time. The corresponding time in seconds is shown in the color legend.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.