Efficiency of thermal relaxation by radiative processes in protoplanetary discs: constraints on hydrodynamic turbulence

TL;DR

This study investigates how radiative thermal relaxation influences the development of hydrodynamic turbulence in protoplanetary discs, providing constraints on where different instabilities can occur based on disc properties.

Contribution

We estimate radiative relaxation timescales using updated opacities and apply them to localize regions of hydrodynamic instabilities in protoplanetary discs.

Findings

Most of a typical T Tauri disc is unstable to hydrodynamic instabilities.

Convective overstability occurs between 2 au and 40 au, and beyond 140 au.

Vertical shear instability develops between 15 au and 180 au.

Abstract

Hydrodynamic, non-magnetic instabilities can provide turbulent stress in the regions of protoplanetary discs, where the MRI can not develop. The induced motions influence the grain growth, from which formation of planetesimals begins. Thermal relaxation of the gas constrains origins of the identified hydrodynamic sources of turbulence in discs. We estimate the radiative relaxation timescale of temperature perturbations and study the dependence of this timescale on the perturbation wavelength, the location within the disc, the disc mass, and the dust-to-gas mass ratio. We then apply thermal relaxation criteria to localise modes of the convective overstability, the vertical shear instability, and the zombie vortex instability. Our calculations employed the latest tabulated dust and gas mean opacities and we account for the collisional coupling to the emitting species. The relaxation criterion defines the bulk of a typical T Tauri disc as unstable to the development of linear hydrodynamic instabilities. The midplane is unstable to the convective overstability from at most $2\mbox{ au}$ and up to $40\mbox{ au}$, as well as beyond $140\mbox{ au}$. The vertical shear instability can develop between $15\mbox{ au}$ and $180\mbox{ au}$. The successive generation of (zombie) vortices from a seeded noise can work within the inner $0{.}8\mbox{ au}$. Dynamic disc modelling with the evolution of dust and gas opacities is required to clearly localise the hydrodynamic turbulence, and especially its non-linear phase.

Equations

$$t_{\mathrm{relax},k}=\frac{\left|\delta T_k\right|}{\delta {\dot{T}}_k}.$$

$$t_{\mathrm{relax}}=\max \left(t_{\mathrm{thick}},t_{\mathrm{thin}}\right).$$

$$\frac{{\rho }_{\mathrm{gas}}}{{\rho }_{\mathrm{dust}}}={\eta }^{-1}\sim 10^2.$$

$$t_{\mathrm{emit}}=\frac{C_{\mathrm{v}}}{16\max \left({\kappa }_{\mathrm{P}}^{\mathrm{d}},{\kappa }_{\mathrm{P}}^{\mathrm{g}}\right)\sigma T^3}$$

$$t_{\mathrm{coll}}=\frac{1}{n\sigma v}.$$

$$t_{\mathrm{relax}}^{d/g}=\max \left(t_{\mathrm{emit}}^{d/g},t_{\mathrm{coll}}^{d/g}\right).$$

$$t_{\mathrm{thin}}=\min \left(t_{\mathrm{relax}}^{\mathrm{d}},t_{\mathrm{relax}}^{\mathrm{g}}\right).$$

$$l_{\mathrm{diff}}=\frac{1}{{\kappa }_{\mathrm{R}}\rho },$$

$$t_{\mathrm{thick}}=\frac{1}{\tilde{D}{k}^2}\propto {\tau }_{\mathrm{R}}\lambda \rho T^{-3},$$

$${\tau }_{\mathrm{R}}={\int }_{-\lambda /2}^{+\lambda /2}{\kappa }_{\mathrm{R}}\rho dl.$$

