© ESO 2020

1 Introduction

Protoplanetary disks are an important ingredient of star and planet formation. They form during the gravitational contraction of rotating pre-stellar cloud cores thanks to the conservation of angular momentum of the infalling material. It is thought that most of the core material is processed by the disk before it lands on the growing protostar or leaves the system through protostellar jets and outflows. This processing includes alterations in the chemical composition and growth of submicron particles to centimeter-sized pebbles, providing building blocks for planets. These processes critically depend on the thermal balance in the disk and knowing the properties of protoplanetary disks is therefore of prime importance for our understanding of star and planet formation.

The studies of protoplanetary disks have traditionally followed two separate pathways. Global simulations follow the collapse of pre-stellar cores to the protostellar stage characterized by the formation of a star and protostellar disk (e.g., Machida et al. 2010; Joos et al. 2013; Seifried et al. 2013; Tsukamoto et al. 2015). The forming disks usually show a very complex behavior depending on the mass of the core, amount of initial rotation in the core, and the strength of magnetic fields (e.g., Bate 2018; Wurster & Bate 2019). The interaction with the environment in the form of jets, outflows, and infalling material makes the interpretation of numerical simulations a challenging task. An alternative approach is to look into the evolution of isolated, so to say, already-formed disks (e.g., Kley 1999; Boss 2002; Stamatellos et al. 2007; Mayer et al. 2007). Such an approach allows us to focus on particular aspects of disk evolution and usually permits a better numerical resolution, making this approach particularly valuable for our understanding of the subtleties of star and planet formation.

One of the important aspects of disk evolution is gravitational instability and fragmentation, which is thought to be a possible gateway for the formation of giant planets and brown dwarfs (e.g., Boss 2002; Mayer et al. 2007; Boley et al. 2010; Vorobyov 2013; Meru 2015; Nayakshin 2017; Mercer & Stamatellos 2017). Gravitational instability and fragmentation are particularly sensitive to the disk mass, but also to the thermal balance in the disk controlled largely by disk cooling, and viscous and stellar heating. Starting from Gammie (2001), it has become increasingly popular to employ the so-called β-parameterization to describe the cooling processes in the disk when studying the disk propensity to gravitational fragmentation (e.g., Rice et al. 2003; Cossins et al. 2009; Meru & Bate 2011; Boss 2017; Deng et al. 2017). In this approach the rate of disk cooling is parameterized in terms of the β-parameter, which is the product of the local cooling time t_c and local angular velocity Ω. The popularity of this approach can be explained by its simplicity, which avoids the complicated physics and numerics often involved with solving for the full energy balance equation and at the same time provides valuable insights into an important physical process: disk gravitational fragmentation.

The often used approach is to adopt a uniform β-parameter throughout the disk, but vary its value to determine its effect on disk gravitational fragmentation. There is, however, no justification that the β-parameter should be uniform throughout the disk. Moreover, it is not clear if the β-parameterization accurately describes the thermal balance in the disk. Several numerical studies address the disk propensity to fragment depending on whether the simplified β-cooling or a more sophisticated cooling–heating scheme is used (e.g., Gammie 2001; Johnson & Gammie 2003), but a systematic study of the temperature distribution in the disk on global evolutionary timescales has not been performed yet. It is therefore important to understand the risks that are involved with using the simplified β-approximation.

In this paper we consider three different approaches to describing the thermal balance in the disk. First, we consider the scheme that takes radiative cooling, stellar and viscous heating, and PdV work into account in the limit of equal dust and gas temperatures. This approach has been extensively used in one- and two-dimensional disk dynamics simulations (Johnson & Gammie 2003; Rice & Armitage 2009; Vorobyov & Basu 2010; Zhu et al. 2012) and also adapted to three-dimensional smoothed-particle simulations (Stamatellos et al. 2007). Second, we introduce a new cooling–heating scheme that allows a separate calculation of the gas and dust temperatures in protoplanetary disks. This scheme is similar in methodology (but not exactly the same) to the methods earlier presented by Pavlyuchenkov et al. (2015) and Bate & Keto (2015). We particularly search for disk regions where the gas temperature can deviate notably from that of dust. Finally, we consider the simplified β-cooling and determine the applicability of this simplified approach in describing disk evolution.

We use disk models in the thin-disk limit to study the effect of different cooling–heating schemes. A simplified disk dynamics allows us to focus on the thermal properties of the disk, and to run disk simulations for a much longer time than in full three-dimensional simulations. Nevertheless, we believe that our results regarding the applicability of the β-approximation and importance of separate dust and gas thermal evolution will remain valid in fully three-dimensional disk models. More importantly, our thermal model can find applications in simulations of low-metallicity disks, where decoupling of dust and gas temperatures is expected to be significant.

The paper is organized as follows. In Sect. 2 we describe in detail different cooling–heating schemes that we used in our disk models. In Sect. 3 we describe the disk evolution using the cooling–heating scheme with separate dust and gas temperatures. In Sect. 4 we compare the disk evolution in models with separate and similar dust and gas temperatures. In Sect. 5 we consider models with a simplified β-cooling. The main results are summarized in Sect. 6.

2 Model description

In this section, we describe the main aspects of our model regarding the gas dynamics computations, while the subsequent subsections elaborate on computations of the thermal balance in the disk. We use numerical hydrodynamics simulations in the thin-disk limit to compute the formation and global evolution of young circumstellar disks. To avoid time steps that are too small, we set a dynamically inactive sink cell in the center of our computational domain with a radius of r_sc = 5 au.

The starting point of each simulation is the gravitational collapse of a pre-stellar core. In the adopted thin-disk approximation, the core has the form of a flattened pseudo-disk, a spatial configuration that can be expected in the presence of rotation and large-scale magnetic fields (e.g., Basu 1997). As the collapse proceeds, the inner regions of the core spin up and a centrifugally balanced circumstellar disk forms when the inner infalling layers of the core hit the centrifugal barrier near the sink cell. The material that has passed to the sink before the instance of circumstellar disk formation constitutes a seed for the central star, which grows further through accretion from the circumstellar disk. The infalling core continues to land at the outer edge of the circumstellar disk until the core depletes. The infall rates on the circumstellar disk are in agreement with what can be expected from the free-fall collapse (Vorobyov 2010). Computations continue up to 0.5 Myr, thus covering the entire embedded phase and the early T Tauri phase of disk evolution.

We take into account turbulent viscosity described via the Shakura & Sunyaev α-parameterization and disk self-gravity. The forming protostar is not just a source of gravity. Its characteristics, such as the radius and photospheric luminosity, are calculated in line with the disk evolution using the stellar evolution tracks obtained with the STELLAR code (Yorke & Bodenheimer 2008). These characteristics are then used to calculate the total stellar luminosity and the radiation flux impinging the surface of the disk and contributing to its heating in models where detailed disk cooling and heating are taken into account.

The equations of mass and momentum in the thin-disk limit are $$\frac{\partial \Sigma }{\partial t}=-{\mathbf{\nabla }}_p\cdot \left(\Sigma {\mathbf{v}}_p\right),$$(1) $$\frac{\partial }{\partial t}\left(\Sigma {\mathbf{v}}_p\right)+{\left[\mathbf{\nabla }\cdot \left(\Sigma {\mathbf{v}}_p\otimes {\mathbf{v}}_p\right)\right]}_p=-{\mathbf{\nabla }}_{p}P+\Sigma {\mathbf{g}}_p+{\left(\mathbf{\nabla }\cdot \mathbf{\Pi }\right)}_p,$$(2)

where the subscripts p and p′ refer to the planar components (r, ϕ) in polar coordinates, Σ is the gas mass surface density, P is the vertically integrated gas pressure calculated via the ideal equation of state as P = (γ − 1)e, γ is the ratio of specific heats, ${\mathbf{v}}_p=v_{\text{r}}\widehat{\mathbf{r}}+v_{\varphi }\widehat{\mathbf{\varphi }}$ is the velocity in the disk plane, ${\mathbf{g}}_p=g_{\text{r}}\widehat{\mathbf{r}}+g_{\varphi }\widehat{\mathbf{\varphi }}$ is the gravitational acceleration in the disk plane (including that of the disk and the star), and ${\mathbf{\nabla }}_{\text{p}}=\widehat{\mathbf{r}}\partial /\partial r+\widehat{\mathbf{\varphi }}{r}^{-1}\partial /\partial \varphi$ is the gradient along the planar coordinates of the disk. Turbulent viscosity enters the basic equations via the viscous stress tensor Π, and we calculate the magnitude of kinematic viscosity ν using the α-parameterization with a spatially uniform α-parameter. Two limiting cases were considered: an MRI-active disk with α = 0.01 and an MRI-suppressed disk with α = 10⁻⁴.

2.1 β-cooling

In its simplest form, the energy balance equation can be described as $$\frac{\partial e}{\partial t}+{\mathbf{\nabla }}_{\text{p}}\cdot \left(e{v}_p\right)=-P\left({\mathbf{\nabla }}_p\cdot v_p\right)-\frac{e}{t_{\text{c}}}+{\left(\mathbf{\nabla }\cdot v\right)}_{p{p}^{\prime }}:{\mathbf{\Pi }}_{p{p}^{\prime }},$$(3)

where e is the internal energy of gas per surface area. The characteristic cooling time is related to the β-parameter as t_c = β∕Ω. This approach, hereafter referred to in the text as the β-cooling scheme, takes into account the advection of internal energy with the gas flow, heating and cooling through adiabatic compressionand expansion of gas flows (first term on the right-hand side), radiative cooling approximated by the β-parameter, and viscous heating (last term on the right-hand side). We note that in many studies with the β-cooling scheme, but not in our paper, viscous heating is neglected. To avoid a catastrophic overcooling of the disk we turn off the β-cooling term as soon as the gas temperature drops below a threshold value set equal to 4 K. This β-cooling approach was employed in studies of disk fragmentation in, for example, Rice et al. (2003), Meru & Bate (2011), Boss (2017), and Rice & Nayakshin (2018). In our study, the energy Eq. (3) is closed with the ideal equation of state P = e(γ − 1) for a perfect gas, for which the ratio of specific heats is set equal to a constant value of γ = 1.4.