$$\Sigma ={\Sigma }_{\mathrm{c}}{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{-p}\exp \left[-{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{2-p}\right],$$

$${\Sigma }_{\mathrm{c}}=\left(2-p\right)\frac{M_{\mathrm{disc}}}{2\pi R_{\mathrm{c}}^2}.$$

$$\rho =\frac{\Sigma }{\sqrt{2\pi }H}\exp \left(\frac{-z^2}{2H^2}\right),$$

$$T=T_0{\left(\frac{R_{\mathrm{in}}}{R}\right)}^{1/2},$$

$${\kappa }^{\mathrm{d}}=\frac{{\eta }^{\prime }}{\eta }{\kappa }^{\mathrm{d}},$$

$${\int }_{-\lambda /2}^{+\lambda /2}{\kappa }_{\mathrm{R}}\rho dl={\tau }_{\mathrm{R}}=\frac{{\left(2\pi \right)}^2}{3}\frac{⟨{\kappa }_{\mathrm{R}}⟩}{{\kappa }_{\mathrm{P}}},$$

$${\tau }_{\mathrm{R},\lambda /2}={\int }_{-\lambda /4}^{+\lambda /4}{\kappa }_{\mathrm{R}}\rho dl.$$

$${\tau }_{\mathrm{P}}^{\ast }=\int \left({\kappa }_{\mathrm{P}}^{\mathrm{g}}\left(T_{\ast }\right)+{\kappa }_{\mathrm{P}}^{\mathrm{d}}\left(T_{\ast }\right)\right)\rho dR=\int {\kappa }_{\mathrm{P}}^{\mathrm{g}}\left(T_{\ast }\right)\rho dR.$$

$${\tilde{\tau }}_{\mathrm{P}}^{\ast }=\int {\kappa }_{\mathrm{P}}^{\mathrm{g}}\left(T_{\ast },T_{\mathrm{g}}\right)\rho dR.$$

$$t_{\mathrm{thick}}\propto {\tau }_{\mathrm{R}}\lambda \propto {\lambda }^2,$$

$$t_{\mathrm{thick}}\propto \lambda ,$$

$$\Omega t_{\mathrm{c}}=\frac{h|q|}{\gamma -1},$$

$${\partial }_t\left(E_{\mathrm{int}}+E_{\mathrm{R}}\right)=-\nabla \text{@vec}{F}_{\mathrm{R}}.$$

$$E_{\mathrm{int}}=\rho C_{\mathrm{v}}T,$$

$$E_{\mathrm{R}}=a{T}^4,$$

$${\partial }_t\left(E_{\mathrm{int}}+E_{\mathrm{R}}\right)=4a{T}^3\left(1+\frac{\rho c_{\mathrm{V}}}{4a{T}^3}\right){\partial }_{t}T=4a{T}^3\left(\frac{3}{4}+\frac{1}{4}{\eta }^{-1}\right){\partial }_{t}T.$$

$$\text{@vec}{F}_{\mathrm{R}}=-D\nabla E_{\mathrm{R}}.$$

$$D=\frac{\varkappa c}{{\kappa }_{\mathrm{R}}\rho },$$

$$-\nabla \text{@vec}{F}_{\mathrm{R}}=4a{T}^{3}D\times \left(△T+\frac{3{\left(\nabla T\right)}^2}{T}+\text{@vec}A\cdot \nabla T\right),$$

$$\text{@vec}A=\frac{\nabla {\kappa }_{\mathrm{R}}}{{\kappa }_{\mathrm{R}}}+\frac{\nabla \rho }{\rho }=\nabla \ln {\kappa }_{\mathrm{R}}+\nabla \ln \rho =\nabla \ln \rho \left[1+\rho {\partial }_{\rho }\ln {\kappa }_{\mathrm{R}}\right].$$

Full text

*Aims. * We aim to estimate the radiative relaxation timescale of temperature perturbations in protoplanetary discs. We study the dependence of the thermal relaxation on the perturbation wavelength, the location within the disc, the disc mass, and the dust-to-gas mass ratio. We then apply thermal relaxation criteria to localise modes of the convective overstability, the vertical shear instability, and the zombie vortex instability.

*Methods. * For a given temperature perturbation, we estimated two timescales: the radiative diffusion timescale $t_{\mathrm{thick}}$ and the optically thin emission timescale $t_{\mathrm{thin}}$. The longest of these timescales governs the relaxation: $t_{\mathrm{relax}}=\max\left(t_{\mathrm{thick}},t_{\mathrm{thin}}\right)$. We additionally accounted for the collisional coupling to the emitting species. Our calculations employed the latest tabulated dust and gas mean opacities.

*Results. * The relaxation criterion defines the bulk of a typical T Tauri disc as unstable to the development of linear hydrodynamic instabilities. The midplane is unstable to the convective overstability from at most $2\mbox{ au}$ and up to $40\mbox{ au}$, as well as beyond $140\mbox{ au}$. The vertical shear instability can develop between $15\mbox{ au}$ and $180\mbox{ au}$. The successive generation of (zombie) vortices from a seeded noise can work within the inner $0{.}8\mbox{ au}$.

*Conclusions. * A map of relaxation timescale constrains the origins of the identified hydrodynamic turbulence-driving mechanisms in protoplanetary discs. Dynamic disc modelling with the evolution of dust and gas opacities is required to clearly localise the hydrodynamic turbulence, and especially its non-linear phase.

1 Introduction

Despite the ample observational manifestation, accretion discs around young stars remain a puzzle. In such environments, the major mode of planet formation takes place (Johansen et al. 2014; Raymond et al. 2014; Chabrier et al. 2014). The exact mechanism of turbulence generation in circumstellar discs is still debated (Turner et al. 2014) and there may be multiple mechanisms operating simultaneously in different locations of a disc yet interacting with each other (Armitage 2015). From the theoretical perspective, it is clear that turbulence is central both for angular momentum transport and planetesimal aggregation (Turner et al. 2014).

One of the most promising models of turbulence in protoplanetary discs is the magneto-rotational instability (MRI, Balbus & Hawley 1991). The MRI operates in weakly magnetised yet well-coupled shearing flows. The major localising constraint on the MRI is a low ionisation degree of the gas (Gammie 1996; Dzyurkevich et al. 2013). Accretion in the poorly ionised shearing magnetic flows is dominated by the non-ideal magnetohydrodynamic (MHD) effects with the accretion flow itself being not necessarily turbulent. Lesur et al. (2014); Bai (2015); Bai & Stone (2017) have shown in three-dimensional MHD modelling that the Hall effect can largely control the dynamics of the laminar midplane flow yielding efficient wind accretion. Still, turbulent motions in MRI-inactive regions (’dead zones’) can develop via hydrodynamic or gravitational instabilities. These require specific thermodynamic conditions, which can be used to localise their origins in discs. The requirement for an appropriate thermal relaxation time is a necessary but not sufficient condition for the hydrodynamic instabilities to develop: the growth rates are affected by the stresses associated with the non-ideal MHD effects. There are currently three identified purely hydrodynamic111we do not discuss the gravitational instability here because of the controversy over the cooling time criterion(Meru & Bate 2011, 2012; Paardekooper 2012; Hopkins & Christiansen 2013; Baehr & Klahr 2015; Takahashi et al. 2016) sources of turbulent motions in protoplanetary discs: the vertical shear instability222also known as the Goldreich-Schubert-Fricke instability in the context of differentially rotating stars (Goldreich & Schubert 1967; Fricke 1968)(Urpin & Brandenburg 1998; Urpin 2003; Nelson et al. 2013; Barker & Latter 2015; Lin & Youdin 2015), the convective overstability in its linear (Klahr & Bodenheimer 2003; Klahr & Hubbard 2014; Lyra 2014; Latter 2016) and non-linear state (subcritical baroclinic instability, Klahr & Bodenheimer 2003; Petersen et al. 2007; Lesur & Papaloizou 2010), and the zombie vortex instability (Marcus et al. 2013, 2015). The first two are linear instabilities: they can drive turbulence out of an initial equilibrium. The zombie vortex instability is subcritical: it demands initial, finite-amplitude noise. The onset, linear growth, and non-linear phase of these hydrodynamic perturbations depend on thermal relaxation because of the buoyant nature of the restoring forces.

A careful localisation of the transport mechanisms is required for modelling disc-planet interaction (Lobo Gomes et al. 2015), core migration (Matsumura et al. 2007, 2009; Mordasini et al. 2015), and disc transport studies. It is further required to make use of the growing body of high-resolution, sensitive observations (Pohl et al. 2015), which are especially expanded by the Atacama Large Millimeter Array and extreme adaptive optics systems like Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE) at the VLT.

Furthermore, during turbulent transport, both the induced small-scale motions (Brauer et al. 2008; Zsom et al. 2010) and the gas drag in local density enhancements (Whipple 1972; Barge & Sommeria 1995; Klahr & Lin 2001; Klahr & Bodenheimer 2006) can affect the dust grain growth (Birnstiel et al. 2012) and influence the orbital migration of protoplanets (Nelson & Papaloizou 2004; Ida et al. 2008). Dust evolution in a turbulent disc with an ongoing planetesimal formation changes the dust density and its contribution to the opacity of the disc (Johansen et al. 2006; Birnstiel et al. 2012). The migration of protoplanets depends on the ability of the surrounding disc material to liberate the excess thermal energy, that is, the thermal relaxation (Lyra et al. 2010; Dittkrist et al. 2014). For example, the so-called type I migration of low-mass embedded cores in locally isothermal discs leads to their fast (on a timescale shorter than the disc lifetime, Korycansky & Pollack 1993) depletion. A finite relaxation rate, however, can lead to a reversal of the migration direction (Paardekooper & Mellema 2006; Kley & Crida 2008; Kley et al. 2009) by enhancing the corotation torque (Baruteau & Masset 2008).

The computational challenge of incorporating radiative transfer into dynamical modelling as well as the uncertainties in opacities and their secular change precludes the accurate treatment of energy evolution in accretion discs. A sensible parametrisation for numerical studies is the thermal relaxation timescale in the energy evolution equation. 3D hydrodynamic simulations with full radiative transport necessitate a realistic a priori estimate of the relaxation timescale as a function of disc parameters. Both the set-up of the numerical experiments and their physical interpretation require effective thermal relaxation timescale.

In this paper, we present a simple method to estimate the relaxation timescale of linear temperature perturbations (Sect. 2). The formulae were derived from non-compressive linear analysis of the energy evolution equation. Given the mean opacities and a predefined disc structure (Sect. 3), the technique allows mapping of the radiative relaxation timescale over the disc interiors (Sect. 4). Such maps directly constrain the origins of the hydrodynamic turbulence. We discuss the restrictions of the method and measures for further improvement in Sect. 5. A brief summary in Sect. 6 closes the paper.

2 Radiative thermal relaxation

Thermal relaxation timescale is a characteristic evolution timescale of linear temperature perturbations. In such dilute media as protoplanetary discs (PPDs), the relaxation is chiefly supplied by means of radiative transport. It is convenient to express the thermal relaxation in terms of a dimensionless parameter $\beta=t/\Omega^{-1}$, where $\Omega^{-1}$ is the sound crossing time over one pressure scale height $H$.

To isolate the characteristic timescale of the relaxation, we have considered the evolution of only the temperature (actually, thermal and radiation energy densities, see Appendix A.3). We neglected the associated changes in pressure and velocity field (cf. Lyra 2014). Linearising the energy evolution equation, we find the characteristic decay timescale of a Fourier mode $\delta T_{k}=(T-T_{\mathrm{0}})_{k}$:

[TABLE]

Here $T$ is the perturbed temperature, $T_{\mathrm{0}}$ the equilibrium temperature ($\dot{T}_{\mathrm{0}}=0$) and $\delta\dot{T}$ designates the time derivative of the perturbation.

To estimate the radiative relaxation timescale, we calculated the characteristic timescales of the underlying processes. The relevant processes differ depending on the optical thickness of the mode (the optical thickness of matter encompassed within a spatial Fourier wavelength). An appropriate frequency-averaged measure of this optical thickness is the Rosseland optical depth (see Sect. 2.2). This distinguishes optically thick and optically thin relaxation regimes with characteristic timescales $t_{\mathrm{thick}}$ and $t_{\mathrm{thin}}$, respectively. Because the underlying processes operate simultaneously, one can approximate the relaxation timescale by the longest timescale amongst the two (cf. Hubeny 1990):

[TABLE]

Below we introduce each of these timescales.

2.1 Optically thin relaxation regime

The thermal energy of the gas in a PPD dominates the thermal energy of the dust due to the high gas-to-dust mass fraction:

[TABLE]

However, the dust emissivity, as measured by the Planck mean opacity, is comparable or exceeds that of the gas at low temperatures. The radiative efficiency of the gas with the specific thermal capacity $C_{\mathrm{v}}$ is expressible in terms of the thermal emission timescale (see Appendix B.1)

[TABLE]

with $\sigma$ being the Stefan-Boltzmann constant, $\kappa^{\mathrm{d/g}}_{\mathrm{P}}$ the Planck mean opacity of the dust (per unit gas mass) or the gas, respectively. Timescale (4) is evaluated locally, and, therefore, does not depend on perturbation wavelength. This is a sensible approximation for optically thin modes. The optically thin emission timescale as given by Eq. (4) is the shortest achievable relaxation timescale for the given material in local thermodynamic equilibrium (LTE).

At dust temperatures in the range $650<T<1\,300\mbox{ K}$, when the abundant volatile materials are already evaporated, the values of the Planck mean opacities of the dust and the gas are comparable (Malygin et al. 2014; Semenov et al. 2003). This means that the gas dominates the thermal emission. Here, we have approximated the medium as being made of three sorts of species: first – the most abundant poorly emitting gas (thermal carriers: H2, He); second – the remaining gas (including emitters, like H2O, TiO, SiO, CO, etc.), and third – the dust. Collision with electrons are taken into account in the equilibrium radiative processes via the Planck opacities. The bulk of thermal energy of the gas is stored in the kinetic energy of the most abundant H2 molecules, which have no permanent dipole moment333The less abundant HD isotopologue, though, has a permanent dipole moment that permits its detection (Bergin et al. 2013) and can transfer their energy to the species with higher emissivity (both dust and gas, see Appendix B.2) via collisions. The speed of this coupling was measured by respective collision timescale

[TABLE]

In a linear approximation (small temperature perturbations) thermal capacity of neither the dust nor the emitting gas is exhausted. Thus, the gas and the dust emitters each set a characteristic relaxation timescale

[TABLE]

If $t^{\mathrm{d}}_{\mathrm{coll}}>t^{\mathrm{d}}_{\mathrm{emit}}$ (or $t^{\mathrm{g}}_{\mathrm{coll}}>t^{\mathrm{g}}_{\mathrm{emit}}$), we say that the dust (or gas) is collisionally decoupled. Relaxation via dust and gas channels operates simultaneously, the fastest channel determining the optically thin relaxation of the whole material:

[TABLE]

In the so-called molecular layers above the midplane, one has to account for photoelectric heating, mutual irradiation of the dust and the gas, photochemistry, etc. (Kamp & van Zadelhoff 2001; Akimkin et al. 2013; Henning & Semenov 2013). Those are not taken into account in this study, where we have focussed on the bulk amount of disc gas in the midplane.

2.2 Optically thick relaxation regime

A perturbation with wavelength exceeding the diffusion length scale

[TABLE]

can not relax as fast as the local thermal emission timescale because the thermal photons get absorbed and re-emitted multiple times before escaping the perturbed medium. Relaxation timescale of mode $\lambda$ can be approximated by the corresponding diffusion timescale (see Appendix A)

[TABLE]

with the effective opacity in the frequency-averaged diffusion coefficient $\tilde{D}\propto(\kappa\rho)^{-1}$ being the Rosseland opacity $\kappa_{\mathrm{R}}$ and

[TABLE]

being the Rosseland thickness over $\lambda$. One expects the optically thick timescale (9) to approximate the relaxation timescale (2), if the medium encompassed within the perturbation wavelength is optically thick, $\tau_{\mathrm{R}}>1$ (see Sect. 4.1 for a refinement).

3 Set-up

3.1 Density model

We have adopted an axisymmetrical disc density structure from Dzyurkevich et al. (2013) to allow comparison with the MHD studies on the magnetically dead zone. The central source is a solar-like star with $M_{\mathrm{*}}=1M_{\odot}$. The disc has a surface density profile

[TABLE]

with $\Sigma_{\mathrm{c}}$ being the surface density at the cut-off radius $R_{\mathrm{c}}$:

[TABLE]

The surface density is normalised to $\Sigma=1700\mbox{ g cm}^{-2}$ at $1\mbox{ au}$. Index $p$ of $0{.}9$ and the cut-off radius of $40\mbox{ au}$ is chosen following Andrews et al. (2009) to match the minimum-mass solar nebula at $1$ and $100\mbox{ au}$. $M_{\mathrm{disc}}$ is the total mass of the disc. The volumetric density is readily unfolded for a vertically isothermal case

[TABLE]

where $H$ is the local pressure scale height.

3.2 Temperature model

We adopt a vertically isothermal radially stratified temperature profile following Dzyurkevich et al. (2013)

[TABLE]

with $T_{\mathrm{0}}=280\mbox{ K}$ at $R_{\mathrm{in}}=1\mbox{ au}$. Since we focus on the bulk disc matter in the midplane region, we ignore upper disc layers and hence vertical $T$-gradient. A lower limit of $10\mbox{ K}$ stops the power law decrease at $\sim 250\mbox{ au}$. This is a typical lowest temperature in outer midplane of the discs around Sun-like stars predicted by detailed radiative transfer models, such as RADMC-3D (Dullemond 2012).

The profiles of surface density, volumetric density, and the midplane temperature are given in Fig. 1. The highest value of $M_{\mathrm{disc}}=0{.}1M_{\odot}$ considered for this model assures the configuration is stable against self-fragmentation, according to the Toomre $Q$-criterion.

3.3 Opacities

We used tabulated dust and gas mean opacities (Fig. 2). The dust opacity is that of iron-poor homogeneous silicate (IPS) spheres from Semenov et al. (2003), the gas opacity is that of the solar-mixture gas from Malygin et al. (2014). The total opacity is the direct sum of the gas and the dust component. Figures 2a,c display the combined data. Gas chemistry and opacity in outer disc regions are poorly constrained observationally (see Dutrey et al. 2014, and references therein). In calculations with significant dust depletion, we considered two cases of pure gas opacity: (i) data from Malygin et al. (2014) with the $650\mbox{ K}$ cut-off at low temperatures; (ii) data from Malygin et al. (2014) complemented with scaled-down (in density) atmospheric opacities from Freedman et al. (2008) with the lower boundary of $75\mbox{ K}$. For yet lower temperatures, the value at the lowest available temperature is always used. Data from Freedman et al. (2008) accounts for the settling of condensates in a gravitational field (e.g. in cool dwarves) but neither Freedman et al. (2008) nor Malygin et al. (2014) account for the freeze out on grain surfaces. However, the data from Semenov et al. (2003) takes into account freeze out of volatiles such as organics, H2O, and CO, the main absorbers at cold temperatures in the gas phase.

Planck opacities from Freedman et al. (2008) have relatively high values in the interval $180-320\mbox{ K}$ and a steep decline below $\approx 120\mbox{ K}$ (Fig. 2b). Rosseland opacities from Freedman et al. (2008) exceed those from Malygin et al. (2014) in the interval $480-700\mbox{ K}$ and then drop by $\sim 3$ orders of magnitude at $75\mbox{ K}$ (Fig. 2d). We will refer to the data complemented with Freedman et al. (2008) as the small gas opacity case and to the data extrapolated from Malygin et al. (2014) as the large gas opacity case. The two models can be treated as rough borders of the poorly-constrained parameter space. One of the motivations beyond the large gas opacity case can be an increase in the gas-phase opacity (e.g. a release of molecules not taken into account in Freedman et al. (2008)) after local dust evaporation.

When varying the dust-to-gas mass ratio (Sect. 4.4), we modified the dust opacity by

[TABLE]

where $\eta^{\prime}$ is the new dust-to-gas mass ratio and $\eta=0{.}014$ the nominal dust-to-gas mass ratio, for which the dust opacity tables were calculated (see Table 1 in Semenov et al. 2003, case of IPS). The gas opacity stays unchanged.

4 Results

In this section, the values for the relaxation rate are dimensionless numbers $\Omega t$, with $\Omega$ being the local Keplerian frequency. The radiative diffusion timescale is calculated for diffusion in the vertical direction unless explicitly stated otherwise. The spatial perturbation wavelength is conveniently measured in local pressure scale heights: $\tilde{\lambda}=\lambda/H$, or $\tilde{k}=2\pi\tilde{\lambda}^{-1}=kH$.

4.1 Relaxation regimes: verification

One expects the transition from the optically thin to the optically thick relaxation to happen when the perturbation wavelength $\lambda$ becomes of the order the diffusion length scale $l_{\mathrm{diff}}$. Equating (4) and (9) while letting $\lambda=l_{\mathrm{diff}}$ gives

[TABLE]

where $\langle\kappa_{\mathrm{R}}\rangle$ is the spatial average of the Rosseland opacity over $\lambda$. A relevant tracer of the transition $t_{\mathrm{thin}}=t_{\mathrm{thick}}$ is the Rosseland optical depth over the region of the same sign of the perturbation amplitude (between two subsequent nodes)444For dust continuum, $\tau_{\mathrm{R}}\approx 1$ at the transition but for pure gas this value is lower due to the larger difference between the Planck and the Rosseland mean opacities.,

[TABLE]

Figure 3 presents a verification of the maximum time criterion introduced in Sect. 2. In all panels of Fig. 3 we have marked $\tau_{\mathrm{R,\lambda/2}}=1$, which closely (but not exactly) traces the transition between the optically thick and the optically thin relaxation (2). In an opaque inner disc, even the high-frequency modes are optically thick and relax on the radiative diffusion timescale. This can be either dynamically slow, $\Omega t>1$ (up to $\Omega t\sim 10^{3}$), or dynamically fast, $\Omega t<1$ (down to $\Omega t\sim 10^{-3}$). The fastest achievable radiative relaxation, $\Omega t\sim 10^{-4}-10^{-3}$, is that of the optically thin modes in the yet coupled regions ($t_{\mathrm{coll}}<t_{\mathrm{emit}}$, see Fig. 3).

Absorption of the stellar photons ($T_{*}=5777\mbox{ K}$) is shown in Fig. 3 for the two gas opacity models. In the small gas opacity case, the single-temperature Planck means from Freedman et al. (2008) are used to calculate the optical depth:

[TABLE]

In the large gas opacity case, the two-temperature555The temperature of the radiation field $B_{\nu}\left(T_{*}\right)$ differs from the gas temperature $T_{\mathrm{g}}$ Planck means are available from Malygin et al. (2014):

[TABLE]

In both cases, the integration is done along the ray of constant polar angle in spherical coordinates between $0{.}3\mbox{ au}$ and $300\mbox{ au}$. The contribution from innermost regions ($<0{.}3\mbox{ au}$) is neglected.

Figure 4 displays a comparison between LTE emission timescale (4), collision timescale (5), and radiative diffusion timescale (9) for midplane modes and the fiducial dust-to-gas mass ratio of $0{.}014$. The two panels correspond to the two gas opacity models. The optically thick relaxation, being dominated by the dust, is not sensitive to the gas opacity model. Relaxation timescale in this regime drops with radius, its absolute value depends on the perturbation wavelength. Low-frequency modes $\tilde{\lambda}=1$ relax dynamically slow ($1<\Omega t<10^{3}$) within first $\approx 8\mbox{ au}$. The farther out the faster they relax: the relaxation timescale drops to $\approx 10^{-3}\Omega^{-1}$ at $\approx 50\mbox{ au}$. High-frequency yet optically thick modes $\tilde{\lambda}\leq 10^{-2}$ always relax dynamically fast ($\Omega t<1$). Optically thick relaxation timescale (9) significantly underestimates the relaxation timescale if a mode becomes optically thin (the long-dashed lines in Fig. 4). This is an artefact of the free-streaming limit in the flux-limited diffusion approximation (see Sect. 5).

The gas opacity model becomes important for optically thin modes. In the large gas opacity case (Fig. 4b), the optically thin midplane modes relax on the gas thermal emission timescale, $\Omega t\approx 2{.}5\times 10^{-4}$, from $15\mbox{ au}$ to $200\mbox{ au}$ in the midplane.

Figure 5 maps the timescales of the processes that contribute to the thermal relaxation all over the disc. The map of the radiative diffusion timescale (Fig. 5a) is plotted for $\tilde{\lambda}=2$. Figure 5b maps the optically thin relaxation timescale. This is a combination of the collision and the thermal emission timescales according to (7). Figures 5c and d also show the dust-to-gas collision timescale and the dust LTE emission timescale, respectively. The assumed vertically unstratified temperature structure is captured in Fig. 5d. The size distribution of the grains (see Appendix B.2) as well as the dust-to-gas mass ratio are uniform over the disc. A wedge-like structure at around $3\mbox{ au}$ in Figs. 5b,d is due to the water ice line.

4.2 Perturbation wavelength

At a given moment of time, the density, temperature, and opacities are defined at any point in the disc. Then the regime of radiative relaxation – optically thick or thin – is decided by the perturbation wavelength.

Figure 6 shows relaxation timescale of vertical modes as a function of their wavelength at different radial locations. Relaxation via thermal emission puts a wavelength-independent lower limit on relaxation timescale (indicated in Fig. 6). In the optically thick regime, the dependence of the relaxation timescale on $\lambda$ changes itself with $\lambda$: for $\tilde{\lambda}\la 1$,

[TABLE]

and for $\tilde{\lambda}>1$,

[TABLE]

because the Rosseland depth saturates (see the dotted lines in Fig. 6).

In the bulk of the disc, the relaxation of optically thin modes via thermal emission is dynamically fast ($\Omega t\ll 1$, see Figs. 6b-d), optically thick modes can be both isothermal ($\Omega t\ll 1$) and adiabatic ($\Omega t\gg 1$). High-frequency modes $\tilde{\lambda}\ll 10^{-1}$ remain isothermal even if they are optically thick (Fig. 6a, Fig. 4). A transition between the two relaxation regimes occurs at larger radii for larger $\tilde{\lambda}$ such that the Rosseland depth of the mode $\tau_{\mathrm{R}}\approx 2$ at the transition (Figs. 6b,c, see also Sect. 4.1).

4.3 Anisotropic relaxation

Optically thick relaxation in non-isotropic media can be anisotropic itself: diffusion timescale in different directions is, in general, different. Figure 7 maps a relative difference between relaxation timescales of vertical and radial modes of the same wavelength (anisotropy, hereafter). Vertical diffusion is up to $60\%$ faster than the radial diffusion for $\tilde{\lambda}=2$ modes in the midplane (Fig. 7a). This is opposite above the midplane: diffusion in the radial direction is faster than in the vertical, and the difference is larger. The anisotropy is the largest when the radial direction becomes optically thin while the vertical direction remains optically thick (due to the exponential increase of the density towards the midplane). This is typically at $2-3H$ above the midplane. The maximum anisotropy is $\approx 100\%$ for $\tilde{\lambda}=2$, $\approx 1\%$ for $\tilde{\lambda}=0{.}1$ and $\approx 0{.}5\%$ for $\tilde{\lambda}=0{.}02$.

4.4 Dust-to-gas mass ratio

Models of Birnstiel et al. (2012) suggest that the dust growth during planet formation leads to a significant dust opacity decrease. In this section, we show how the radiative relaxation time changes with dust depletion degree. We homogeneously reduce the dust-to-gas mass ratio while leaving both the effective collision cross section per particle (the size distribution) and the dust opacity per unit dust mass (composition) unchanged for simplicity. Dust thermal radiation timescale (4) and dust-to-gas collision timescale (5) both increase with dust reduction. This is also expected to be so due to depletion of small grains caused by the dust growth.

Under typical PPD conditions, the Rosseland mean is dominated by the dust continuum below some $1\,500\mbox{ K}$, but the Planck means of the gas and the dust compare down to at least $650\mbox{ K}$ owing to molecular ro-vibrational bands (Malygin et al. 2014). The optically thick relaxation is determined by Rosseland opacity and hence is sensitive to the dust evolution effects. The gas opacity is negligible in the optically thick relaxation up until significant ($\times 10^{-2}$) dust depletion. The optically thin regime is, in contrast, sensitive to the gas opacity model, and already at the fiducial mass ratio (Fig. 4).

The optically thick relaxation in Fig. 8 is calculated for radial modes $\tilde{\lambda}_{r}=0{.}628$ ($\tilde{k}_{r}=10$, with an application for the VSI, see Sect. 4.6.1), the qualitative trend is similar for other wavelengths. Dust amount reduction yields speed-up in optically thick relaxation. With fewer absorbers, the Rosseland optical depth is smaller, hence, the radiative diffusion timescale is shorter. The relaxation timescale drops with dust reduction until it equals the optically thin relaxation timescale.

In what way the optically thin relaxation changes with dust depletion, depends on how the Planck means of the dust and the gas compare. In the small gas opacity case, the optically thin relaxation slows down with the dust amount decrease (Fig. 8a). This is due to the shrinkage of the net dust emitting surface. In the case of zero dust (the dotted line in Fig. 8a), the disc is optically thin in the vertical direction and the relaxation is restrained by the gas emissivity. In the large gas opacity case, however, optically thin modes do not “feel” the dust depletion. This is because the gas determines the relaxation as a more efficient (larger Planck mean) yet collisionally faster-coupled emitter than the dust already at the fiducial density ratio of $0{.}014$.

Figures 9 and 10 map the relaxation timescale of $\tilde{k}_{r}=10$ modes (see Sect. 4.6.1) under different dust depletion degrees for the small and the large gas opacity case, respectively. The dust-to-gas mass ratio is varied evenly throughout the disc. Therefore, the maps do not describe the relaxation in a realistic, planet-forming environment, where the ratio can be unevenly distributed and evolve with time. However, these highlight the dependence of the relaxation on gas opacity and dust depletion, and indicate the necessity of taking into account the mean gas opacity in numerical modelling with dust evolution.

4.5 Disc mass

The outcome of varying disc mass is shown in relaxation maps in Fig. 11. The effect is due to the density change. The optically thick relaxation is slower in a more massive disc than in a lighter one because of the larger Rosseland depth (17) across the mode in the denser medium. This is responsible for the change in relaxation in the inner opaque parts of the disc: bottom left area in Figs. 11a-d. The collision timescales vary with density as $n^{-1}$, so the density change affects relaxation in collisionally-decoupled outer disc: top right area in Figs. 11a-d. The thermal emission timescale does not change with the disc mass because it depends only weakly on density (through the Planck mean opacity). The net result is: in more massive discs, layers of fast relaxation are narrower and located farther out from the central object, comparing to lighter discs.

4.6 Hydrodynamic transport in circumstellar discs

To generate turbulence in an axisymmetric and vertically unstratified disc flow, a purely hydrodynamic mechanism must circumvent the strong centrifugal stability imposed by a positive radial specific angular momentum gradient (Rayleigh’s stability criterion Balbus et al. 1996). A general requirement for the dynamic stability of a stratified fluid, the Solberg-Høiland conditions (Rüdiger et al. 2002), allow development of an instability for Rayleigh-stable (if unstratified) velocity fields. The circumvention of the stabilising radial stratification is either via vertical stratification of angular frequency $\Omega(z)$ or negative entropy gradient $\partial S/\partial r<0$. Each of the identified instabilities described below requires conditions that can be used to constrain its location within a disc (see Fig. 12).

4.6.1 Vertical shear instability

The vertical shear instability (VSI) can develop in a disc with a global temperature profile, where the angular frequency is a function of both radial and vertical coordinates (Urpin & Brandenburg 1998; Umurhan et al. 2013). The angular momentum transport by the VSI can be viewed as such a swap of the fluid elements along the curved iso-line of the specific angular momentum, which yields a state with lower total energy (Barker & Latter 2015). The weak vertical stratification in the orbital motion ($\Omega\left(z\right)$) permits circumvention of the centrifugal stability via driving vertically elongated (and slightly tilted) but radially narrow flows. Vertically elongated disturbances $(k_{r}/k_{z}\gg 1)$ can overcome the Coriolis stabilisation. However, motions of this kind are impeded by vertical buoyancy in discs with a stable vertical stratification. The VSI requires sufficiently fast thermal relaxation to overcome the stabilising vertical buoyancy (Urpin 2003; Lin & Youdin 2015). The angular momentum transfer rates resulting from the VSI are modest: $\alpha\approx 10^{-4}$ (as found so far in numerical modelling by Arlt & Urpin 2004; Nelson et al. 2013; Stoll & Kley 2014).

Nelson et al. (2013) performed a detailed numerical study of the VSI for two classes of equilibrium disc structures with controversial results on the cooling criterion. The vertically isothermal models indeed demonstrated the requirement of very rapid thermal relaxation, $\Omega t\la 10^{-2}$, for the instability to develop. In contrast, the vertical polytropic models pushed the critical relaxation limit as far as to $\Omega t\approx 10$ for specific temperature and density profiles. In analytic investigation of the VSI including the radiative cooling, Lin & Youdin (2015) confirmed a generic utility of the critical cooling timescale of the form

[TABLE]

with $h$ being the disc aspect ratio, $q$ the radial power-law temperature gradient, and $\gamma$ the adiabatic index. Lin & Youdin (2015) showed that time (22) is a robust critical relaxation criterion for modes $5\leq\tilde{k}_{r}\leq 10$. For low-frequency modes $\tilde{k}_{r}=1$, the threshold is at slightly longer times. For high-frequency modes $\tilde{k}_{r}\ga 30$, the linear growth carries on for cooling times considerably longer than (22), but the growth rate drops significantly (see Fig.13 in Lin & Youdin 2015). In Fig. 12, the zone of onset of VSI modes $\tilde{\lambda}_{r}=0{.}628$, predicted by $t<t_{\mathrm{c}}=6{.}25\times 10^{-2}\Omega^{-1}$, is hashed with grey lines. This occupies a sizeable fraction of the disc volume, and spans from $\approx 15\mbox{ au}$ to $\approx 180\mbox{ au}$ in the midplane. Above the midplane, the VSI can develop closer to the star.

4.6.2 Convective overstability

Jumps in opacity, surface density, ionisation, or stellar irradiation at close radii can introduce a negative radial entropy gradient, at least in some specific locations in a disc. This negative entropy gradient makes the protoplanetary discs baroclinic666Baroclinic means having misaligned isobars and isopycnals, which provides another route to instability (Klahr & Bodenheimer 2001; Raettig et al. 2013). Recently, (Klahr & Hubbard 2014) has shown that radially stratified baroclinic flows are linearly unstable to non-axisymmetric perturbations given a finite thermal relaxation time. During an oscillatory epicyclic motion, the parcel of gas heats up when displaced radially inwards, and cools down when displaced radially outwards. When crossing its starting position, the gas has a different temperature from its surroundings if the relaxation time is non-zero. The related density difference enhances the restoring buoyant acceleration (overstability), and the initial oscillation is amplified. This phenomenon obtained the name ‘convective overstability’ (Klahr & Hubbard 2014). The finite-amplitude noise resultant from the convective overstability (COV) can trigger a non-linear state of vortex amplification coined as subcritical baroclinic instability (SBI, Klahr & Bodenheimer 2003; Petersen et al. 2007; Lesur & Papaloizou 2010; Lyra & Klahr 2011). Lyra (2014) confirmed the phenomenon via compressible linear analysis and demonstrated in both 2D and 3D shearing box runs that the SBI is indeed the saturated state of the convective overstability. The lifetime of the resultant vortex is also a non-monotonic function of thermal relaxation rate. The exact nonlinear outcome of the COV in discs with radial and vertical stratification is yet to be determined in numerical modelling (Latter 2016). The SBI eventually yields growing vortices, which can induce Reynolds stresses rateable to $\alpha=10^{-4}-10^{-2}$ (Raettig et al. 2013, 2015; Lobo Gomes et al. 2015). These anticyclonic vortices can also be triggers for fast and efficient planetesimal formation (Klahr 2003).

The growth rate of linear $r-z$ mode of the COV peaks at $\Omega t\gamma=1$ (where $\gamma$ is the gas adiabatic index, Klahr & Hubbard 2014). Yet the growth rate of the non-linear regime, the self-similar amplification of vortices, peaks more likely at $\Omega t/\chi=1$ ($\chi$ being the aspect ratio of the disc). Relaxation rates of $10^{-2}<\Omega t<10^{2}$ are gained by $\tilde{\lambda}_{z}=2$ modes in the midplane from $\approx 2\mbox{ au}$ to $\approx 40\mbox{ au}$, and beyond $\approx 130\mbox{ au}$. Away from the midplane, the suitable relaxation zones bend towards the central object (the light grey areas in Fig. 12) forming two arc-shaped regions, where the COV can develop. The contour lines of $\Omega t=1$ are approximately in the middle of each of those regions. For higher-frequency modes, $\tilde{\lambda}<2$, the inner arc region is even closer to the star: its inner edge is at $\approx 0{.}6\mbox{ au}$ for $\tilde{\lambda}=0{.}5$.

4.6.3 Zombie vortices

Marcus et al. (2015, 2016) showed that in 3D stably stratified baroclinic (Rayleigh-stable if unstratified) adiabatic shearing flows, an initial Kolmogorov turbulence undergoes a reverse energy cascade due to the background shear. In an environment with negative rotation shear (background Keplerian shear), anticyclonic regions merge to form a bigger anticyclonic vortex, while cyclonic regions get stretched into thin cyclonic vortex layers. At some point, wave solutions for vorticity become singular because of the shear, and a baroclinic critical layer is formed (Maslowe 1986). This baroclinic critical layer then induces two vortex layers with the opposite sense of rotation. Anticyclonic layers are linearly unstable, roll-up into individual vortices, which then merge into one larger vortex (Marcus et al. 2013). The resultant vortex excites inertial gravity waves, which propagation is affected by the Coriolis force and the shear. The instability is the result of a resonant interaction between a Rossby wave and a gravity wave (Umurhan et al. 2016). Each generation of vortices excites new baroclinic critical layers, and the process repeats. The non-linear phase of self-reproduction of such structures obtained a name of zombie mode instability or zombie vortex instability (ZVI).

The estimates of the angular momentum transport rate from the ZVI are not yet conclusive due to small local domains considered in numerical studies. The speculations are that the acoustic waves shed by zombie vortices can provide $\alpha\sim 10^{-3}$ (Marcus et al. 2015). Also, the ZVI could contribute to turbulent mixing (Turner et al. 2007).

The inertial gravity waves, being driven by buoyancy, require long thermal relaxation times, $N^{2}t^{2}\gg 1$ ($N$ is the Brunt-Väisälä frequency). Lesur & Latter (2016) found that the ZVI requires $\Omega t>16$ and simultaneously Peclet numbers777a dimensionless number comparing nonlinear advection to thermal diffusion of $\mbox{Pe}>10^{4}$. This second condition constrains the volume, where the ZVI can operate, to the regions within $0{.}8\mbox{ au}$ in the radial direction and one local pressure scale height around the midplane (grey area in Fig. 12). These zone is indeed MRI inactive according to Dzyurkevich et al. (2013). Thus, the ZVI is a viable model for driving turbulence in inner regions of protoplanetary discs.

5 Discussion

Hydrodynamic means of angular momentum transport in discs need to be studied because they can provide a turbulent stress in MRI-inactive regions and aid planetesimal formation. In numerical modelling for self-gravitating discs (Johnson & Gammie 2003; Boley et al. 2006), a usable estimate of the disc cooling results from vertically integrated azimuthally-averaged thermal energy balance (see, however, Takahashi et al. 2016). In other identified hydrodynamic sources of disc turbulence, a resolution of the vertical direction and the high-frequency part of the perturbation spectrum is important (Nelson et al. 2013; Klahr & Hubbard 2014). For these, the cooling timescale of a disc annulus is not a viable estimate of the thermal relaxation (Nero & Bjorkman 2009; Nero 2010). The relaxation timescale refers to a particular perturbation mode, and, in general, depends on its wavelength $\lambda$, location and orientation.

In this work, we map radiative relaxation timescale of linear temperature perturbations of various wavelengths over disc interiors of given $\left(\rho,T\right)$ structures. We distinguish between two different regimes of relaxation: the optically thick and the optically thin (see Fig. 4). Optically thick (in the Rosseland sense) perturbations relax on the radiative diffusion timescale, $t_{\mathrm{thick}}\propto\rho\tau_{\mathrm{R}}\lambda T^{-3}$. The optically thick relaxation is, in general, anisotropic. In the midplane, optically thick vertical modes $\lambda_{z}=2H$ ($H$ is the pressure scale height) relax up to $60\%$ faster than radial modes $\lambda_{r}=2H$. Above the midplane, the radial relaxation can be up to $100\%$ faster. Optically thin perturbations relax on the thermal emission timescale, $t_{\mathrm{thin}}\propto\kappa_{\mathrm{P}}^{-1}T^{-3}$, provided LTE. In our calculations, we account for the collisional decoupling from the emitters in the optically thin relaxation. Because of this, our relaxation timescales in low-density atmospheric layers and cold distant midplane regions are longer than those found by Lin & Youdin (2015). However, in the case of non-LTE, other radiative processes, not included in this study, should be taken into account, so the collisional timescale is, in general, not a reliable estimate of the relaxation timescale.

The two relaxation regimes depend differently on relevant opacity: $t_{\mathrm{thick}}\propto\kappa_{\mathrm{R}}$, $t_{\mathrm{thin}}\propto\kappa_{\mathrm{P}}^{-1}$. The antipodal relation to opacity results in the opposite response to the reduction of dust amount (Cai et al. 2006). We do not consider here any time evolution of the dust, which is required to accurately determine the relaxation timescales in the inner disc, as pointed out in Nelson et al. (2000). Our midplane temperatures, though, never exceed those of severe dust depletion, yet evaporation and condensation of the species is taken into account in the dust opacity tables (Semenov et al. 2003).

We show that the diffusion timescales calculated in the flux-limited diffusion approximation (as in Stoll & Kley 2014) underestimate the relaxation timescale in the optically thin medium (see Fig. 4). This is because the flux-limited diffusion formalism (28) does not account for the finite timescale of the intrinsic radiative processes but causal restriction on the energy density propagation (the free-streaming limit). A finite timescale of the collisional coupling (5) and finite ratio of the material’s thermal capacity to its emissivity (4) both have to be taken into account. We report the bottom limit on the radiative relaxation of $\Omega t\approx 10^{-4}$ (for the considered models) set by thermal emission timescale (4).

5.1 The method: thin and thick

In this section, we discuss the limits of the current method and highlight the measures for further improvement.

In linear analysis, relaxation timescale (2) is independent of the perturbation amplitude $\delta T_{\mathrm{0}}$ and is solely defined by the location within the disc. A value of the relaxation timescale can be assigned to each point $x$ of an optically thin Fourier mode $\delta T(\@vec{x})$. The maximum of these values can be taken as the relaxation timescale of the mode. We make a simplification by taking the relaxation time at the point of extremal $\delta T=\delta T_{\mathrm{min/max}}=\delta T_{\mathrm{m}}$. This is a viable approximation given smooth Planck opacity within the perturbation wavelength.

The optically thick relaxation timescale demands non-local calculations: the photons emitted at the point of $\delta T_{\mathrm{0}}$ subsequently diffuse through the medium with $0<|\delta T/\delta T_{\mathrm{m}}|<1$. Variations in the diffusion length scale at all those points alter the relaxation timescale. We use an effective diffusion coefficient harmonically averaged over the perturbation wavelength $\lambda$ (see Appendix A.2). For an arbitrary optically thick perturbation $\left(k_{z},k_{r},m\right)$, the bulk of the excess thermal energy escapes through the direction of the least Rosseland depth (cf. the flashlight effect, Nakano 1989). The optically thick relaxation is, therefore, sensitive to the disc density structure and especially its non-isotropy or time fluctuations. Numerical modelling with time evolution is required here.

The radiative diffusion and thermal emission timescales both imply fast – on timescales shorter than the relaxation timescale itself – thermal coupling between the radiating (dust and gas) species and the thermal carriers (the gas). If thermal emission timescale (4) is shorter than the time to set LTE, we limit the relaxation by the appropriate collision timescale. Non-LTE cooling and non-equilibrium thermal feedback processes between the dust, the gas and the radiation field (external irradiation) are the necessary further refinements.

The optical depths and emission rates are calculated using the joint dust and gas mean opacity tables obtained by direct summation of the corresponding mean values. At low temperatures ($\la 1\,500\mbox{ K}$), the Rosseland opacity of the dust exceeds that of the solar-mixture gas by several orders of magnitude. Hence, adding the two means to each other does not introduce a considerable error far from the dust sublimation temperature: the opacity is either dominated by the dust or by the gas. Further, data from Semenov et al. (2003) accounts for the freeze out of volatile organics, water and carbon monoxide, the main gas-phase absorbers at low temperatures. In the case of considerable (by a factor of $\approx 10^{2}$ or more) dust depletion, this approach is invalid because the Rosseland averaging is not additive.

Dust opacity depends not only on the dust chemistry but on the grain size distribution as well (Pollack et al. 1994; Henning & Stognienko 1996). As pointed out in Nelson et al. (2000), the grain size evolution proceeds on short (relative to the local orbital period) timescales in the inner disc. Calculating the radiative relaxation there necessitates considering time-dependent size and composition of the grains. In the atmospheric layers imposed to an external radiation field, the photochemistry breaks the assumptions of equilibrium chemistry and LTE, under which the opacities were calculated. The resolution of this issue requires running a frequency-dependent radiative transport on top of a non-equilibrium chemical reactions network and disc hydrodynamics. However, gas chemistry in dusty environments of discs is poorly constrained: many important reaction rates and transitions have not yet been measured in laboratories (Dutrey et al. 2014). Our two constructed gas opacity models probe the parameter space. The analysis could benefit from the yet missing consistent dust and gas opacity calculations.

6 Summary

The radiative cooling timescale is important in the context of disc stability, associated transport, and disc observational appearance. A simple recipe to compute the radiative relaxation timescale in protoplanetary discs is introduced. The relaxation timescale as a function of perturbation wavelength, temperature and density is obtained by means of non-compressible linear analysis. There are two regimes of relaxation depending on how the thermal energy evolves: the optically thick and the optically thin. First is determined by radiative diffusion, second – by thermal emission. In comparison to previous studies, collisional decoupling is taken into account in the optically thin relaxation, both the Planck and the Rosseland tabulated opacities of the dust and the gas are used, different dust depletion degrees are considered. The optically thick regime can be both dynamically fast and slow ($\Omega t\leogr 1$), the optically thin one is dynamically fast if not limited by collisional coupling. The fastest, LTE relaxation is $\Omega t\sim 10^{-4}$ in a typical T Tau disc.

With this recipe, one can evaluate the relaxation timescale of a grid cell or an SPH particle based on its actual optical depth. This has applications for studying hydrodynamic and gravitational instabilities in the context of planet formation and transport in accretion discs, core migration, and planet-disc interaction. Calculated maps of the relaxation timescale constrain the origins of the identified hydrodynamic instabilities in discs (Fig. 12). A large volume of a typical T Tau disc is unstable to the development of linear hydrodynamic instabilities: the convective overstability (COV, Klahr & Hubbard 2014) and the vertical shear (a.k.a. Goldreich-Schubert-Fricke) instability (VSI, Nelson et al. 2013). The zone of operation of the VSI spans from $15\mbox{ au}$ to $180\mbox{ au}$ in the midplane. The COV can develop both closer in and farther out than that but there is a COV-dead layer from $40\mbox{ au}$ to $140\mbox{ au}$ in the midplane, where relaxation is too fast ($\Omega t<10^{2}$) even for the low-frequency modes $\lambda=2H$. A subcritical hydrodynamic instability, the zombie vortex replication (Marcus et al. 2013), can work inside $\approx 0{.}8\mbox{ au}$ and within one pressure scale height, where thermal ionisation is yet inefficient to drive the MRI. Hydrodynamic modelling with time-evolution of the opacities is required for further study.

Acknowledgements.

This research has been supported by the International Max-Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). M.M. thanks Hubert Klahr for valuable insights on hydrodynamic stability of rotating fluids, Kai-Martin Dittkrist for fruitful discussions and hints. The authors thank the anonymous referee for useful comments on the scope of the presented approach and various suggestions, which significantly improved the appearance of the material.

Appendix A Optically thick relaxation

A.1 Basic equations

To obtain a characteristic damping time of optically thick thermal perturbations, we start with the thermal energy evolution equation considering only radiative processes (see Kuiper et al. 2010)

[TABLE]

Here $E_{\mathrm{int}}$ is the internal energy density of the gas, $E_{\mathrm{R}}$ the radiation energy density, $\@vec{F}_{\mathrm{R}}$ the flux of the radiation energy density. At high optical depths, the temperature of the gas and the ambient thermal radiation field are equal. The internal energy density of the gas with mass density $\rho$ and temperature $T$ is given by

[TABLE]

with $C_{\mathrm{v}}$ being the specific thermal capacity at constant volume. The radiation energy density is defined by

[TABLE]

where $a=4\sigma/c$ is the radiation constant. Use of (24) and (25) in the left-hand side of (23) yields

[TABLE]

Here, $\eta=E_{\mathrm{R}}/(E_{\mathrm{R}}+E_{\mathrm{int}})$. The flux of the thermal radiation energy reads

[TABLE]

In the framework of the flux-limited diffusion, the diffusion coefficient $D$ is

[TABLE]

with $\kappa_{\mathrm{R}}$ being the Rosseland mean opacity, $c$ the speed of light, $\rho$ the mass density, and $\varkappa$ the flux limiter (Levermore & Pomraning 1981). The gradient of the thermal radiation flux unveils