One drawback of the above approach is that it does not take stellar irradiation into account, although this heating mechanism can be important once the star has formed. There are various modifications to the standard β-cooling scheme (see, e.g., Baehr & Klahr 2015), but in this study we adopt the form $$\frac{\partial e}{\partial t}+{\mathbf{\nabla }}_{\text{p}}\cdot \left(e{v}_p\right)=-P\left({\mathbf{\nabla }}_p\cdot v_p\right)-\frac{e-e_{\text{irr}}}{t_{\text{c}}}+{\left(\mathbf{\nabla }\cdot v\right)}_{p{p}^{\prime }}:{\mathbf{\Pi }}_{p{p}^{\prime }},$$(4)

where e_irr is the internal energy per surface area defined exclusively by stellar and background irradiation as e_irr = ΣℛT_irr∕μ, with the mean molecular weight set equal to μ = 2.33. Here, T_irr is the irradiation temperature calculated as $$T_{\text{irr}}^4=T_{\text{bg}}^4+\frac{F_{\text{irr}}(r)}{\sigma },$$(5)

where T_bg is the uniform background temperature set equal to the initial temperature of the natal cloud core and F_irr (r) is the radiation flux absorbed by the disk surface at a radial distance r from the central star. The flux is calculated as $$F_{\text{irr}}(r)=\frac{L_{\ast }}{4\pi r^2}\cos {\gamma }_{\text{irr}},$$(6)

where γ_irr is the incidence angle of radiation arriving at the disk surface (with respect to the normal) at radial distance r. The incidence angle is calculated using a flaring disk surface, as described in Vorobyov & Basu (2010). The stellar luminosity L_* is the sum of the accretion and stellar photospheric luminosities. The modified β-cooling term works as a relaxation process on a timescale t_c toward the thermal state defined by T_irr. The stronger the mismatch between e and e_irr, the faster the system strives to attain the thermal state defined by irradiation (for a fixed value of β). We note that the other terms on the right-hand side of Eq. (4) can act to push the thermal balance away from that established by stellar and background irradiation.

2.2 Similar thermal evolution of gas and dust

A more sophisticated approach to computing the thermal balance in the disk involves solving the energy equation considering the effects of viscous heating, disk radiative cooling, and stellar heating via irradiation. The equation reads $$\frac{\partial e}{\partial t}+{\mathbf{\nabla }}_{\text{p}}\cdot \left(e{v}_p\right)=-P\left({\mathbf{\nabla }}_p\cdot v_p\right)-\Lambda +\Gamma +{\left(\mathbf{\nabla }\cdot v\right)}_{p{p}^{\prime }}:{\mathbf{\Pi }}_{p{p}^{\prime }},$$(7)

where Λ and Γ are the cooling and heating rates due to dust cooling and stellar (and background) irradiation, respectively.

This approach, hereafter referred to as the thermal evolution scheme 1 (or ThES1), was employed to study disk fragmentation in, for example, Johnson & Gammie (2003), Vorobyov & Basu (2010), and Zhu et al. (2012). The difference between these studies lies in the degree of sophistication in calculating the radiative cooling and disk heating, in neglecting or taking viscous heating into account. For instance, Johnson & Gammie (2003) considered only a cooling term and neglected disk heating through stellar irradiation and viscosity. Their cooling term reads $$\Lambda =\frac{16}{3}\sigma T_{\text{mp}}^4\frac{{\tau }_{\text{R}}}{1+{\tau }_{\text{R}}^2},$$(8)

where τ_R is the mean Rosseland optical depth, T_mp is the midplane temperatureof dust (and gas), and σ is the Stefan–Boltzmann constant. Zhu et al. (2012) added irradiation heating by the central star (and also background) in the form $$\Gamma =\frac{16}{3}\sigma T_{\text{irr}}^4\frac{{\tau }_{\text{R}}}{1+{\tau }_{\text{R}}^2}.$$(9)

Vorobyov & Basu (2010) also considered viscous heating due to turbulence via α-parameterization.