[TABLE]

with

[TABLE]

Equilibrium temperature $T_{\mathrm{0}}$ satisfies

[TABLE]

Using (26) and (29) in (23), one finds the temperature evolution equation

[TABLE]

with

[TABLE]

following Kuiper et al. (2010).

A.2 Linear analysis

Linearisation $T=T_{\mathrm{0}}+\delta T$ of (32) with the use of (31) gives

[TABLE]

Let

[TABLE]

be a wavenumber of the Fourier mode with spatial alignment along direction $x$ of the fastest radiative heat transport (not necessarily the vertical direction). A dissipation time of this mode is then found by employing the Fourier ansatz $\delta T_{k}=\delta T_{k}^{0}\exp\left(-i\@vec{k}\@vec{x}+i\omega t\right)+c.c.$:

[TABLE]

Here $\@vec{A}^{\prime}=\partial_{T}\left[\nabla\ln\rho\left(1+\rho\partial_{\rho}\ln\kappa_{\mathrm{R}}\right)\right]$ is a partial derivative of (30) with respect to temperature, and

[TABLE]

the effective diffusion coefficient. If the wavenumber is small (or the density fluctuations within it are significant: $\left|\nabla_{z}\rho\right|\rho^{-1}>k$), we perform averaging over $\lambda$

[TABLE]

Because the Rosseland mean for dust is a moderately varying function of density, $|\partial_{\rho}\kappa_{\mathrm{R}}|\ll\rho^{-1}\kappa_{\mathrm{R}}$, we approximate $\@vec{A}^{\prime}=\partial_{T}\nabla\ln\rho$. For an adiabatic equation of state, $\rho\propto T^{\gamma}$, $\@vec{A}^{\prime}\nabla T=-\gamma\left(\nabla\ln T\right)^{2}+\nabla\gamma\nabla\ln T$. Separating the terms proportional to first and second powers of $\nabla\ln T$ in the denominator, one finds

[TABLE]

Here, $k_{\mathrm{T}}=\left(3+\gamma\right)^{1/2}\nabla\ln T$ is the equilibrium temperature wavenumber, and $k_{\mathrm{A}}=\left(\nabla\gamma\nabla\ln T\right)^{1/2}$ the adiabatic wavenumber. In a hydrostatic picture, we argue, the factor $\nabla\gamma$ makes $k_{\mathrm{A}}$ negligible in favour of the thermal wavenumber $k_{\mathrm{T}}$. The equation reduces to

[TABLE]

The thermal wavenumber $k_{\mathrm{T}}$ is determined by the characteristic length scale of equilibrium temperature gradient. In a vertically isothermal radially stratified disc with $T\propto r^{-q}$ this term may be estimated as $\nabla T/T\sim-qh\mbox{H}^{-1}$, with H being the local pressure scale height, and $h$ the local aspect ratio usually of the order of $0{.}05$. Thus, we set