The form of the cooling–heating terms may vary depending on the degree of sophistication in calculating radiative cooling of dust from the disk surface. In this work we use the expressions derived in Dong et al. (2016), $$\Lambda =\frac{8{\tau }_{\text{P}}\sigma T_{\text{mp}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(10) $$\Gamma =\frac{8{\tau }_{\text{P}}\sigma T_{\text{irr}}^4}{1+2{\tau }_{\text{P}}+\frac{3}{2}{\tau }_{\text{R}}{\tau }_{\text{P}}},$$(11)

where τ_P is the Planck optical depth. We note that the cooling and heating rates in Dong et al. (2016) were written for one side of the disk and need to be multiplied by a factor of 2. The energy Eq. (7) is closed with the ideal equation of state P = e(γ − 1), where the ratio of specific heats is set equal to a constant value of γ = 1.4. We also note that the cooling and heating terms of the form similar to Eqs. (8) and (9) are often used together with the viscous equation of Pringle (1981) for the gas surface density to compute the thermal balance in the disk (e.g., Rice et al. 2010; Kimura 2016).

2.3 Different thermal evolution of gas and dust

The thermal evolution scheme considered in Sect. 2.2 makes no difference between the gas and dust temperatures. This is a valid approximation at high densities when collisions between gas molecules and dust particles are sufficiently frequent to establish a thermal equilibrium between these two disk subsystems on timescales much shorter than the dynamical one. However, it is not clear a priori if this condition is fulfilled throughout the entire extent of a protostellar or protoplanetary disk. In the outer disk regions densities may be too low to provide strict thermal coupling between gas and dust. In this section we present a new cooling–heating scheme, referred to as the thermal evolution scheme 2 or ThES2, which is designed to lift the limitation of equal gas and dust temperatures. In ThES2 we also have a spatially and temporally varying ratio of specific heats γ, thus lifting another limitation of the β-cooling and ThES1 schemes, for which a perfect gas with constant γ was assumed.

In ThES2 we do not make a clear distinction between the cooling and heating rates, as was done with Λ and Γ in the previous sections, and we introduce the integrated rate of energy loss or gain per surface area Q_tot. The evolution equation for the gas internal energy per surface area in ThES2 reads as $$\frac{\partial e}{\partial t}+{\mathbf{\nabla }}_{\text{p}}\cdot \left(e{v}_p\right)=-P\left({\mathbf{\nabla }}_p\cdot v_p\right)-Q_{\text{tot}}+{\left(\mathbf{\nabla }\cdot v\right)}_{p{p}^{\prime }}:{\mathbf{\Pi }}_{p{p}^{\prime }},$$(12)

where Q_tot is defined as $$Q_{\text{tot}}=\left(Q_{\text{cont}}+Q_{\text{H2}}+Q_{\text{HD}}+Q_{\text{chem}}+Q_{\text{metal}}\right)2H,$$(13)

where H is the vertical scale height calculated assuming a local hydrostatic equilibrium in the gravitational field of the star and disk (see Vorobyov & Basu 2009), Q_cont is the rate of radiative energy loss (or gain) in the infrared continuum, $Q_{{\text{H}}_{\text{2}}}$ is the H₂ line cooling rate, Q_HD is the HD line cooling rate, Q_chem is the chemical cooling–heating rate (through chemical reactions), and Q_metal is the metal line cooling–heating rate. All constituents of Q_tot are volumetric cooling or heating rates and, for simplicity, they are assumed to be independent of the vertical distance from the disk midplane. This assumption allows us to convert the volumetric rates to the rates per surface area by means of vertical integration and multiplication by the disk thickness 2H. We describe how to calculate these individual cooling rates below.

The net rate of continuum cooling by energy transport from gas to radiation per unit volume is $$Q_{\text{cont}}=4\pi \left(\eta -{\chi }_{\text{a}}J\right),$$(14)

where η is the emission coefficient, J is the mean intensity, and χ_a is the absorption coefficient, given by $${\chi }_{\text{a}}=({\kappa }_{\text{P},\text{d}}+{\kappa }_{\text{P},\text{g}})\rho ,$$(15)

with the mass density ρ. We calculate the Planck mean opacities using the tables from Semenov et al. (2003) for the dust and Mayer & Duschl (2005) for the gas. We note that Semenov opacities are defined per unit gas mass assuming a dust-to-gas mass ratio of 1:100. The emission coefficient is $$\eta =\frac{\sigma \rho }{\pi }\left({\kappa }_{\text{P},\text{g}}{T}_{\text{g}}^4+{\kappa }_{\text{P},\text{d}}{T}_{\text{d}}^4\right),$$(16)

where T_g and T_d are the gas and dust temperatures, respectively. The gas temperature is determined from the ideal equation of state $P=\Sigma R{T}_{\text{g}}/\mu$, where μ is the mean molecular weight, σ is the Stefan–Boltzmann constant, and $R$ is the universal gas constant.

The dust temperature is determined in the steady-state limit by the energy balance on dust grains due to the thermal emission, absorption, and collision with gas (Omukai et al. 2010), $${\kappa }_{\text{P},\text{d}}B(T_{\text{d}})={\kappa }_{\text{P},\text{d}}J+{\Gamma }_{\text{coll}},$$(17)

where B(T) is the Planck function given by $$B(T)=\frac{\sigma }{\pi }{T}^4,$$(18)

where T is the temperature of dust or stellar irradiation (see Eq. (20) below) and Γ_coll is the heating rate of dust through collisions with gas particles. The collisional heating rate is (Hollenbach & Mckee 1979) $${\Gamma }_{\text{coll}}=4.4\times {10}^{-6}{\left(f/\rho \right)}_{\text{dust}}{n}_{\text{H}}{\left(\frac{T_{\text{g}}}{1000 \text{K}}\right)}^{1/2}\left(T_{\text{g}}-T_{\text{d}}\right),$$(19)

where ${\left(f/\rho \right)}_{\text{dust}}$ is the total volume of dust per unit gas mass and f is the mass fraction of dust grains, both taken from Pollack et al. (1994). We note that the steady-state assumption for the dust temperature allowed us to eliminate Γ_coll from the gas internal energy equation by rewriting Γ_coll in terms of the Planck function and the mean intensity.

The mean intensity used in Eqs. (14) and (17) is $$J=\frac{1}{1+x}\left(B(T_{\text{irr}})+x\frac{\eta }{{\chi }_{\text{a}}}\right),$$(20)

where x is the function that smoothly connects the optically thin and thick limits, written by Tanaka & Omukai (2014) as $$x={\tau }_{\text{P}}+\frac{3}{4}{\tau }_{\text{P}}{\tau }_{\text{R}},$$(21)

with the Planck and Rosselanck mean optical depths τ_P and τ_R, respectively.The Planck (Rosseland) mean optical depth is calculated as $${\tau }_{\text{P}(\text{or} \text{R})}=\frac{1}{2}\left({\kappa }_{\text{P}(\text{or} \text{R}),d}+{\kappa }_{\text{P}(\text{or} \text{R}),\text{g}}\right)\Sigma .$$(22)

We obtain the Rosseland mean opacities in the similar manner to the Planck mean opacities.

The H₂ and HD-linecooling rates are calculated by the following similar form: $$Q_{{\text{H}}_2(\text{HD})}={\overline{\beta }}_{\text{esc},{\text{H}}_{\text{2}}(\text{HD})}{Q}_{{\text{H}}_2(\text{HD}),\text{thin}}{\text{e}}^{-\sqrt{{\tau }_{\text{P}}{\tau }_{\text{R}}}}.$$(23)

Here $Q_{{\text{H}}_{\text{2}}(\text{HD}),\text{thin}}$ is the cooling rate in the optically thin regime given by the fitting function for H₂ from Glover (2015) and for HD from Flower et al. (2000). We take into account the line-averaged escape probability $\overline{\beta }$ to consider the effect of photon trapping in the large column density case. The values of line-averaged escape probabilities for H₂ and HD are obtained by using the fitting functions in Fukushima et al. (2018) and Eq. (A.2).

The chemical cooling and heating are the processes associated with chemical reactions. We follow the chemical evolution of eight species (H, H₂, H⁺, H⁻, D, HD, D⁺, and e) and take into account the 21 hydrogen and 6 deuterium reactions summarized in Table B.1. We consider H ionization–recombination and H₂ dissociation–formation as the chemical cooling–heating processes. The chemical cooling rate is $$Q_{\text{chem}}=\left({\varepsilon }_{\text{H}}\frac{\text{d}y({\text{H}}^{+})}{\text{d}t}-{\varepsilon }_{{\text{H}}_2}\frac{\text{d}y({\text{H}}_2)}{\text{d}t}\right)n_{\text{H}},$$(24)

where ɛ_H = 13.6 eV and ${\varepsilon }_{{\text{H}}_2}$ = 4.48 eV are the binding energies. The chemical fraction of species i is defined using the number density of species i, n(i), and that of hydrogen nuclei n_H as $$y(i)=\frac{n(i)}{n_{\text{H}}}.$$(25)

The numberdensity of hydrogen nuclei is $$n_{\text{H}}=\frac{\rho }{(1+4y_{\text{He}})m_{\text{H}}},$$(26)

where y_He is the number fraction of He relative to hydrogen nuclei and m_H is the hydrogen nuclei mass.

We consider the atomic fine-structure line emission of CII and OI as the metal line cooling Q_metal. We model CII as a two-level system and OI as a three-level system and count level populations from the statistical balance among each level. We take the level energies, the spontaneous radiative decay rates, and the collisional deexcitation rate coefficients from Hollenbach & Mckee (1989). The metal line cooling–heating rate can be divided into the line cooling–heating rates of CII and OI as follows: $$Q_{\text{metal}}=Q_{\text{CII}}+Q_{\text{OI}},$$(27) $$\begin{array}{lll}{Q}_{\text{CII}(\text{OI})}\hfill & =\hfill & y_{\text{CII}(\text{OI})}{n}_{\text{H}}\left(Z/Z_{\text{local}}\right)\hfill \\ \hfill & \hfill & \times {\sum }_{\text{ul}}h{\nu }_{\text{ul}}{\beta }_{\text{esc},\text{ul}}{A}_{\text{ul}}{f}_{\text{u}}\frac{S({\nu }_{\text{ul}})-B({\nu }_{\text{ul}};T_{\text{rad}})}{S({\nu }_{\text{ul}})}.\hfill \end{array}$$(28)

Here the chemical fractions of CII and OI are y_CII = 9.27 × 10⁻⁵ and y_OI = 3.568 × 10⁻⁴, Z∕Z_local is the metallicity relative to solar one, hν_ul is the energy difference between the upper level u and the lower level l, β_esc,ul is the line escape probability, A_ul is the spontaneous radiative decay rate, f_u is the occupancy of upper level, S(ν_ul) is the source function, and B(ν_ul;T_rad) is the Planck function. We note that the metal lines heat gas if the gas temperature is lower than the irradiation temperature. In this work, Z∕Z_local is unity. The line escape probability is $${\beta }_{\text{esc},\text{ul}}=\left(\frac{1-{\text{e}}^{-{\tau }_{\text{ul}}}}{{\tau }_{\text{ul}}}\right){\text{e}}^{-\sqrt{{\tau }_{\text{P}}{\tau }_{\text{R}}}},$$(29)

where the optical depth for line emission τ_ul is given by $${\tau }_{\text{ul}}=\frac{c^3}{8{\pi }^{3/2}{\nu }_{\text{ul}}^3}{A}_{\text{ul}}\left(\frac{g_{\text{u}}}{g_{\text{l}}}{f}_{\text{u}}-f_{\text{l}}\right)\frac{N_{\text{CII}(\text{OI})}}{v_{\text{th}}},$$(30)

where g_u,l is the statistical weight of upper level u and lower level l, N_CII(OI) is the column density of CII (OI), and v_th is the thermal velocity. The column density of CII (OI) is $$N_{\text{CII}(\text{OI})}=2H{n}_{\text{H}}{y}_{\text{CII}(\text{OI})}.$$(31)

The thermal velocity is $$v_{\text{th}}=\sqrt{\frac{2k_{\text{B}}{T}_{\text{g}}}{\mu m_{\text{H}}}},$$(32)

where k_B is the Boltzmann constant. The source function is calculated by $$S({\nu }_{\text{ul}})=\frac{2h{\nu }_{\text{ul}}^3}{c^2}{\left[\frac{g_{\text{u}}{f}_{\text{l}}}{g_{\text{l}}{f}_{\text{u}}}-1\right]}^{-1}.$$(33)

Our thermal model is based on the minimum model of Omukai et al. (2005), which includes only CII and OI line cooling (without solving C and O chemistry) and dust cooling in addition to the primordial gas thermal and chemical processes. This model can reproduce the temperature evolution calculated by more elaborate models relatively well. We also note that line cooling is only important at low densities (≤ 10⁴ cm⁻³).

2.4 Initial and boundary conditions

In this work we considered five model cores, the parameters of which are provided in Table 1. The initial radial profile of the gas surface density Σ and angular velocity Ω of the pre-stellar core has the form $$\Sigma =\frac{r_0{\Sigma }_0}{\sqrt{r^2+r_0^2}},$$(34) $$\Omega =2{\Omega }_0{\left(\frac{r_0}{r}\right)}^2\left[\sqrt{1+{\left(\frac{r}{r_0}\right)}^2}-1\right],$$(35)

where Σ₀ and Ω₀ are the angular velocity and gas surface density at the center of the core, and r₀ is the radius of the central plateau. This radial profile is typical of pre-stellar cores with a supercritical mass-to-flux ratio that are formed through ambipolar diffusion, with the specific angular momentum remaining constant during axisymmetric core collapse (Basu 1997). All pre-stellar cores are initially unstable to gravitational collapse, but differ in the amount of mass and angular momentum. In particular, model 1 is the most massive, while model 3 is the least massive. In addition, model 3 is distinguished by a factor of 2 higher initial ratio of rotational to gravitational energy. The initial gas and dust temperatures are set equal to 10 K.

The initial chemical composition of the cores in the ThES2 is as follows. We calculate the time evolution of the central density, temperature, and chemical composition of the collapsing cloud core with the one-zone treatment as in Omukai et al. (2005) until the central density reaches 10⁶ cm⁻³. The values of the chemical fractions of the eight species at that time are y(H) = 3 × 10⁻¹⁰, y(H₂) = 0.5, y(H⁺) = y(e) = 10⁻⁸, y(D) = 2 × 10⁻¹⁶, y(HD) = 3 × 10⁻⁵, and y(H⁻) = y(D⁺) = 0.

We distinguish between different cooling and heating schemes by adding the corresponding prefix. For example (ThES1)-model 1 would correspond to model 1 with the thermal evolution scheme 1. In addition, we put the letter “v” after the model number to denote the models with an increased value of the viscous α-parameter, thus simulating a fully MRI-active disk.

The inner boundary condition located at r_sc should be chosen with a certain care. If the inner boundary allows for matter to flow only in one direction from the active disk to the sink cell, then any wave-like motions near the inner boundary, such as those triggered by spiral density waves in the disk, would resultin a disproportionate flow through the sink–disk interface. As a result, an artificial depression in the gas density near the inner boundary develops over the course of time because of the lack of compensating back flow from the sink to the disk. A solution to this problem was proposed in Vorobyov et al. (2018), where a free inflow–outflow boundary condition was introduced, allowing matter to flow freely from the disk to the central sink cell and vice versa according tothe computed mass transport rate through the sink–disk interface. In particular, the mass of material ΔM_flow (always positive definite) that passes through the sink–disk interface is further split into two components, ΔM_* and ΔM_s.c., which are used to update the gas surface density in the sink cell Σ_s.c. and the stellar mass M_* according tothe following algorithm: $$\begin{array}{l}\text{if}{\Sigma }_{\text{s}.\text{c}.}^n<{\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_r(r_{\text{s}.\text{c}.})<0\text{then}\\ {\Sigma }_{\text{s}.\text{c}.}^{n+1}={\Sigma }_{\text{s}.\text{c}.}^n+\Delta M_{\text{s}.\text{c}.}/S_{\text{s}.\text{c}.}\\ M_{\ast }^{n+1}=M_{\ast }^n+\Delta M_{\ast }\\ \text{if}{\Sigma }_{\text{s}.\text{c}.}^n<{\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_r(r_{\text{s}.\text{c}.})\ge 0\text{then}\\ {\Sigma }_{\text{s}.\text{c}.}^{n+1}={\Sigma }_{\text{s}.\text{c}.}^n-\Delta M_{\text{flow}}/S_{\text{s}.\text{c}.}\\ M_{\ast }^{n+1}=M_{\ast }^n\\ \text{if}{\Sigma }_{\text{s}.\text{c}.}^n\ge {\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_r(r_{\text{s}.\text{c}.})<0\text{then}\\ {\Sigma }_{\text{s}.\text{c}.}^{n+1}={\Sigma }_{\text{s}.\text{c}.}^n\\ M_{\ast }^{n+1}=M_{\ast }^n+\Delta M_{\text{flow}}\\ \text{if}{\Sigma }_{\text{s}.\text{c}.}^n\ge {\overline{\Sigma }}_{\text{in}.\text{disk}}^n\text{and}{v}_r(r_{\text{s}.\text{c}.})\ge 0\text{then}\\ {\Sigma }_{\text{s}.\text{c}.}^{n+1}={\Sigma }_{\text{s}.\text{c}.}^n-\Delta M_{\text{flow}}/S_{\text{s}.\text{c}.}\\ M_{\ast }^{n+1}=M_{\ast }^n.\end{array}$$

Here ${\overline{\Sigma }}_{\text{in}.\text{disk}}$ is the averaged surface density of gas in the inner active disk (the averaging is usually done over one au immediately adjacent to the sink cell), S_s.c. is the surface area of the sink cell, and v_r(r_s.c.) is the radial component of velocity at the sink–disk interface. We note that v_r (r_s.c.) < 0 when the gas flows from the active disk to the sink cell and v_r(r_s.c.) > 0 in the opposite case. The superscripts n and n + 1 denote the current and the updated (next time step) quantities. The exact partition between ΔM_* and ΔM_s.c. is usually set to 95%:5%, meaning that most of the mass lands directly on the star and only a small fraction is retained by the sink. This corresponds to fast mass transport through the sink. The effect of the ΔM_* : ΔM_s.c. partition onthe disk evolution is studied in Vorobyov et al. (2019). The calculated values of ${\Sigma }_{\text{s}.\text{c}.}^{n+1}$ are used at the next time step as the inner boundary values for the gas surface density. The radial velocity and internal energy at the inner boundary are determined from the zero gradient condition, while the azimuthal velocity is extrapolated from the active disk to the sink cell assuming a Keplerian rotation.

The rate may be both negative, meaning the flow of mass from the disk to the sink, and positive, meaning the opposite flow from the sink to the disk. The mass transport rate through the sink–disk interface is also used to calculate the net mass of gas in the sink andin the star (for details, see Kadam et al. 2019).

The known gas mass in the sink cell is then used as the inner boundary values for the surface density in the disk. The radial velocity and internal energy at the inner boundary are determined from the zero gradient condition, while the azimuthal velocity is extrapolated from the active disk to the sink cell assuming a Keplerian rotation. These inflow–otuflow boundary conditions enable a smooth transition of the surface density and angular momentum between the inner active disk and the sink cell, preventing (or greatly reducing) the formation of an artificial drop in the surface density near the inner boundary. Finally, we note that the outer boundary condition is set to a standard free outflow, allowing material to flow out of the computational domain, but not allowing any material to flow in.

2.5 Solution procedure

The continuity and momentum Eqs. (1) and (2) and also the energy Eqs. (3), (7), and (12), depending on the adopted cooling–heating scheme, are solved on the polar grid (r, ϕ) using the operator-split solution procedure similar in methodology to the ZEUS-2D code (Stone & Norman 1992). The computational domain extends from the sink cell boundary at r_sc = 5 au to the initial cloud core radius at r_out (see Table 1). The star (once formed) is located at the coordinate origin and the stellar motion in response to the disk potential is not taken into account in this study. The adopted resolution is 512 × 512 grid cells, which on the logarithmically spaced grid corresponds to a spatial resolution of 0.1 au at a radial distance of 7 au and 1.0 au at 70 au. To correctly simulate disk fragmentation, the local Jeans length must be resolved by at least four numerical cells (Truelove et al. 1998). In the thin-disk limit, the Jeans length can be expressed as (Vorobyov 2013) $$R_{\text{J}}=\frac{c_{\text{s}}^2}{\pi G\Sigma },$$(36)

where c_s is the sound speed and G is the gravitational constant. Fragments usually condense out of the densest sections of spiral arms at a typical distance of 100 au and then either migrate inward or scatter outward. The typical surface densities and temperatures in spiral arms do not exceed 100 g cm⁻² and 100 K. Adopting these values, the corresponding Jeans length is R_J ≈ 20 AU. The numerical resolution at 100 au is 1.4 au, thus fulfilling the Truelove criterion.

The solution is split into the transport and source steps. In the transport step, the update of hydrodynamic quantities due to advection is done using the third-order piecewise parabolic interpolation scheme of Colella & Woodard (1984). In the source step, the update of hydrodynamic quantities due to gravity, turbulent viscosity, cooling, and heating is performed. The gravitational potential of the matter in the computational domain is found by solving for the Poisson integral (Binney & Tremaine 1987): $$\begin{array}{l}\Phi (r,\varphi )\\ =-G{\displaystyle {\int }_{{\text{r}}_{\text{sc}}}^{{\text{r}}_{\text{out}}}{r}^{\prime }}\text{d}{r}^{\prime }{\displaystyle {\int }_0^{2\pi }\frac{\Sigma (r^{\prime },{\varphi }^{\prime })\text{d}{\varphi }^{\prime }}{\sqrt{r^{\prime }{}^2+r^2-2r{r}^{\prime }\cos ({\varphi }^{\prime }-\varphi )}}}.\end{array}$$(37)

We note that we do not introduce an explicit smoothing length when calculating the integral (Eq. (37)), as was advocated in Müeller et al. (2012) because our method for calculating the integral already includes an implicit smoother set equal to the size of the grid cell in which the potential is calculated (see Eq. (2-206) in Binney & Tremaine 1987). Since the size of the cell and the disk scale height are both linearly proportional to radial distance in our model, our implicit smoothing length is also linearly proportional to the disk scale height, in agreement with Müeller et al. (2012), but the coefficient of proportionality may be different.

We use an explicit integrator to compute the viscous force and heating (the last terms on the right-hand side of the momentum and internal energy equations). This is found to be adequate as long as the α-parameter does not greatly exceed 0.01. The update of the internal energy per surface area in the β-cooling scheme is done using an analytic solution, while in ThES1 and ThES2 the update due to cooling and heating is done implicitly using the Newton–Raphson method of root finding, complemented by the bisection method where the Newton–Raphson iterations fail to converge. The implicit solution is applied to avoid time steps that are too small that may emerge in regions of fast heating or cooling. A small amount of artificial viscosity is added to smooth out the shocks, which may occur in the gas flow, but the associated torques are much smaller than those due to turbulent viscosity.

We solve non-equilibrium kinetic equations for H, H₂, H⁺, D, HD, D⁺, and e, while the H⁻ fraction is calculated from the equilibrium of reactions 3, 4, 11, 12, 15, and 16 in Table B.1. The method of calculating the rate coefficients of the reverse reactions using the rate coefficients of the forward reactions (summarized in Table B.1) is explained in Appendix C in Matsukoba et al. (2019). We assume that helium is always neutral and its fractional abundance is y_He = 8.333 × 10⁻².

We also assume that our species are collisionally coupled with gas, which eliminates the need for solving separate equations of motion for each species. The remaining continuity equation for the surface density (Σ_i) of each of the species is written as $$\frac{\partial {\Sigma }_i}{\partial t}+{\mathbf{\nabla }}_p\cdot \left({\Sigma }_i{\mathbf{v}}_p\right)=k_{j,k}{\Sigma }_j{\Sigma }_k-k_{k,i}{\Sigma }_k{\Sigma }_i,$$(38)

where the right-hand terms are the sources and sinks due to chemical reactions. The set of Eqs. (38) is solved in two steps. First, Σ_i are updated by solving implicitly the set of non-equilibrium kinetic equations taking chemical reactions into account. This step is performed between the source and transport steps of the hydrodynamic part. Then the chemical species are advected with the gas flow using the same third-order-accurate scheme of Colella & Woodard (1984).

3 Disk evolution in ThES2

We start with describing the disk evolution in the framework of the most elaborate thermal evolution scheme ThES2 with separate gasand dust temperatures. Figure 1 presents the gas surface density distribution in the inner 2000 × 2000 au² box for the five considered models. The most massive model (in terms of the pre-stellar core) is shown in the top row, while the least massive model with different α-values (10⁻² and 10⁻⁴) is shown in the bottom two rows. The intermediate-mass model is shown in the second and third rows, also for the two values of α-parameter.

Clearly, the mass of the pre-stellar core determines the properties of the disk that form as a result of gravitational collapse. In the most massive (ThES2)-model 1v the disk is strongly fragmented during the considered evolution period (up to 0.41 Myr), while in the least massive (ThES2)-model 3v and (ThES2)-model 3 the disk only shows signatures of fragmentation in the very early stages (< 0.1 Myr) and becomes virtually axisymmetric in the later evolution. Low turbulent viscosity in (TheS2)-model2 and (ThES2)-model3 helps gravitational instability last longer, in agreement with the recent findings of Rice & Nayakshin (2018). This occurs because viscosity acts to smooth out local inhomogeneities, and also reduces the net disk mass due to an elevated mass transport. The disks in the low-viscosity models are also more compact due to the lack of viscous spreading. An increased rate of pre-stellar core rotation, as indicated by a higher β_c -value in (ThES2)-model 3v and (ThES2)-model 3, does not offset the effect of a decreased initial core mass. These two models show weaker and less persistent signatures of gravitational instability and fragmentation. Although higher β_c models can form more extended disks (gravitational instability and fragmentation are strongest at large distances), the higher M_core models formmore massive disks and this factor appears to be decisive for the development and sustainability of gravitational instability and fragmentation.

Figure 2 presents the spatial distribution of dust temperature in the inner 2000 × 2000 au² for the five considered models. Overall, the higher mass models are warmer than their lower mass counterparts, which can be explained by a higher stellar luminosity feedback in the models that form from more massive cores. The terminal stellar masses in (ThES2)-model 1v and (ThES2)-model 2v are 0.67 M_⊙ and 0.43 M_⊙, respectively. Stars with such masses have photospheric luminosities of 3.3 L_⊙ and 1.8 L_⊙, according to the adopted stellar evolution tracks from Yorke & Bodenheimer (2008). The low-mass (ThES2)-model 3v has the terminal stellar mass of 0.18 M_⊙ and its photospheric luminosity is only 0.3 L_⊙. A similar trend is found for the accretion luminosities: higher mass models have higher accretion luminosities thanks to higher accretion rates driven by more massive (and more gravitationally unstable) disks. Viscous models also have higher disk temperatures, which can be explained by additional viscous heating. Finally, we note that the gaseous clumps formed via disk fragmentation are distinguished by higher dust temperatures compared to the immediate disk environment. This temperatureincrease is caused by compressional heating (PdV work) of gravitationally bound and optically thick clumps.

One interesting feature that can be noted in Fig. 2 is a slight temperature increase at the outer edge of the disk. This effect is most pronounced in Fig. 3, however, which shows the gas temperature distribution in the inner 2000 × 2000 au² box for the five considered models. Clearly, the gas temperature distribution is strongly non-monotonic: the gas temperaturegenerally declines with radius, but there is a high-temperature rim in the disk outer regions where the gas temperaturecan exceed 100 K. We note that the gas temperature in the immediate surroundings is just a few tens of Kelvin.

To better illustrate the origin of the jump in the gas temperature distribution, we plot in Fig. 4 the gas velocity field superimposed on the gas surface density distribution in (ThES2)-model 2v at t = 0.16 Myr. The black contour line defines the disk regions where the gas surface density is equal to 0.1 g cm⁻², a value below which protoplanetary disks usually have a sharp outer edge (see Fig. 8 in Andrews et al. 2009). Clearly, the jump in the gas temperature occurs near the disk outer edge, where the gas surface density is low and where the infalling matter from the envelope meets the rotating disk. The converging gas flows produce additional compressional heating to the gas component, but the low surface densities of gas and dust prevent the gas from quickly attaining thermal equilibrium with the dust through mutual collisions. As a result, the gas temperature decouples from that of dust. This means that the gas temperature jump is expected to be most pronounced in the embedded stages of disk evolution, which seems to be the case in Fig. 3. The strength of the gas temperature jump diminishes with time in the intermediate- and low-mass models 2 and 3, for which the embedded phase ends at t = 0.19 Myr and t = 0.13 Myr, respectively (the end of the embedded phase is set to the time instance when less than 5% of the initial pre-stellar core mass still resides in the infalling envelope). The high-mass model 1 remains in the embedded phase for the entire duration of our simulations.

The reason for decoupled gas and dust temperatures can be understood from the following analysis. Assuming that radiative cooling of dust is balanced by collisional heating with gas, Eq. (17) can be expressed as $$T_{\text{d}}\simeq 120 \text{K}{\left(\frac{T_{\text{g}}}{100\text{K}} \right)}^{0.3}{\left(\frac{n_{\text{g}}}{{10}^{10} {\text{cm}}^{-3}}\right)}^{0.2},$$(39)

where n_g is the number density of gas. By setting T_g = T_d we define the threshold temperature above which gas and dust thermally decouple from each other. This threshold temperature can be written as $$T_{\text{crit}}\simeq 130 \text{K}{\left(\frac{n_{\text{g}}}{{10}^{10} {\text{cm}}^{-3}}\right)}^{0.3}.$$(40)

To illustrate the effect of threshold temperature, we take (ThES2)-model 2v at t = 0.16 Myr and plot T_crit as a function of radial distance in Fig. 5 (thick black line). We used the azimuthally averaged gas number density when calculating T_crit. The red thin and thick solid lines show the minimum and maximum azimuthal variations in the gas temperature, respectively, while the red thin and thick dashed lines present the corresponding variations for the dust temperature. We note that gas and dust temperatures coincide inside 200 au where both temperatures are lower than the threshold value. This is a thermally coupled region of the disk. The variations in gas and dust temperatures begin to deviate from each other beyond 200 au. In particular, the gas temperature becomes systematically higher than the threshold temperature T_crit, meaning that the disk is now in a thermally decoupled state.

To elaborate further on the cause for the gas temperature jump, we considered a test model in which weartificially increased the rate of collisional energy exchange between dust and gas (Γ_coll) by a factor of 50. This exercise mimics an increase in the density of both material species without affecting the integrity of the disk. If the temperature jump is due to slow exchange of energy between compressionally heated gas and radiatively cooled dust, then the jump should diminish as we increase Γ_coll. Figure 6 demonstrates that this is indeed the case. The top row presents the gas surface density, and gas and dust temperatures for the standard (ThES2)-model 2v at t = 0.16 Myr. The middle row shows the same quantities in a test model with Γ_coll increased artificially by a factor of 50. Clearly, the gas and dust temperatures in this test case are similar and the gas temperaturejump near the outer disk edge is greatly reduced. The bottom row presents the resulting distributions for (ThES2)-model 2 with a reduced rate of viscous heating. As can be seen, the rate of viscous heating does not affect notably the strength of the gas temperature jump. We conclude that gas-to-dust energy exchange defined by Eq. (19) is the most important mechanism to capture the effect of temperature decoupling. It sets the dust temperature through Eq. (17), and the resulting dust temperature enters the dominant Q_cont term in Eq. (13). The second in importance is the Q_metal term, but its effect is notable only near the disk outer edge. We note that ThES2 can be applied to a wide range of metallicities, and at lower metallicities the other terms in Eq. (13) can become important.

What could be the consequences of decoupling between the gas and dust temperatures? We show in Sect. 4 that this decoupling does not have a notable effect on the disk structure and propensity to gravitational instability and fragmentation. However, an increase in the gas temperature near the disk outer edge may have important consequences for the chemical processing of gas that flows in from the envelope. An increase in the gas temperature to more than 100 K could launch gas phase reactions that are expected to be dormant in these otherwise cold disk outer regions. For instance, Oya et al. (2016) inferred a local increase in the gas kinetic temperature in the disk outer regions of IRAS 16293-2422 based on the peculiar chemical composition, and explained this feature by a possible shock heating at the disk–envelope interface. The observational detection of a gas temperature jump is however not unambiguous, as more recent observations of IRAS 16293-2422 revealed no such structures (van’t Hoff et al. 2020). Our numerical simulations also suggest that these features are not omnipresent and their occurrence depends on the disk evolution stage.

Decoupling of gas and dust temperatures may also affect the growth rate of small (sub)micron-sized dust particles that flow in with gas from the envelope. If volatile species become oversaturated in the warm gas environment near the disk outer edge, this may facilitate the growth of icy mantles on cold dust particles (a similar effect can be observed in a Turkish bath when water vapor condenses on cold objects that are brought to the bath).

We note that the temperature decoupling between gas and dust may not only be limited to the outer disk regions. Recent studies havealready demonstrated the importance of radiative disk properties on the formation and position of gaps, spirals, and snow lines in protoplanetary disks (Zhang & Zhu 2020; Ziampras et al. 2020). The formation of gaps and rings in the dust density distribution can also lead to a reduced rate of energy transfer between gas and dust in the regions of depressed density (i.e., gaps), possibly resulting in temperature decoupling. This may have important consequences for the disk mass estimates which sensitively depend on the assumed disk temperature. We plan to explore this effect in follow-up studies.

Finally, in Fig. 7 we make a detailed comparison of the azimuthally averaged radial gas, dust, and irradiation temperatureprofiles in all five models considered. In particular, the red and blue curves present the gas and dust temperatures, while the green dotted curve is the temperature of stellar and background irradiation. Let us first consider theblack line which illustrates the maximum deviation of the gas temperature from that of dust (see the right-hand axis). Clearly, this deviation can reach hundreds of Kelvin in the massive and intermediate-mass models 1 and 2, especially in the early stages of disk evolution. The deviation peaks in the disk outer regions in the vicinity of the disk outer edge (marked with the vertical dashed lines) and diminishes in the inner parts of the disk.

We also note that the gas and dust temperatures show considerable azimuthal variations, as illustrated by the shaded areas: blue for the dust temperature variations and pink for the gas temperature variations. These variations reflect the underlying non-axisymmetric distribution of gas in the disk and the circumdisk environment (see Fig. 1). In the inner disk regions variations in the two temperatures are similar (which is why only the blue shaded area is visible). In the outer disk regions, however, the variation in the gas temperature greatly exceeds that of dust. Interestingly, the gas temperature in the regions beyond 1000 au drops systematically below that of dust (which is 10 K, set by the background irradiation). This is an inverse effect compared to that found for the disk outer edge where the gas temperature exceeds that of dust. In the regions beyond 1000 au the only notable heating mechanism is the background irradiation, which directly sets the dust temperature. Extremely low gas densities, however, prevent dust and gas temperatures from equalizing through mutual collisions, leading to progressive decoupling between the two temperatures.

Fig. 1 Gas surface density distributions in the five models considered. Each row presents a specific model, as indicated, and each column corresponds to a specific time starting from disk formation. The scale bar is in log g cm−2.

Fig. 2 Similar to Fig. 1, but for the dust temperature. The scale bar is in log K.

Fig. 3 Similar to Fig. 1, but for the gas temperature. The scale bar is in log K.

Fig. 4 Gas velocity field superimposed on the gas temperature distribution in (ThES2)-model 2v. The black contour outlines a gas surface density of 0.1 g cm−2.

Fig. 5 Decoupling of gas and dust temperatures in the disk and envelope in (ThES2)-model 2v at t = 0.16 Myr. Shown are the threshold temperature (Tcrit, black thick line) above and below which the gas and dust temperatures are thermally decoupled and coupled, respectively. The red thin and thick lines present the minimum and maximum azimuthal variations of the gas temperature, respectively,while the blue dashed thin and thick lines show the corresponding quantities for the dust temperature. The brown circleindicates the position of the disk outer edge.

Fig. 6 Comparison of the gas surface density, and gas and dust temperatures in model 2. Top and bottom rows: (ThES2)-model 2v and (ThES2)-model 2, middle row: results for a test model with the rate of collisional heat exchange between dust and gas Γcoll increased artificially by a factor of 50.

Fig. 7 Azimuthally averaged radial distributions of temperatures for gas (red solid line), dust (blue dash-dotted line), and radiation (green dotted line) in five models considered. The shaded areas indicate the range of azimuthal variations of gas (pink) and dust (blue) temperatures at each radius. The black line is the difference between the maximum temperatures of gas and dust (right-hand axis). The gray dashed line shows the radius of the disk outer edge. The arrangement of panels is the same as in Fig. 1.

Fig. 8 Comparison of gas surface density distributions in model 2v with ThES2 (separate dust and gas temperatures) and ThES1 (equal gas and dust temperatures). The scale bar is in log g cm−2.

4 Comparison of disk evolution in ThES1 and ThES2

In this section, we compare the disk evolution for two different thermal schemes ThES1 and ThES2. Our motivation is to discover whether the disk evolution with separate dust and gas temperatures (ThES2) can be notably different from the disk evolution with equal dust and gas temperatures (ThES1). For this purpose, we have chosen model 2v with the α-value set equal to 10⁻². Figure 8 presents the gas surface density distributions in the inner 2000 × 2000 au² box for (ThES2)-model 2v (top row) and (ThES1)-model v2 (bottom row). Somewhat surprisingly, the overall evolution is similar whether we consider ThES2 or ThES1. The disks in both models are gravitationally unstable and are prone to fragmentation in the initial stages of evolution.

To make a more quantitative analysis, we estimated the strength of gravitational instability by calculating the global Fourier amplitudes defined as $$C_{\text{m}}(t)=\frac{1}{M_{\text{d}}}\left|{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{r_{\text{sc}}}^{R_{\text{d}}}\Sigma }}(r,\varphi ,t)e^{im\varphi }r\text{d}r\text{d}\varphi \right|,$$(41)

where M_d is the disk mass, R_d is the disk’s outer radius set for simplicity to 500 au, and m is the ordinal number of the spiral mode. When the disk surface density is axisymmetric, the amplitudes of all modes that are not equal to zero vanish. When, for example, C_m(t) = 0.1, the perturbation amplitude of spiral density waves in the disk is 10% that of the underlying axisymmetric density distribution. The resulting Fourier amplitudes for (ThES1)-model 2v and (ThES2)-model 2v are shown in Fig. 9. It appears that (ThES2)-model 2v with distinct dust and gas temperatures is slightly more gravitationally unstable, but the difference in the Fourier amplitudes is insignificant.

Furthermore, we calculated the number of fragments in the disk at a given time instance formed through gravitational fragmentation using the fragment-tracking algorithm described in Vorobyov (2013). The results are presented in Fig. 10 for (ThES2)-model 2v (top panel) and (ThES1)-model 2v. An increase in the number of fragments shows recent fragmentation, and a decrease shows either recent tidal destruction or accretion of the fragments on the star. Again, the model with distinct dust and gas temperatures appears to form more fragments, but the model with equal dust and gas temperatures appears to sustain disk fragmentation for a longer time. To summarize, there are slight quantitative differences in the disk evolution for different thermal evolution schemes, but they do not cause qualitative changes in disk dynamics, such as suppression of disk fragmentation or disk stabilization against gravitational instability.

To further explore the difference in the evolution of model 2v with distinct thermal evolution schemes, we show in Fig. 11 the gas and dust spatial temperature distributions. More specifically, the first and second rows present the gas and dust temperature distributions in (ThES2)-model 2v, respectively, while the third row presents the temperature distribution (same for gas and dust) in (ThES1)-model 2v. In addition, in Fig. 12 we compare the azimuthally averaged gas and dust temperatures in the two considered models at the same evolutionary times as in Fig. 11. The solid lines show the gas and dust temperatures for different thermal schemes, while the dashed lines show the temperature of stellar and background irradiation.

In the very early disk evolution, the gas and dust temperatures in (TheS1)-model 2v are systematically higher than in (ThES2)-model v2, but this difference vanishes in the later evolution. The azimuthally averaged gas and dust temperature profiles become almost indistinguishable, except for a well-defined jump in the gas temperature in the outer disk regions of (ThES2)-model 2v (separate gas and dust temperatures). There is also a notable positive deviation of the gas and dust temperatures from that of stellar irradiation in the inner disk regions, which occurs due to additional viscous heating of the disk. The gas and dust temperatures are similar in the bulk of the disk because the collisional exchange of energy between gas and dust is efficient in equalizing the corresponding temperatures. Only in the outer disk regions is this trend broken, thanks to decreased gas and dust densities and increased rates of compressional heating near the disk outer edge where the inflowing envelope lands at the disk (see Figs. 4 and 5). As the compressed gas moves closer to the star, it quickly cools through increased collisions with dust and dust radiative cooling. The spatially limited extent of the disk regions where gas and dust temperatures decouple from each other may explain why the disk evolution weakly depends on the considered thermal evolution schemes.

Fig. 9 Comparison of the global Fourier amplitudes in model 2v with ThES1 (red lines) and ThES2 (black lines). The global amplitudes for four modes (m = 1, 2, 3, and 4) are shown in the four panels. The time is counted from the beginning of the core collapse.

Fig. 10 Comparison of the number of fragments in the disk in model 2v with ThES2 (top panel) and ThES1 (bottom panel). The time is counted from the instance of disk formation.

Fig. 11 Comparison of gas and dust temperature distributions in model v2 with ThES1 and ThES2. First and second rows: gas and dust temperature distributions in ThES2. Third row: temperature distribution (the same for gas and dust) in ThES1. Bottom row: comparison of the azimuthally averaged temperatures in the two considered thermal evolution schemes. The solid lines show the temperatures of gas and dust, while the dashed lines provide the temperatures of stellar irradiation. The scale bar is in log Kelvin.

Fig. 12 Comparison of the azimuthally averaged gas and dust temperature distributions in model v2 with ThES1 and ThES2. Thesolid lines show the temperatures of gas and dust, while the dashed lines provide the temperatures of stellar irradiation. The gas and dust temperatures in the (ThES2)-model v2 coincide everywhere except for the outer regions.

5 Comparison with the β-cooling scheme

In this section, we compare the disk evolution in model 2v using two thermal evolution schemes that are opposite in complexity: the most sophisticated ThES2 and the most simplified β-cooling. Our purpose is to determine whether β-cooling can be used as a valid substitute for the more sophisticated thermal evolution schemes. We considered several values of the β-parameter and distinguish between the β-models by adding a prefix beta to the model. For instance, (beta=3)-model 2v would correspond to model 2v with the β-cooling scheme and β-value equal to 3.0. In addition, we distinguish between the β-models with stellar and external irradiation by adding the suffix “Ir” to the β-value as in (beta=3Ir)-model 2v. We start by considering the case without irradiation and continue with the β-models taking irradiation into account.

5.1 The case without irradiation

Figure 13 presents the gas surface density distributions in the inner 2000 × 2000 au² region for (ThES2)-model 2v (top row) and three models 2v with different β-values: (beta=3)-model 2v (second row), (beta=10)-model 2v (third row), and (beta=30)-model 2v (bottom row). Clearly, the disk evolution in model 2v with ThES2 is notably different from that obtained with the β-cooling scheme. Theβ = 3 model carries some resemblance to the model with ThES2, but it is clearly more prone to gravitational fragmentation. The other two models with β = 10 and β = 30 are conspicuously different from the model with ThES2. The disk in these β-cooling models is much more extended initially and has a flocculent structure which is not typical of circumstellar disks.

The high-density circumstellar structure in models with β-cooling is not a truecentrifugally balanced disk with a near-Keplerian rotation, but rather a pseudo-disk with a significant deviation from circular motion. Figure 14 presents the gas velocity fields superimposed on the gas surface density distributions in (ThES2)-model 2v (top row), (beta=3)-model 2v (middle row), and (beta=10)-model 2v (bottom row) in the early stages of evolution. The red circles outline the disk regions within which the relative deviation of the azimuthally averaged angular velocity v_ϕ from the Keplerian rotation is less than 10% (in most of this inner region it does not exceed 1–2%). When calculating the Keplerian velocity we also took into account the contribution from the enclosed gaseous and dusty material. The green contour lines outline the radial extent beyond which the gas surface density drops below 0.1 g cm⁻².

In the case of ThES2, the radial extent within which the gas surface density is greater than 0.1 g cm⁻² agrees quite closely with the regions within which the deviation from the Keplerian rotation is smaller than 10%. Some mismatch is seen at t = 0.16 kyr to the south, but this is caused by a pronounced lopsidedness of the disk at this time instance. This means that Σ = 0.1 g cm⁻² may be regarded as the disk outer edge, as was already noted for protoplanetary disks in Ophiuchus by Andrews et al. (2009). When we turn to the β-cooling models, however, the mismatch between the dense regions with Σ > 0.1 g cm⁻² and regions with nearly Keplerian rotation becomes much more pronounced. This means that the dense structure outlined by the green contour lines is in fact a pseudo-disk with a significant deviation (tens of percent) from the Keplerian motion.

Figure 15 presents a comparison of the number of fragments in the disk of (ThES2)-model 2v (top panel), (beta=3)-model 2v (middle panel), and (beta=10)-model 2v (bottom panel). The disk in the β =30 model did not fragment. Clearly, the β = 3 model produces too many fragments and the β = 10 model too few fragments compared to the ThES2 model. The number of fragments decreases in the β-models with increasing β-value, as was also found in other studies of disk fragmentation using the β-cooling scheme (e.g., Meru & Bate 2011). The general trend of decreasing strength of gravitational instability with increasing β-value is expected: slower cooling leads to warmer disks and reduced gravitational instability. Overall, neither of the considered simplified β-cooling models can reproduce the strength of gravitational instability and fragmentation found in models with a more sophisticated thermal evolution scheme.

Finally, in the first and second rows of Fig. 16 we present the spatial distribution of gas temperatures in (ThES2)-model 2v and (ThES1)-model 2v, respectively. The other three rows show the gas temperatures in the β-cooling models. Each column corresponds to a specific age of the disk. The gas temperatures in the β-models are strikingly different from those of gas and dust in the ThES2 model. The inner disk regions in the former models are often colder than the periphery, which is a direct consequence of decreasing cooling time with decreasing radial distance for a spatially constant β-parameter. This trend is corroborated in Fig. 17, which shows the azimuthally averaged gas temperature profiles for the same models and at the same evolutionary times as in Fig. 16. The mismatch between the β-models and more sophisticated thermal evolution schemes is significant. We note that we had to impose an absolute lower limit on the gas temperature (4 K) in the β-models to avoid overcooling in the inner disk regions.

It may appear that we can achieve a better agreement with the ThES2 model by choosing the right β-value. This is unlikely because the characteristic cooling time t_c can be highly variable both in time and space, meaning that β is also variable. We estimated the cooling time in (ThES1)-model 2 as t_c = e∕Λ and confirmed that β = t_cΩ varies by more than an order of magnitude both radially and azimuthally within the disk extent, as Fig. 18 demonstrates.

Fig. 13 Comparison of gas surface density distributions in model v2 with ThES2 and β-cooling. First row:(ThES2)-model 2v, second row: (beta=3)-model 2v, third row: (beta=10)-model 2v, and bottom row: (beta=30)-model 2v. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.

Fig. 14 Gas velocity field superimposed on the gas surface density distribitions in (ThES2)-model 2v (top row), (beta=3)-model 2v (middle row), and (beta=10)-model 2v (bottom row) at two evolutionary times, as indicated in each panel. The red circles outline the radial distance beyond which the azimuthally averaged angular velocity deviates from the Keplerian rotation by more than 10%. The green contours outline the regions where the gas surface density drops to 0.1 g cm−1. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.

Fig. 15 Comparison of the number of fragments at a given time instance in the disks of (ThES2)-model 2v (top panel), (beta=3)-model 2v (middle panel), and (beta=10)-model 2v (bottom panel). The time is counted from the instance of disk formation.

5.2 The case with irradiation

In this section, we present the results for the β-models that take stellar and background irradiation into account. Compared to the previous section, we dropped the model with β =30, because the models with increasingly higher β-values demonstrated progressively worse agreement with the ThES2 thermal scheme. Instead, we considered lower values of β, which showed better agreement.

Figure 19 shows the gas surface density distributions in the inner 2000 × 2000 au² box for (ThES2)-model 2v (top row) and three β-models with the β-values set equal to 0.5, 3, and 10. The β-models have the same initial parameters of pre-stellar cores and the same value of the viscous α-parameter (10⁻²) as in (ThES2)-model 2v. A comparison of Figs. 13 and 19 shows that the β-models with irradiation can better reproduce the disk evolution obtained in the ThES2 thermal scheme. In particular, models with β =3 and especially with β = 0.5 possess the disk structure that is similar to what was obtained in (ThES2)-model 2v. As the red circles demonstrate, the near-Keplerian rotation is established in the β = 0.5 model throughout most of the disk extent and at all considered evolutionary times, as in (ThES2)-model 2v. Nevertheless, some differences can still be found, for example, a more compact disk in the β-models at later evolutionary stages and a more diffused disk in the β = 3.0 model in theearly evolution. The β = 10 model (and higher-β models) cannot reproduce the disk structure obtained with the most sophisticated thermal evolution scheme, whether or not we take irradiation into account.

To corroborate our conclusions, we show in Fig. 20 the number of fragments formed via gravitational fragmentation in the disk of (ThES2)-model 2v and in the disk of β =0.5 and β =1.0 models with stellar and background irradiation. We note that models with higher values of β did not show disk fragmentation, including the β = 3.0 model considered previously. Regarding the number of fragments, the β =0.5 model can best reproduce disk fragmentation in (ThES2)-model 2v. The β = 1.0 model forms too few fragments compared to (ThES2)-model 2v. Nevertheless, in both β-models disk fragmentation starts earlier and ends later than in the model with the most sophisticated thermal evolution scheme.

Finally, Fig. 21 presents the azimuthally averaged radial profiles of gas temperature in (ThES2)-model 2v and in three β-models with irradiation. Clearly, the β = 0.5 model can best reproduce the radial temperature profile in the most sophisticated thermal scheme. Nevertheless, the inner disk regions are colder in this β-model compared to (TheS2)-model 2v. Strong β-cooling in the innerfast-rotating disk regions completely overwhelms disk viscous heating in the lowest-β model. Interestingly, all β-models show a localincrease in the gas temperature in the outer disk regions and in the inner envelope, a feature that is also present in the ThES2 scheme, but is absent in the ThES1 scheme (see Fig. 17). This may be related to strong compressional heating due to PdV work of infalling envelope material and reduced β-cooling in the outerregions with slow rotation. The amplitude of the temperature jump in the β-models is a factor of several higher than in (ThES2)-model 2v.

We conclude that β-models with irradiation can more closely reproduce the thermal state obtained in the ThES2 scheme than the β-models without irradiation, but the agreement is still not acceptable. We note that the required computational resources for the ThES2 scheme are significantly higher than those for the ThES1 and β-schemes. Therefore, the former scheme is advisable for simulations where gas and dust temperature decoupling is expected to be of particular importance. In other situations, the ThES1 scheme may be a method of choice due to its relatively easy coding and the moderate computational resources required.

Fig. 16 Comparison of gas temperature distributions in models with different thermal evolution schemes. First and second rows: gas temperatures in (ThES2)-model 2v and (ThES1)-model 2v, respectively. Third, fourth, and last rows: gas temperatures in (beta=3)-model 2v, (beta=10)-model 2v, and (beta=30)-model 2v, respectively. In the β-models the gas and dust temperatures are the same. The time is counted from the instance of disk formation. The scale bar is in log Kelvin.

Fig. 17 Comparison of the azimuthally averaged gas temperatures in models with different thermal evolution schemes. The red and black solid lines present the gas temperature profiles for ThES2 and ThES1 in model 2v, respectively. The blue solid, dashed, and dotted lines show the profiles for model 2v with β = 3, β = 10, and β = 30, respectively.

Fig. 18 Values of the β-parameter calculated in (ThES1)-model 2 at 2.8 kyr after the disk formation instance. Each red dot corresponds to a β-value in an individual grid cell. Significant azimuthal and radial variations are evident.

Fig. 19 Comparison of gas surface density distributions in model v2 with ThES2 and in β-models that takestellar and background irradiation into account. First row: (ThES2)-model 2v, second row: (beta=0.5Ir)-model 2v, third row: (beta=3Ir)-model 2v, and bottom row: (beta=10Ir)-model 2v. The red circles outline the radial distance beyond which the azimuthally averaged angular velocity deviates from the Keplerian rotation bymore than 10%. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.

Fig. 20 Comparison of the number of fragments at a given time instance in the disks of (ThES2)-model 2v (top panel), (beta=0.5Ir)-model 2v (middle panel), and (beta=1.0Ir)-model 2v (bottom panel). The time is counted from the instance of disk formation.

Fig. 21 Comparison of the azimuthally averaged gas temperatures in (TheS2)-model 2v (red lines) and in β-models with irradiation. The blue solid, dashed, and dotted lines show the profiles for model 2v with β = 0.5, β = 3, and β = 10, respectively.

6 Summary and discussion

Here we explored numerically the global long-term evolution of protostellar disks with different cooling and heating schemes. For this purpose, we used three approaches to describe the thermal balance in the disk: a simplistic β-cooling scheme with and without stellar and background irradiation considered, a more realistic scheme with similar gas and dust temperatures (ThES1), and a sophisticated scheme with separate gas and dust temperatures (ThES2). The last scheme can also be applied to low-metallicity protostellar disks.

The adopted thermal schemes were tested on global disk models computed in the thin-disk limit without and with turbulent viscosity using the Shakura & Sunyaev α-parameterization. The main results can be summarized as follows:

• In the ThES2 scheme, the gas and dust temperatures are similar in the inner high-density regions of the disk, but may significantly deviate from each other in the vicinity of the disk outer edge. In this low-density region a pronounced decoupling between the dust and gas temperatures develops thanks to additional compressional heating of the gas caused by the infalling envelope material and because of slow collisional energy exchange between gas and dust in the low-density environment. The gas temperature may exceed that of dust by tens or even hundreds of Kelvin.

• In the outer circumdisk environment occupied by a rarefied envelope the gas temperature also decouples from that of dust, but in the opposite direction: the gas temperature drops below that of dust. This effect was also reported in Pavlyuchenkov et al. (2015) and Bate (2018).

• The global disk evolution is weakly sensitive to decoupling of gas and dust temperatures. Gravitational instability in the case of separate gas and dust temperatures (ThES2) is only slightly stronger than in the case of similar gas and dust temperatures (ThES1), asindicated by higher Fourier amplitudes, and gravitational fragmentation is slightly more frequent. Overall, separate gas and dust temperatures do not cause qualitative changes to the disk evolution, which appears to be more sensitive to the presence or absence of turbulent viscosity in the disk.

• Decoupling of gas and dust temperatures may nevertheless be of significance for the chemical evolution and dust growth. An increase in the gas temperature to more than 100 Kelvin could launch gas phase reactions in the disk outer regions that otherwise may be dormant. Moreover, decoupling of gas and dust temperatures may also facilitate the growth of icy mantles on cold dust particles if volatile species become oversaturated in the warm gas environment in the vicinity of the disk outer edge. The decoupled gas and dust temperatures may also affect the disk mass estimates.

• Simplistic constant-β models without irradiation fail to reproduce the disk evolution in more sophisticated thermal models with or without separate gas and dust temperatures. We attribute this to the intrinsic variability of the β-parameter both in time and space. β-cooling models with stellar and background irradiation taken into account can better match the dynamical and thermal evolution obtained in the ThES2 thermal scheme, particularly for β ≈ 0.5−1.0, but the agreement is still incomplete.

In the future, it will be interesting to consider cases when the gas and dust temperatures decouple through the bulk of the disk, and not only in the disk outermost regions and in the envelope. This may occur either at metallicities lower than solar or when the dust-to-gas ratio drops below the canonical 1:100 value in otherwise solar metallicity disks.

Acknowledgements

We are thankful to the anonymous referee for constructive comments that helped to improve the manuscript. E.I.V. and M.G. acknowledge support from the Austrian Science Fund (FWF) under research grant P31635-N27. K.O and R.M acknowledge support work by MEXT/JSPS KAKENHI Grant Number17H01102, 17H02869, 17H06360. The simulations were performed on the Vienna Scientific Cluster.

Appendix A Line-averaged escape probability of HD

The HD line cooling rate is reduced by the effect of self-absorption when the column density of HD in the disk is large. The photon escape probability for self-absorption in individual transitions is given by $${\beta }_{\text{HD},\text{ul}}=\frac{1-{\text{e}}^{-{\tau }_{\text{ul}}}}{{\tau }_{\text{ul}}},$$(A.1)

where τ_ul is the optical depth for line emission given by the same formula of Eq. (30). The photon escape probability is a function of the column density of HD and the gas temperature because the optical depth for line emission depends on these values, as you can see from Eqs. (30) and (32). We treat four levels in the range of rotational levels 0 < J < 3. The level energies and the spontaneous radiative decay rates are from Dalgarno & Wright (1972) and Abgrall et al. (1982), respectively.We only consider the helium impact as the collisional deexcitation process and its rate coefficients are taken from Galli & Palla (1998). The cooling rate of HD line emission is calculated by counting level populations from the statistical balance among four levels. We calculate the cooling rate in the range of the column density of HD 10¹⁵   cm⁻² < N_HD < 10²⁵ cm⁻² and the gas temperature30 K < T_g < 3000 K. We fit the line-averaged escape probability obtained from the results in a functional form as $${\overline{\beta }}_{\text{esc},\text{HD}}=\frac{1}{{\left(1+N_{\text{HD}}/N_{\text{c}}\right)}^{\alpha }},$$(A.2)

where $$\alpha =a_0+a_1{T}_{\text{g}}+a_2{T}_{\text{g}}^2+a_3{T}_{\text{g}}^3+a_4{T}_{\text{g}}^4+a_5{T}_{\text{g}}^5,$$(A.3)

and $$\begin{array}{lll}\text{log} N_{\text{c}}\hfill & =\hfill & b_0+b_1\text{log} T_{\text{g}}+b_2{\left(\text{log} T_{\text{g}}\right)}^2+b_3{\left(\text{log} T_{\text{g}}\right)}^3\hfill \\ \hfill & \hfill & +b_4{\left(\text{log} T_{\text{g}}\right)}^4+b_5{\left(\text{log} T_{\text{g}}\right)}^5.\hfill \end{array}$$(A.4)

The fitting coefficients in Eqs. (A.3) and (tabfloatlinking A.4) are given in Tables A.1 and A.2.

Appendix B Chemical reactions

In our new cooling–heating scheme described in Sect. 2.3, we follow the non-equilibrium chemical evolution. Our chemical network is composed of 21 hydrogen reactions and 6 deuterium reactions. We describe the treated chemical reactions and their rate coefficients in Table B.1. In reactions from 1 to 20, only the rate coefficients of forward reaction are shown. The calculation method of the rate coefficients of reverse reactions is from Matsukoba et al. (2019).

References

All Tables

All Figures

	Fig. 1 Gas surface density distributions in the five models considered. Each row presents a specific model, as indicated, and each column corresponds to a specific time starting from disk formation. The scale bar is in log g cm−2.
In the text

	Fig. 2 Similar to Fig. 1, but for the dust temperature. The scale bar is in log K.
In the text

	Fig. 3 Similar to Fig. 1, but for the gas temperature. The scale bar is in log K.
In the text

	Fig. 4 Gas velocity field superimposed on the gas temperature distribution in (ThES2)-model 2v. The black contour outlines a gas surface density of 0.1 g cm−2.
In the text

Fig. 5 Decoupling of gas and dust temperatures in the disk and envelope in (ThES2)-model 2v at t = 0.16 Myr. Shown are the threshold temperature (Tcrit, black thick line) above and below which the gas and dust temperatures are thermally decoupled and coupled, respectively. The red thin and thick lines present the minimum and maximum azimuthal variations of the gas temperature, respectively,while the blue dashed thin and thick lines show the corresponding quantities for the dust temperature. The brown circleindicates the position of the disk outer edge.

In the text

	Fig. 6 Comparison of the gas surface density, and gas and dust temperatures in model 2. Top and bottom rows: (ThES2)-model 2v and (ThES2)-model 2, middle row: results for a test model with the rate of collisional heat exchange between dust and gas Γcoll increased artificially by a factor of 50.
In the text

Fig. 7 Azimuthally averaged radial distributions of temperatures for gas (red solid line), dust (blue dash-dotted line), and radiation (green dotted line) in five models considered. The shaded areas indicate the range of azimuthal variations of gas (pink) and dust (blue) temperatures at each radius. The black line is the difference between the maximum temperatures of gas and dust (right-hand axis). The gray dashed line shows the radius of the disk outer edge. The arrangement of panels is the same as in Fig. 1.

In the text

	Fig. 8 Comparison of gas surface density distributions in model 2v with ThES2 (separate dust and gas temperatures) and ThES1 (equal gas and dust temperatures). The scale bar is in log g cm−2.
In the text

	Fig. 9 Comparison of the global Fourier amplitudes in model 2v with ThES1 (red lines) and ThES2 (black lines). The global amplitudes for four modes (m = 1, 2, 3, and 4) are shown in the four panels. The time is counted from the beginning of the core collapse.
In the text

	Fig. 10 Comparison of the number of fragments in the disk in model 2v with ThES2 (top panel) and ThES1 (bottom panel). The time is counted from the instance of disk formation.
In the text

Fig. 11 Comparison of gas and dust temperature distributions in model v2 with ThES1 and ThES2. First and second rows: gas and dust temperature distributions in ThES2. Third row: temperature distribution (the same for gas and dust) in ThES1. Bottom row: comparison of the azimuthally averaged temperatures in the two considered thermal evolution schemes. The solid lines show the temperatures of gas and dust, while the dashed lines provide the temperatures of stellar irradiation. The scale bar is in log Kelvin.

In the text

	Fig. 12 Comparison of the azimuthally averaged gas and dust temperature distributions in model v2 with ThES1 and ThES2. Thesolid lines show the temperatures of gas and dust, while the dashed lines provide the temperatures of stellar irradiation. The gas and dust temperatures in the (ThES2)-model v2 coincide everywhere except for the outer regions.
In the text

	Fig. 13 Comparison of gas surface density distributions in model v2 with ThES2 and β-cooling. First row:(ThES2)-model 2v, second row: (beta=3)-model 2v, third row: (beta=10)-model 2v, and bottom row: (beta=30)-model 2v. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.
In the text

Fig. 14 Gas velocity field superimposed on the gas surface density distribitions in (ThES2)-model 2v (top row), (beta=3)-model 2v (middle row), and (beta=10)-model 2v (bottom row) at two evolutionary times, as indicated in each panel. The red circles outline the radial distance beyond which the azimuthally averaged angular velocity deviates from the Keplerian rotation by more than 10%. The green contours outline the regions where the gas surface density drops to 0.1 g cm−1. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.

In the text

	Fig. 15 Comparison of the number of fragments at a given time instance in the disks of (ThES2)-model 2v (top panel), (beta=3)-model 2v (middle panel), and (beta=10)-model 2v (bottom panel). The time is counted from the instance of disk formation.
In the text

Fig. 16 Comparison of gas temperature distributions in models with different thermal evolution schemes. First and second rows: gas temperatures in (ThES2)-model 2v and (ThES1)-model 2v, respectively. Third, fourth, and last rows: gas temperatures in (beta=3)-model 2v, (beta=10)-model 2v, and (beta=30)-model 2v, respectively. In the β-models the gas and dust temperatures are the same. The time is counted from the instance of disk formation. The scale bar is in log Kelvin.

In the text

	Fig. 17 Comparison of the azimuthally averaged gas temperatures in models with different thermal evolution schemes. The red and black solid lines present the gas temperature profiles for ThES2 and ThES1 in model 2v, respectively. The blue solid, dashed, and dotted lines show the profiles for model 2v with β = 3, β = 10, and β = 30, respectively.
In the text

	Fig. 18 Values of the β-parameter calculated in (ThES1)-model 2 at 2.8 kyr after the disk formation instance. Each red dot corresponds to a β-value in an individual grid cell. Significant azimuthal and radial variations are evident.
In the text

Fig. 19 Comparison of gas surface density distributions in model v2 with ThES2 and in β-models that takestellar and background irradiation into account. First row: (ThES2)-model 2v, second row: (beta=0.5Ir)-model 2v, third row: (beta=3Ir)-model 2v, and bottom row: (beta=10Ir)-model 2v. The red circles outline the radial distance beyond which the azimuthally averaged angular velocity deviates from the Keplerian rotation bymore than 10%. The time is counted from the instance of disk formation. The scale bar is in log g cm−2.

In the text

	Fig. 20 Comparison of the number of fragments at a given time instance in the disks of (ThES2)-model 2v (top panel), (beta=0.5Ir)-model 2v (middle panel), and (beta=1.0Ir)-model 2v (bottom panel). The time is counted from the instance of disk formation.
In the text

	Fig. 21 Comparison of the azimuthally averaged gas temperatures in (TheS2)-model 2v (red lines) and in β-models with irradiation. The blue solid, dashed, and dotted lines show the profiles for model 2v with β = 0.5, β = 3, and β = 10, respectively.
In the text