[TABLE]

For high-frequency modes, $\lambda\ll\mbox{H}$, the equilibrium thermal wavenumber is negligible in favour of the wavenumber of the mode itself. For low-frequency modes, $\lambda\ga\mbox{H}$, (41) still holds true for geometrically thin discs (small aspect ratios). Hence, the optically thick relaxation time evaluates the diffusion time of the excess thermal photons over the mode wavelength. For the application cases of vertically isothermal disc, it is legitimate to permanently set

[TABLE]

In this case, the relaxation time of an optically thick mode depends on its wavenumber $k$, the radiative diffusion coefficient $D$, and the ratio $f$, which is smaller than unity in this regime (see Sect. A.3).

A.3 Energy or temperature evolution time?

Dividing (26) by the total thermal energy $E_{\mathrm{tot}}=E_{\mathrm{R}}+E_{\mathrm{int}}$ gives a relation between the energy and temperature evolution times

[TABLE]

If $\eta\ll 1$, the two times are equal up to $O(3\eta)$. The “critical” density $\rho_{\mathrm{cr}}=aT^{3}/c_{\mathrm{V}}$ at which the radiation energy density equals the internal energy density is as low as

[TABLE]

For instance, $\rho_{\mathrm{cr}}\approx 10^{-15}\mbox{ g cm}^{-3}$ at $300\mbox{ K}$, that is several orders of magnitude below the typical midplane values for a PPD in T Tau phase. We find $\eta$ between $10^{-5}$ and $10^{-3}$ within the optically thick regions in the fiducial disc model.

Appendix B Optically thin relaxation

There are two important timescales in the optically thick relaxation: the timescale of the collisional coupling to the emitters, and the LTE emission timescale (see Sect. 2.1).

B.1 Thermal emission time

Consider the following energy balance equation:

[TABLE]

Here, the terms $Q^{+/-}$ designate the specific (per unit mass) heating or cooling terms due to the perturbation on top of an equilibrium temperature field. The emission rate of matter at equilibrium temperature $T$ is proportional to the Planck mean opacity $\kappa_{\mathrm{P}}$ at that temperature (Kirchhoff’s law). The perturbed internal energy then evolves as

[TABLE]

where $\delta T$ is the amplitude of the temperature perturbation. Substituting a Fourier mode $(T-T_{0})_{k}=\delta T_{k}=\delta T_{k}^{0}\exp\left(-i\@vec{k}\@vec{x}-i\omega t\right)$, we linearise

[TABLE]

This implies the evolutionary time

[TABLE]

independent of the perturbation wavenumber.

B.2 Collisional coupling

In PPDs, relative velocities in both dust-to-gas and gas-to-gas collisions are attributed to the thermal motions of the gas particles

[TABLE]

where $\mu$ is the mean molecular weight, $k_{\mathrm{B}}$ the Boltzmann constant. The mean time between two subsequent collisions is

[TABLE]

where $\sigma_{\mathrm{c}}$ is the collisional cross section with the targets of number density $n$.

B.2.1 Dust-to-gas collisions

The underlying grain size distribution in Semenov et al. (2003) is an MRN function (Mathis et al. 1977) with a spectral slope of $3{.}5$ and the modified maximum grain size of $5\mbox{ }\mu\mbox{m}$ (Pollack et al. 1985):

[TABLE]

This results in the effective collisional cross section $\sigma\approx\pi a_{\mathrm{+}}a_{\mathrm{-}}\approx 1{.}5\times 10^{-9}\mbox{ cm}^{2}$.

B.2.2 Gas-to-gas collisions at low temperatures

The emissivity of the gas is supplied by different species at different temperatures. In LTE, the most efficiently radiating species can be found by comparing their relative contributions to the total Planck mean

[TABLE]

Here $i$ labels the species (e.g. CO, H2O, TiO, NH3, etc.), and

[TABLE]

is a contribution to the Planck mean from specie $i$. The effective coupling time to radiative gas species can be estimated as

[TABLE]

At the temperatures of interest, the most abundant species, H2 and He, are not efficient emitters in terms of (53), but they store the bulk of thermal energy. The rest of the matter can be modelled as a sort of efficient emitters satisfying

[TABLE]

Coupling timescale (54) is then approximated as collisional timescale

[TABLE]

with the effective neutral-neutral collisional cross section $\sigma_{\mathrm{nn}}$. We estimate

[TABLE]

where $f=n^{\prime}/n_{\mathrm{gas}}$ is the number fraction of the radiating species. In our calculations, we put $f=10^{-2}$ and $\sigma_{\mathrm{nn}}=8{.}5\times 10^{-17}\mbox{ cm}^{2}$.