Gravito-turbulence in irradiated protoplanetary discs

TL;DR

This study uses radiation hydrodynamics simulations to investigate gravito-turbulence in irradiated protoplanetary discs at 50 AU, revealing a finite surface density range for turbulence, the nature of turbulent motions, and the effects of irradiation and cooling.

Contribution

First detailed simulation study of gravito-turbulence in irradiated protoplanetary discs with realistic physics, identifying the turbulence range and energy transport mechanisms.

Findings

Gravito-turbulence exists between 80 and 250 gcm$^{-2}$ surface density.

Turbulent motions are supersonic at all heights, dissipating via shocks and compression.

Irradiation has minimal impact on midplane turbulence unless the grazing angle exceeds 0.32.

Abstract

Using radiation hydrodynamics simulations in a local stratified shearing box with realistic equations of state and opacities, we explored the outcome of self-gravity at 50 AU in a protoplanetary disc irradiated by the central star. We found that gravito-turbulence is sustained for a finite range of the surface density, from $\sim 80$ to $\sim$ 250 gcm$^{-2}$. The disk is laminar below the range while fragments above it. In the range of gravito-turbulence, the Toomre parameter decreases monotonically from $\sim 1$ to $\sim 0.7$ as the surface density increases while an effective cooling time is almost constant at $\sim 4$ in terms of the inverse of the orbital frequency. The turbulent motions are supersonic at all heights, which dissipates through both shock waves and compressional heating. The compressional motions, occurring near the midplane, create upward flows, which not only contribute to supporting the disc but also to transporting the dissipated energy to the disc surfaces. The irradiation does not affect much the gravito-turbulence near the midplane unless the grazing angle is larger than 0.32. We also show that a simple cooling function with a constant cooling time does not approximate the realistic cooling.

Equations

$$\beta \equiv t_{\text{cool}}\Omega <3,$$

$$Q\equiv \frac{c_{\text{s}}\Omega }{\pi G\Sigma },$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot (\rho \mathbf{v}\mathbf{v})=-\nabla p-\rho \nabla \Phi +\frac{{\kappa }_{\text{R}}\rho }{c}\mathbf{F},$$

$$\frac{\partial e}{\partial t}+\nabla \cdot (e\mathbf{v})=-(\nabla \cdot \mathbf{v})p-\left(4\pi B(T)-cE\right){\kappa }_{\text{P}}\rho ,$$

$$\frac{\partial E}{\partial t}+\nabla \cdot (E\mathbf{v})=-\nabla \mathbf{v}:\mathsf{P}+\left(4\pi B(T)-cE\right){\kappa }_{\text{P}}\rho -\nabla \cdot \mathbf{F},$$

$${\nabla }^2\Phi =4\pi G\rho ,$$

$$\frac{\partial e}{\partial t}=-{\kappa }_{\text{P}}\rho \left(4\pi B(T)-cE\right)+q_{\text{irr}},$$

$$\frac{\partial E}{\partial t}={\kappa }_{\text{P}}\rho \left(4\pi B(T)-cE\right)+\nabla \cdot \left(\frac{c\lambda (R)}{{\kappa }_{\text{R}}\rho }\nabla E\right),$$

$$0=-{\kappa }_{\text{P∗}}\rho I-\frac{dI}{ds},$$

$$⟨f⟩(z,t)\equiv \frac{\int \int f(x,y,z,t)dxdy}{\int \int dxdy},$$

$$\overline{⟨f⟩}\equiv \frac{\int ⟨f⟩dt}{\int dt},$$

$$⟨⟨f⟩⟩\equiv \int ⟨f⟩dz.$$

$${⟨⟨f⟩⟩}_{\rho }\equiv \frac{\int ⟨f⟩⟨\rho ⟩dz}{\int ⟨\rho ⟩dz}.$$

$$⟨to0.0pt\text{hss}⟨E_{\text{k}}⟩to0.0pt\text{hss}⟩\equiv ⟨to0.0pt\text{hss}⟨\frac{1}{2}\rho {\left(\mathbf{v}-{\mathbf{v}}_{\text{K}}\right)}^2⟩to0.0pt\text{hss}⟩,$$

$$⟨to0.0pt\text{hss}⟨E_{\text{t}}⟩to0.0pt\text{hss}⟩\equiv ⟨to0.0pt\text{hss}⟨e⟩to0.0pt\text{hss}⟩,$$

$$⟨to0.0pt\text{hss}⟨E_{\text{g}}⟩to0.0pt\text{hss}⟩\equiv ⟨to0.0pt\text{hss}⟨\frac{1}{2}\rho \Phi ⟩to0.0pt\text{hss}⟩,$$

$$⟨to0.0pt\text{hss}⟨q^{+}⟩to0.0pt\text{hss}⟩\equiv ⟨to0.0pt\text{hss}⟨-p\nabla \cdot \mathbf{v}⟩to0.0pt\text{hss}⟩+⟨to0.0pt\text{hss}⟨q_{\text{shock}}⟩to0.0pt\text{hss}⟩+⟨to0.0pt\text{hss}⟨q_{\text{irr}}⟩to0.0pt\text{hss}⟩+⟨to0.0pt\text{hss}⟨q_{\text{num}}⟩to0.0pt\text{hss}⟩,$$

$$⟨to0.0pt\text{hss}⟨q^{-}⟩to0.0pt\text{hss}⟩\equiv ⟨to0.0pt\text{hss}⟨\frac{\partial F_z}{\partial z}⟩to0.0pt\text{hss}⟩+⟨to0.0pt\text{hss}⟨\frac{\partial e{v}_z}{\partial z}⟩to0.0pt\text{hss}⟩,$$

$$Q_{\text{mid}}\equiv \frac{{⟨⟨c_{\text{s}}⟩⟩}_{\rho }\Omega }{\pi G\Sigma }.$$

$${\overline{\beta }}_{\text{eff}}\equiv \frac{\overline{⟨to0.0pt\text{hss}⟨e⟩to0.0pt\text{hss}⟩}{\Omega }^{-1}}{\overline{⟨to0.0pt\text{hss}⟨-(4\pi B(T)-cE){\kappa }_{\text{P}}\rho ⟩to0.0pt\text{hss}⟩}}=3.17$$

$${\overline{\beta }}_{\text{mid}}\equiv \frac{{\overline{⟨to0.0pt\text{hss}⟨e⟩to0.0pt\text{hss}⟩}}_{\rho }{\Omega }^{-1}}{{\overline{⟨to0.0pt\text{hss}⟨-(4\pi B(T)-cE){\kappa }_{\text{P}}\rho ⟩to0.0pt\text{hss}⟩}}_{\rho }}=\mathrm{4.37.}$$

$$\frac{1}{\overline{⟨\rho ⟩}}\frac{d}{dz}\overline{⟨p⟩}+\frac{1}{\overline{⟨\rho ⟩}}\overline{⟨\frac{\partial }{\partial z}\left(\rho v_z^2\right)⟩}=-{\Omega }^{2}z-\frac{1}{\overline{⟨\rho ⟩}}\overline{⟨\rho \frac{\partial \Phi }{\partial z}⟩}.$$

$$Q=\frac{c_{\text{s}}\Omega }{\pi G\Sigma }\sim \frac{{\Omega }^{2}H}{\pi G\Sigma }\sim 2\frac{{\Omega }^{2}z}{d\Phi /dz},$$

$$\alpha \equiv \frac{W_{xy}}{\overline{⟨to0.0pt\text{hss}⟨p⟩to0.0pt\text{hss}⟩}},$$

$$\frac{3}{2}\Omega W_{xy}=\frac{3}{4\pi }{\Omega }^2\dot{M}.$$

$$0=\frac{3}{2}\Omega \overline{⟨\frac{g_x{g}_y}{4\pi G}⟩}-\overline{⟨\frac{\partial }{\partial z}\left(\rho \Phi v_z+\frac{g_z}{4\pi G}\frac{\partial \Phi }{\partial t}\right)⟩}+\overline{⟨\rho \mathbf{v}\cdot \nabla \Phi ⟩}.$$

$$\overline{⟨\frac{\partial e{v}_z}{\partial z}⟩}+\overline{⟨\frac{\partial F_z}{\partial z}⟩}=\overline{⟨q_{\text{shock}}⟩}+\overline{⟨q_{\text{num}}⟩}-\overline{⟨p\nabla \cdot \mathbf{v}⟩}+\overline{⟨q_{\text{irr}}⟩}.$$

$$\frac{N^2}{{\Omega }^2}\equiv \frac{d\ln \overline{⟨\rho ⟩}}{d\ln z}-\frac{1}{\overline{⟨{\Gamma }_1⟩}}\frac{d\ln \overline{⟨p⟩}}{d\ln z},$$

$$\frac{\partial e}{\partial t}=-(\nabla \cdot \mathbf{v})p-\frac{e}{\beta {\Omega }^{-1}}+q_{\text{irr}},$$

$${\overline{\beta }}_{\text{mid}}\sim \frac{\Sigma {\overline{⟨to0.0pt\text{hss}⟨T⟩to0.0pt\text{hss}⟩}}_{\rho }}{{\overline{⟨to0.0pt\text{hss}⟨T⟩to0.0pt\text{hss}⟩}}_{\rho }^4/\left({\overline{⟨to0.0pt\text{hss}⟨\kappa ⟩to0.0pt\text{hss}⟩}}_{\rho }\Sigma \right)}\sim {\Sigma }^{4b/3-2},$$

$$\beta =\frac{\overline{⟨to0.0pt\text{hss}⟨p⟩to0.0pt\text{hss}⟩}/({\Gamma }_1-1)}{3/2\alpha \overline{⟨to0.0pt\text{hss}⟨p⟩to0.0pt\text{hss}⟩}}=\frac{1}{\alpha },$$

$$W_{xy}\sim \alpha \Sigma {\overline{⟨to0.0pt\text{hss}⟨T⟩to0.0pt\text{hss}⟩}}_{\rho }\sim {\Sigma }^{7/3}.$$

$$\xi (\mathbf{r})\equiv \frac{\int \delta \rho (\mathbf{r}+{\mathbf{r}}^{\prime })\delta \rho ({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }}{\int \delta {\rho }^2({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }},$$

$$\frac{\lambda }{a}$$

$$=0.53Q\left(\frac{\gamma }{\mu }\right){\left(\frac{T}{15\text{ K}}\right)}^{\frac{1}{2}}{\left(\frac{a}{50\text{ AU}}\right)}^{\frac{1}{2}}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-\frac{1}{2}}.$$

Full text

††pubyear: 2016††pagerange: Gravito-turbulence in irradiated protoplanetary discs–C

1 Introduction

Angular momentum transport in protoplanetary discs controls their time evolution and thus strongly affects the planet formation process within them. In accretion discs where gas and magnetic field are well coupled, magnetic turbulence driven by the magneto-rotational instability (MRI) transports angular momentum quite efficiently (Balbus & Hawley, 1991). However, that would not be the case in protoplanetary discs, where temperatures are too low for thermal ionization to operate and thus MRI is generally suppressed by non-ideal magnetohydrodynamic effects (e.g. Turner et al., 2014).

In some cold and massive protoplanetary discs, angular momentum can be transported by shear stresses associated with the gravitational instability (GI) (see Kratter & Lodato, 2016, for a recent comprehensive review). A natural consequence of the long-range nature of gravity is formation of spiral arms as a result of GI, which globally transport angular momentum. On the other hand, Gammie (2001) showed another nonlinear outcome of GI, called gravito-turbulence, in which angular momentum transport can be described locally as in the $\alpha$ disc model (Shakura & Sunyaev, 1973). He used razor-thin, local shearing box simulations with a cooling function that has a constant cooling time $t_{\text{cool}}$ in terms of $\Omega^{-1}$, which is called the $\beta$ cooling prescription. He showed that fragmentation occurs when the cooling is rather fast as

[TABLE]

while quasi-steady gravito-turbulence of $Q\sim 1$ is sustained otherwise. Here, $Q$, defined as

[TABLE]

is called Toomre parameter, and $Q<1$ is the condition for the linear axisymmetric GI in the Keplerian disc, where $c_{\text{s}}$, $\Omega$, and $\Sigma$ are, respectively, the sound velocity, the orbital frequency, and the surface density at the radius considered (Toomre, 1964). At this gravito-turbulent phase where dissipation of the turbulence balances with the cooling, the nominal $\alpha$ parameter and $\beta$ are simply related through $\alpha=(9/4)\gamma(\gamma-1)\beta$ (Gammie, 2001).

Since then, the fragmentation criteria have been extensively studied, especially in protoplanetary discs, by many authors adopting various types of numerical methods and cooling prescriptions in both local and global simulations (e.g. Johnson & Gammie, 2003; Stamatellos & Whitworth, 2009; Cossins et al., 2010; Baehr & Klahr, 2015; Riols & Latter, 2016), mostly motivated by the idea of forming gas giants via GI (Boss, 1997, 1998; Durisen et al., 2007; Zhu et al., 2012).

However, the exact value of the critical $\beta$ for fragmentation is still an open question. The non-convergence of the fragmentation criterion may arise from numerical artifacts (Meru & Bate, 2010; Lodato & Clarke, 2011; Meru & Bate, 2012), inherent stochasticity of fragmentation (Paardekooper, 2012; Hopkins & Christiansen, 2013), the dimension (i.e. 2D vs. 3D) (Young & Clarke, 2015), or the fact that there is no physical temperature floor in the $\beta$ cooling prescription (Lin & Kratter, 2016). The irradiation can be a main heating source in cool protoplanetary discs subject to GI, and thus may affect the fragmentation criterion (Rice et al., 2011). It has also been suggested that a maximum $\alpha\sim 0.1$ might be a more general criterion than the critical $\beta$ criterion (Rice et al., 2005). Alternatively, Takahashi et al. (2016) found that fragmentation boundary determined by the Toomre parameter using global simulations with realistic thermodynamics.

Here, we present local 3D simulations of irradiated, self-gravitating discs with realistic opacities and self-consistent equation of states (EOSs). This is an extended work from Shi & Chiang (2014), who firstly performed 3D local shearing box simulations, using the $\beta$ cooling as well as simple optically-thin cooling prescriptions. We used the flux-limited diffusion approximation (FLD) in the transfer of the disc’s own radiation field while we used a ray-tracing method to calculate the irradiation heating by the central star. The aim of this paper is to explore the physics involved in the gravito-turbulence and the condition for it to be sustained for a fixed radius. We also clarify the effect of the irradiation as well as how different the realistic radiative transfer is than the simple cooling function.

This paper is organized as follows. After we describe our numerical methods in section 2, we present our numerical results in section 3. In Section 4, we mainly discuss validity of our simulations, and we conclude in section 5.

2 Methods

2.1 Basic equations and numerical schemes

The basic equations solved in our simulations are hydrodynamics equations with the Poisson equation for self-gravity and the frequency-integrated angular-moment equations of radiative transfer:

[TABLE]

where $\rho$ is the gas density, $e$ the gas internal energy, $p$ the gas pressure, $T$ the gas temperature (assumed to be the same as the dust temperature), $E$ the radiation energy density, $\mathsf{P}$ the radiation pressure tensor, $\bm{F}$ the radiation energy flux, $\bm{v}$ the velocity field vector, $B(T)=\sigma_{\text{B}}T^{4}/\pi$ the Planck function ($\sigma_{\text{B}}$, the Stefan-Boltzmann constant), and $c$ the speed of light. Under the FLD approximation, $\bm{F}$ and $\mathsf{P}$ are related to $E$ as $\bm{F}=-(c\lambda(R)/\kappa_{\text{R}}\rho)\nabla E$ and $\mathsf{P}=\mathsf{f}(R)E$. Here $\lambda(R)\equiv(2+R)/(6+3R+R^{2})$ is a flux limiter with $R\equiv|\nabla E|/(\kappa_{\text{R}}\rho E)$, and $\mathsf{f}(R)\equiv(1/2)(1-f(R))\mathsf{I}+(1/2)(3-f(R))\bm{n}\bm{n}$ is the Eddington tensor with $f(R)\equiv\lambda(R)+\lambda(R)^{2}R^{2}$ and $\bm{n}\equiv\nabla E/|\nabla E|$ (Turner & Stone, 2001).

The EOSs ($p=p(e,\rho)$ and $T=T(e,\rho)$) and the Rosseland-mean and the Planck-mean opacities ($\kappa_{\text{R}}(\rho,T)$ and $\kappa_{\text{P}}(\rho,T)$) were tabulated beforehand. The EOS tables are updated versions of those used in Tomida et al. (2013). The opacity tables are the same as those used in Hirose (2015), where dust opacities are taken from Semenov et al. (2003) while low-temperatures gas opacities are taken from Ferguson et al. (2005). The EOS and opacity tables are plotted in figure 20 in Appendix C.

We used the shearing box approximation to model a local patch of an accretion disc as a co-rotating Cartesian frame $(x,y,z)$ with a linearized Keplerian shear flow $\bm{v}_{\text{K}}\equiv-(3/2)\Omega x\hat{\bm{y}}$, where the $x$, $y$, and $z$ directions correspond to the radial, azimuthal, and vertical directions, respectively, and $\hat{\bm{y}}$ is the unit vector in the $y$ direction (Hawley et al., 1995). The inertial forces in the co-rotating frame and the vertical component of the central star’s gravity, $-2\Omega\hat{\bm{z}}\times\bm{v}+3\Omega^{2}x\hat{\bm{x}}-\Omega^{2}z\hat{\bm{z}}$, are added in the right hand side (RHS) of the equation of motion (4), where $\hat{\bm{x}}$ and $\hat{\bm{z}}$ are the unit vectors in the $x$ and $z$ direction, respectively. Shearing-periodic, periodic, and outflow boundary conditions are applied to the boundaries in the $x$, $y$, and $z$ direction, respectively. The outflow boundary condition is described in (Hirose et al., 2006).

To solve the above equations, we employed ZEUS (Stone & Norman, 1992), where we also implemented an orbital advection algorithm (Stone & Gardiner, 2010) for accurate calculations in a wide shearing box. The Poisson equation with the vacuum boundary condition in the $z$ direction was solved by Fast Fourier Transforms (Koyama & Ostriker, 2009). To test these methods, we solved one of the disc models (the one with the constant cooling time of $10\Omega^{-1}$) in Shi & Chiang (2014) and got a quantitatively consistent result.

The radiative transfer part of the basic equations is extracted as

[TABLE]

which was solved time-implicitly using a multi-grid algorithm with the Gauss-Seidel method as a smoother. During iteration in each time step, the irradiation heating rate $q_{\text{irr}}$ as well as the coefficients $\kappa_{\text{P}}\rho$, $\kappa_{\text{R}}\rho$ and $\lambda(R)$ are fixed.

The irradiation heating rate $q_{\text{irr}}$ is evaluated by solving the time-independent radiative transfer equation that ignores scattering,

[TABLE]

using a ray tracing method, with the coefficient \kappa_{\text{P{}{*}}}\rho being fixed. Here, $I$ and $s$ have the usual meanings while \kappa_{\text{P{}{}}}=\kappa_{\text{P{}_{}}}(\rho,T) is a mean opacity averaged over the Planck function of the stellar effective temperature $T_{*}$ (Hirose & Turner, 2011). See Appendix A for descritization of the equation (10) to compute $q_{\text{irr}}$.

In our simulations, the kinetic energy may dissipates either numerically in the grid-scale or physically in shock waves. In any case, the dissipated energy is captured in the form of gas internal energy, effectively resulting in additional source terms, $q_{\text{num}}$ and $q_{\text{shock}}$ respectively, in the gas energy equation (5). The shock heating rate $q_{\text{shock}}$ is computed as $-q\nabla\cdot\bm{v}$, where $q$ determines the strength of the ZEUS’s shock capturing viscosity (Stone & Norman, 1992). The method to evaluate the grid-scale dissipation rate $q_{\text{num}}$ is described in detail in Appendix in Hirose et al. (2006). Thus the sum of the kinetic and internal energies is conserved in the simulation box.

2.2 Parameters and the initial conditions

Parameters in a stratified shearing box are the orbital frequency $\Omega$ [s*-1*], which appears in the inertial force terms and the shearing periodic boundary condition, and the surface density $\Sigma$ [g cm*-2*], which is the amount of gas in the box. On the other hand, parameters for the irradiation are the energy flux $F_{\text{irr}}$ [erg cm*-2* s*-1*] and the grazing angle $\theta$ at the surfaces of the box. The surface density $\Sigma$ may vary during a simulation due to outflows through the top and bottom boundaries or the density floor described below, but the variation was typically less than $1$ % per one hundred orbits unless noted.

In this paper, a shearing box was placed at $a=50$ AU away from the central star that has the effective temperature of $T_{*}=4000$ K, the mass of $M_{*}=1M_{\odot}$, and the radius of $R_{*}=1R_{\odot}$. Then, the orbital frequency of the shearing box is determined as $\Omega=\sqrt{GM_{*}/a^{3}}=5.63\times 10^{-10}$ s*-1* while the irradiation energy flux is determined as $F_{\text{irr}}=(R_{*}/a)^{2}\sigma_{\text{B}}T_{*}^{4}=1.26\times 10^{2}$ erg cm*-2* s*-1*. The rest two parameters, the surface density $\Sigma$ and the grazing angle $\theta$, are free physical parameters in this study.

The initial disc was assumed to be isothermal and hydrostatic ignoring self-gravity. The isothermal temperature $T_{0}$ was evaluated using Chiang & Goldreich (1997)’s radiative equilibrium disc model (their equation 12a), $T_{0}=\left({\theta}/{4}\right)^{\frac{1}{4}}\left({R_{*}}/{a}\right)^{\frac{1}{2}}T_{*}$. We note that this temperature $T_{0}$ was used only in constructing the initial disc, in which the mean molecular weight $\mu=2.38$ and the adiabatic exponent $\gamma=5/3$ were also temporarily used. The initial radiation field $E_{0}$ was assumed to be in thermal equilibrium with gas, $E_{0}=aT_{0}^{4}$. The initial velocity field was the linearized Keplerian shear flow $\bm{v}_{0}=\bm{v}_{\text{K}}$, whose $x$ and $z$ components were perturbed randomly up to 0.5% of the local sound speed.

From numerical reasons, we had to introduce floors for density, internal energy, and temperature. The density floor was set to $10^{-6}$ of the initial midplane density. The internal energy floor was set basically to avoid negative values (see Appendix A3 in Hirose et al., 2006, for details). The temperature floor was set to $5$ K, which is the lower limit of the temperature range in the EOS and opacity tables. The total artificial energy injection rate due to these floors was typically less than 1% of the physical heating rate when volume-averaged (see the middle panel in figure 2).

3 Results

3.1 Diagnostics

We present results of our simulations using diagnostics based on horizontally-averaged vertical profiles, which were recorded every $0.01$ orbits.111The data for figure 18 was recorded every single orbit. The horizontally-averaged vertical profile of quantity $f$, for example, was computed as

[TABLE]

where the integrations were done over the full extent of the box in the $x$ and $y$ directions. In the diagnostics, it is time-averaged or vertically-integrated as

[TABLE]

The vertical integration is done over the full extent of the box height, and the time averaging is done between $t=20$ and $t=80$ orbits unless noted. In some diagnostics, we also use density-weighted, vertical averaging as

[TABLE]

3.2 Fiducial run

The two physical parameters, the surface density $\Sigma$ and the grazing angle $\theta$, depend on the global modeling of protoplanetary discs. Therefore, in our local shearing box simulations, they may be chosen somewhat arbitrarily. We chose $\Sigma=100$ g cm*-2* and $\theta=0.02$ for the fiducial run, where the gravito-turbulence was sustained.

The box size and the number of cells in the fiducial run were $(L_{x},L_{y},L_{z})=(24H,24H,12H)$ and $(N_{x},N_{y},N_{z})=(128,128,64)$, respectively. Here and hereafter, the scale height of the initial isothermal disc $H\equiv\sqrt{2RT_{0}/(\mu\Omega^{2})}$ is used as the unit length, where $R$ is the gas constant. It reads $H=4.75\times 10^{13}$ cm = $3.18$ AU from the physical parameters chosen in the second paragraph in section 2.2. We keep H as our unit length in all runs in this paper.

First of all, for readers to grasp what is happening in the fiducial run, we show in figure 1 typical snapshots of density, gas temperature, and shock heating rate in the quasi-steady state. It is seen that many shock waves are excited by collisions of non-axisymmetric density waves driven by GI, raising gas temperatures near the midplane.

3.2.1 Time evolution and cooling time

Figure 2 shows time evolution of vertically-integrated versions of energies (kinetic $E_{\text{k}}$, thermal $E_{\text{t}}$, and self-gravitational $E_{\text{g}}$), heating rate $q^{+}$ and cooling rate $q^{-}$ as well as the Toomre parameter near the midplane, $Q_{\text{mid}}$, which are defined as

[TABLE]

We omit the radiation energy $E$ in the thermal energy $E_{\text{t}}$, the term $-\nabla\bm{v}:\mathsf{P}$ in the heating rate $q^{+}$, and the term ${\partial Ev_{z}}/{\partial z}$ in the cooling rate $q^{-}$ since they are negligible. The fluctuations in the first few orbits are due to the deviation of the initial disc from a hydrostatic equilibrium. Then, around $t=5$ orbits, the axisymmetric mode of GI fully developed, which broke down into turbulent density waves. After this initial transient, the disc reached a statistically steady state around $t=10$ orbits.

The upper panel shows that the ratio of \mathopen{{\mathchoice{\hbox{\displaystyle\kern 8.31935pt\hbox{$\displaystyle\left<\vphantom{E_{\text{k}}}\right.$}\kern-8.31935pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{E_{\text{k}}}\right.$}}\kern 8.31935pt}}{\hbox{\textstyle\kern 8.31935pt\hbox{$\textstyle\left<\vphantom{E_{\text{k}}}\right.$}\kern-8.31935pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{E_{\text{k}}}\right.$}}\kern 8.31935pt}}{\hbox{\scriptstyle\kern 5.83989pt\hbox{$\scriptstyle\left<\vphantom{E_{\text{k}}}\right.$}\kern-5.83989pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{E_{\text{k}}}\right.$}}\kern 5.83989pt}}{\hbox{\scriptscriptstyle\kern 4.40508pt\hbox{$\scriptscriptstyle\left<\vphantom{E_{\text{k}}}\right.$}\kern-4.40508pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{E_{\text{k}}}\right.$}}\kern 4.40508pt}}}}E_{\text{k}}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 8.31935pt\hbox{$\displaystyle\left>\vphantom{E_{\text{k}}}\right.$}\kern-8.31935pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{E_{\text{k}}}\right.$}}\kern 8.31935pt}}{\hbox{\textstyle\kern 8.31935pt\hbox{$\textstyle\left>\vphantom{E_{\text{k}}}\right.$}\kern-8.31935pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{E_{\text{k}}}\right.$}}\kern 8.31935pt}}{\hbox{\scriptstyle\kern 5.83989pt\hbox{$\scriptstyle\left>\vphantom{E_{\text{k}}}\right.$}\kern-5.83989pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{E_{\text{k}}}\right.$}}\kern 5.83989pt}}{\hbox{\scriptscriptstyle\kern 4.40508pt\hbox{$\scriptscriptstyle\left>\vphantom{E_{\text{k}}}\right.$}\kern-4.40508pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{E_{\text{k}}}\right.$}}\kern 4.40508pt}}}} to \mathopen{{\mathchoice{\hbox{\displaystyle\kern 8.01796pt\hbox{$\displaystyle\left<\vphantom{E_{\text{t}}}\right.$}\kern-8.01796pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{E_{\text{t}}}\right.$}}\kern 8.01796pt}}{\hbox{\textstyle\kern 8.01796pt\hbox{$\textstyle\left<\vphantom{E_{\text{t}}}\right.$}\kern-8.01796pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{E_{\text{t}}}\right.$}}\kern 8.01796pt}}{\hbox{\scriptstyle\kern 5.62462pt\hbox{$\scriptstyle\left<\vphantom{E_{\text{t}}}\right.$}\kern-5.62462pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{E_{\text{t}}}\right.$}}\kern 5.62462pt}}{\hbox{\scriptscriptstyle\kern 4.1898pt\hbox{$\scriptscriptstyle\left<\vphantom{E_{\text{t}}}\right.$}\kern-4.1898pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{E_{\text{t}}}\right.$}}\kern 4.1898pt}}}}E_{\text{t}}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 8.01796pt\hbox{$\displaystyle\left>\vphantom{E_{\text{t}}}\right.$}\kern-8.01796pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{E_{\text{t}}}\right.$}}\kern 8.01796pt}}{\hbox{\textstyle\kern 8.01796pt\hbox{$\textstyle\left>\vphantom{E_{\text{t}}}\right.$}\kern-8.01796pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{E_{\text{t}}}\right.$}}\kern 8.01796pt}}{\hbox{\scriptstyle\kern 5.62462pt\hbox{$\scriptstyle\left>\vphantom{E_{\text{t}}}\right.$}\kern-5.62462pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{E_{\text{t}}}\right.$}}\kern 5.62462pt}}{\hbox{\scriptscriptstyle\kern 4.1898pt\hbox{$\scriptscriptstyle\left>\vphantom{E_{\text{t}}}\right.$}\kern-4.1898pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{E_{\text{t}}}\right.$}}\kern 4.1898pt}}}} is about $0.67$, which is larger by a factor of a few than that in the $\beta=10$ case in Gammie (2001). This is because the effective cooling time in our simulation is relatively shorter ($\sim 3.17\Omega^{-1}$) as will be discussed in the next paragraph. It also shows that the three energies vary almost in phase. Cross-correlations of these variations reveal that the self-gravitational energy varies first, the kinetic energy second, and the internal energy third, with mutual delays of $1.8\Omega^{-1}$ and $1.1\Omega^{-1}$, respectively. The order is reasonable as a natural consequence of the causality. Similarly, the middle panel shows that the cooling rate quickly follows the heating rate, keeping a thermal balance, but with a delay of $1.5\Omega^{-1}$. The middle panel also shows that the rate of the artificial energy injection mentioned in section 2.2 is negligible compared to the heating or cooling rates in the steady state. The lower panel shows that the Toomre parameter took its minimum of $0.7$ when the axisymmetric mode of GI fully develops, and settled down to $1.05$ on time average in the steady state.

One may be interested in what is the $\beta$ value here, the cooling time in terms of $\Omega^{-1}$ (Gammie, 2001). If we define an effective cooling time in the volume-averaged sense,

[TABLE]

while if we compute a density-weighted, volume-averaged cooling time,

[TABLE]

On the other hand, we can directly examine the cell-by-cell cooling time, defined as the internal energy $e$ divided by the energy exchange rate between gas and radiation $-{\kappa}_{\text{P}}\rho\left(4\pi B(T)-cE\right)$. Figure 3 is a two-dimensional histogram of the cell-by-cell cooling time vs. height, compiled from all cells contained in selected snapshots in the steady state. Not surprisingly, the cell-by-cell cooling time is not single-valued, but is spread over a finite range at each height. Near the midplane, the mode of the cell-by-cell cooling time is similar to the volume-averaged cooling time ($\overline{\beta}_{\text{eff}}=3.17$). On the other hand, in the upper layers that absorb the stellar irradiation, the cell-by-cell cooling time is exceptionally short ($\leq 10^{-3}$). In section 3.5, we will see that numerical results are totally different when the simple cooling function is employed in place of the FLD radiative transfer.

3.2.2 Hydrostatic balance

Figure 4 shows time-averaged vertical profiles of density $\overline{\left<\rho\right>}$, gas temperature $\overline{\left<T\right>}$, and radiation temperature $\overline{\left<(E/a)^{1/4}\right>}$, where $a\equiv 4\sigma_{\text{B}}/c$ is the radiation constant. The gas temperature near the midplane is low around $20$ K with a peak of $\sim 30$ K at the midplane. It is raised up to $60$ K in the upper layers by absorption of the irradiation of visible light by dust grains (see section 3.2.4). The radiation temperature is almost equal to gas temperature near the midplane ($|z|/H<2$), but they are apart in the upper layers where the irradiation is absorbed. The density profile has an exponential decay at $2\leq|z|/H\leq 4$ and extended tails in the irradiated hot layers ($|z|/H>4$).

Figure 5 compares time-averaged profiles of vertical accelerations that appear in the following time-averaged version of the $z$ component of the momentum equation in a steady state,

[TABLE]

We omit here the radiation force $\overline{\left<{\kappa_{\text{R}}\rho F_{z}}/{c}\right>}$ since it is negligible. Pressure gradient terms in the left-hand side (LHS) are shown in the upper panel while vertical gravity terms in the right-hand side (RHS) are shown in the lower panel. A hydrostatic balance holds well since the sum of the LHS terms and that of the RHS terms agree. Some anomalies near the boundaries are due to downflows arising from a hydrostatic imbalance.

The lower panel shows that the self-gravity dominates the external gravity at $|z|/H<4$, confirming that the disc is self-gravitating. At $|z|/H\sim 2$, the self-gravity peaks, roughly twice as large as the external gravity, which is consistent with the midplane Toomre parameter $\overline{Q}_{\text{mid}}=1.05$ through the following relation:

[TABLE]

where the Poisson equation is used as $d\Phi/dz\sim 4\pi G\rho H\sim 2\pi G\Sigma$.

On the other hand, the upper panel shows that the dynamical pressure gradient (the second term in LHS of equation 23) also competes the gravity as well as the thermal pressure gradient. The dynamical pressure here is created by upward flows driven by collisions of density waves near the midplane. Although the thermal pressure gradient is dominant near the midplane, the dynamical pressure gradient has non-negligible contribution, about a quarter of the total at $|z|/H\sim 2$. In the upper layers, the dynamical pressure gradient is even dominant.

3.2.3 Shear stresses

Figure 6 compares time-averaged vertical profiles of works done by the shear stresses. The work done by the gravitational stress $\frac{3}{2}\Omega\overline{\left<g_{x}g_{y}/4\pi G\right>}$ is the injection rate of $\left<E_{\text{g}}\right>$ into the box at each height while the work done by the Reynolds stress $\frac{3}{2}\Omega\overline{\left<\rho v_{x}\delta v_{y}\right>}$ is the injection rate of $\left<E_{\text{k}}\right>$. Here, $\bm{g}=(g_{x},g_{y},g_{z})=-\nabla\Phi$ is the self-gravitational acceleration. The Reynolds stress is slightly larger than the gravitational stress at the midplane, but the gravitational stress has a wider distribution than the Reynolds stress. When vertically-integrated, they are computed as \frac{3}{2}\Omega\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 20.43832pt\hbox{$\displaystyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\textstyle\kern 20.43832pt\hbox{$\textstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\scriptstyle\kern 13.73808pt\hbox{$\scriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-13.73808pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 13.73808pt}}{\hbox{\scriptscriptstyle\kern 10.29901pt\hbox{$\scriptscriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-10.29901pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 10.29901pt}}}}g_{x}g_{y}/4\pi G\mathclose{{\mathchoice{\hbox{\displaystyle\kern 20.43832pt\hbox{$\displaystyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\textstyle\kern 20.43832pt\hbox{$\textstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\scriptstyle\kern 13.73808pt\hbox{$\scriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-13.73808pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 13.73808pt}}{\hbox{\scriptscriptstyle\kern 10.29901pt\hbox{$\scriptscriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-10.29901pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 10.29901pt}}}}}=1.83 and \frac{3}{2}\Omega\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 14.23245pt\hbox{$\displaystyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\textstyle\kern 14.23245pt\hbox{$\textstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\scriptstyle\kern 9.99673pt\hbox{$\scriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-9.99673pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 9.99673pt}}{\hbox{\scriptscriptstyle\kern 7.62663pt\hbox{$\scriptscriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-7.62663pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 7.62663pt}}}}\rho v_{x}\delta v_{y}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 14.23245pt\hbox{$\displaystyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\textstyle\kern 14.23245pt\hbox{$\textstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\scriptstyle\kern 9.99673pt\hbox{$\scriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-9.99673pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 9.99673pt}}{\hbox{\scriptscriptstyle\kern 7.62663pt\hbox{$\scriptscriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-7.62663pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 7.62663pt}}}}}=1.51, in terms of $G\Sigma_{0}^{2}\Omega H$.

The Shakura-Sunyaev’s $\alpha$, defined as the ratio of time-averaged, vertically-integrated shear stress to time-averaged, vertically-integrated thermal pressure, is commonly used to evaluate the shear stress in accretion discs (Shakura & Sunyaev, 1973). In this paper, we define $\alpha$ as

[TABLE]

which is different from the conventional $\alpha\equiv\left|d\ln\Omega/d\ln r\right|^{-1}W_{xy}/(\Sigma c_{\text{s}}^{2})$ in the literature by a factor of $2/3\gamma$.

Here, it is computed as $0.138$ and $0.105$ for the self-gravitational and Reynolds stresses, respectively.

Another way to evaluate the shear stress is to express it in terms of the mass accretion rate. Assuming a time-steady accretion, the vertically-integrated shear stress $W_{xy}$ and the mass accretion rate $\dot{M}$ are connected through

[TABLE]

Using this relation, the mass accretion rate in the fiducial run is computed as $5.53\times 10^{-6}$ $M_{\odot}$ yr*-1*.

Figure 6 also tells about energetics of the self-gravitational energy. The time and horizontally-averaged version of the self-gravitational energy equation in a steady state is written as follows (see Appendix B):

[TABLE]

The third term is the work done by the self-gravitational force, which is the conversion rate of $\left<E_{\text{g}}\right>$ to $\left<E_{\text{k}}\right>$ at each height. As seen in the figure, the first term (black solid) and the third term (black dotted) do not cancel, meaning that the self-gravitational energy injected via the gravitational stress is redistributed via the second term (black dashed) before converted to the kinetic energy.

3.2.4 Thermal balance

The time and horizontally-averaged thermal energy equation in a steady state can be written as

[TABLE]

Here, heating (cooling) terms are, respectively, gathered in RHS (LHS), and are shown in the upper (lower) panel in figure 7. The profiles of the total heating rate and the total cooling rate match, which means that a thermal balance well holds.

First, we examine the heating rates. The gravito-turbulence here dissipates mainly through shock heating $\overline{\left<q_{\text{shock}}\right>}$. However, the compressional heating (or radiation damping), $-\overline{\left<p\nabla\cdot\bm{v}\right>}$, also plays a major role; it is actually comparable to the shock heating near the midplane. We note that these two major dissipation processes are not quite dependent on the grid resolution unlike the grid-scale dissipation $\overline{\left<{q}_{\text{num}}\right>}$, which is minor here as shown in the figure. The irradiation heating $\overline{\left<q_{\text{irr}}\right>}$ only occurs above the photosphere of the visible light at $|z|/H\sim 3$. When vertically integrated, each heating rate is computed as \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 11.2153pt\hbox{$\displaystyle\left<\vphantom{q_{\text{shock}}}\right.$}\kern-11.2153pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{q_{\text{shock}}}\right.$}}\kern 11.2153pt}}{\hbox{\textstyle\kern 11.2153pt\hbox{$\textstyle\left<\vphantom{q_{\text{shock}}}\right.$}\kern-11.2153pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{q_{\text{shock}}}\right.$}}\kern 11.2153pt}}{\hbox{\scriptstyle\kern 7.92578pt\hbox{$\scriptstyle\left<\vphantom{q_{\text{shock}}}\right.$}\kern-7.92578pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{q_{\text{shock}}}\right.$}}\kern 7.92578pt}}{\hbox{\scriptscriptstyle\kern 6.73396pt\hbox{$\scriptscriptstyle\left<\vphantom{q_{\text{shock}}}\right.$}\kern-6.73396pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{q_{\text{shock}}}\right.$}}\kern 6.73396pt}}}}q_{\text{shock}}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 11.2153pt\hbox{$\displaystyle\left>\vphantom{q_{\text{shock}}}\right.$}\kern-11.2153pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{q_{\text{shock}}}\right.$}}\kern 11.2153pt}}{\hbox{\textstyle\kern 11.2153pt\hbox{$\textstyle\left>\vphantom{q_{\text{shock}}}\right.$}\kern-11.2153pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{q_{\text{shock}}}\right.$}}\kern 11.2153pt}}{\hbox{\scriptstyle\kern 7.92578pt\hbox{$\scriptstyle\left>\vphantom{q_{\text{shock}}}\right.$}\kern-7.92578pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{q_{\text{shock}}}\right.$}}\kern 7.92578pt}}{\hbox{\scriptscriptstyle\kern 6.73396pt\hbox{$\scriptscriptstyle\left>\vphantom{q_{\text{shock}}}\right.$}\kern-6.73396pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{q_{\text{shock}}}\right.$}}\kern 6.73396pt}}}}}=0.888, -\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 13.0085pt\hbox{$\displaystyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-13.0085pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 13.0085pt}}{\hbox{\textstyle\kern 13.0085pt\hbox{$\textstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-13.0085pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 13.0085pt}}{\hbox{\scriptstyle\kern 9.10596pt\hbox{$\scriptstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-9.10596pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 9.10596pt}}{\hbox{\scriptscriptstyle\kern 6.50424pt\hbox{$\scriptscriptstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-6.50424pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 6.50424pt}}}}p\nabla\cdot\bm{v}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 13.0085pt\hbox{$\displaystyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-13.0085pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 13.0085pt}}{\hbox{\textstyle\kern 13.0085pt\hbox{$\textstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-13.0085pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 13.0085pt}}{\hbox{\scriptstyle\kern 9.10596pt\hbox{$\scriptstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-9.10596pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 9.10596pt}}{\hbox{\scriptscriptstyle\kern 6.50424pt\hbox{$\scriptscriptstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}\kern-6.50424pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{p\nabla\cdot\bm{v}}\right.$}}\kern 6.50424pt}}}}}=0.510, \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 10.11826pt\hbox{$\displaystyle\left<\vphantom{{q}_{\text{num}}}\right.$}\kern-10.11826pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{{q}_{\text{num}}}\right.$}}\kern 10.11826pt}}{\hbox{\textstyle\kern 10.11826pt\hbox{$\textstyle\left<\vphantom{{q}_{\text{num}}}\right.$}\kern-10.11826pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{{q}_{\text{num}}}\right.$}}\kern 10.11826pt}}{\hbox{\scriptstyle\kern 7.1422pt\hbox{$\scriptstyle\left<\vphantom{{q}_{\text{num}}}\right.$}\kern-7.1422pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{{q}_{\text{num}}}\right.$}}\kern 7.1422pt}}{\hbox{\scriptscriptstyle\kern 5.95036pt\hbox{$\scriptscriptstyle\left<\vphantom{{q}_{\text{num}}}\right.$}\kern-5.95036pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{{q}_{\text{num}}}\right.$}}\kern 5.95036pt}}}}{q}_{\text{num}}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 10.11826pt\hbox{$\displaystyle\left>\vphantom{{q}_{\text{num}}}\right.$}\kern-10.11826pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{{q}_{\text{num}}}\right.$}}\kern 10.11826pt}}{\hbox{\textstyle\kern 10.11826pt\hbox{$\textstyle\left>\vphantom{{q}_{\text{num}}}\right.$}\kern-10.11826pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{{q}_{\text{num}}}\right.$}}\kern 10.11826pt}}{\hbox{\scriptstyle\kern 7.1422pt\hbox{$\scriptstyle\left>\vphantom{{q}_{\text{num}}}\right.$}\kern-7.1422pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{{q}_{\text{num}}}\right.$}}\kern 7.1422pt}}{\hbox{\scriptscriptstyle\kern 5.95036pt\hbox{$\scriptscriptstyle\left>\vphantom{{q}_{\text{num}}}\right.$}\kern-5.95036pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{{q}_{\text{num}}}\right.$}}\kern 5.95036pt}}}}}=0.137, and \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 8.2617pt\hbox{$\displaystyle\left<\vphantom{q_{\text{irr}}}\right.$}\kern-8.2617pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{q_{\text{irr}}}\right.$}}\kern 8.2617pt}}{\hbox{\textstyle\kern 8.2617pt\hbox{$\textstyle\left<\vphantom{q_{\text{irr}}}\right.$}\kern-8.2617pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{q_{\text{irr}}}\right.$}}\kern 8.2617pt}}{\hbox{\scriptstyle\kern 5.81607pt\hbox{$\scriptstyle\left<\vphantom{q_{\text{irr}}}\right.$}\kern-5.81607pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{q_{\text{irr}}}\right.$}}\kern 5.81607pt}}{\hbox{\scriptscriptstyle\kern 4.62425pt\hbox{$\scriptscriptstyle\left<\vphantom{q_{\text{irr}}}\right.$}\kern-4.62425pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{q_{\text{irr}}}\right.$}}\kern 4.62425pt}}}}q_{\text{irr}}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 8.2617pt\hbox{$\displaystyle\left>\vphantom{q_{\text{irr}}}\right.$}\kern-8.2617pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{q_{\text{irr}}}\right.$}}\kern 8.2617pt}}{\hbox{\textstyle\kern 8.2617pt\hbox{$\textstyle\left>\vphantom{q_{\text{irr}}}\right.$}\kern-8.2617pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{q_{\text{irr}}}\right.$}}\kern 8.2617pt}}{\hbox{\scriptstyle\kern 5.81607pt\hbox{$\scriptstyle\left>\vphantom{q_{\text{irr}}}\right.$}\kern-5.81607pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{q_{\text{irr}}}\right.$}}\kern 5.81607pt}}{\hbox{\scriptscriptstyle\kern 4.62425pt\hbox{$\scriptscriptstyle\left>\vphantom{q_{\text{irr}}}\right.$}\kern-4.62425pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{q_{\text{irr}}}\right.$}}\kern 4.62425pt}}}}}=0.281, in terms of $G\Sigma_{0}^{2}\Omega H$.

Here we describe how the temperature structure shown in figure 4 is formed. The heated upper layers emit the blackbody radiation (sometimes called reprocessed irradiation), and a part of it heats the midplane region. The base temperature in the region between the heated upper layers is determined by this reprocessed irradiation. On the other hand, the peak temperature of $30$ K at the midplane is determined by dissipation of the turbulence (described in the previous paragraph). Therefore, if there is no dissipation of the turbulence, the temperature profile near the midplane will be flat at the base temperature. This is actually seen when the surface density is small or the grazing angle is large (see sections 3.3 and 3.4; see also Hirose & Turner, 2011).

Next we examine the cooling rates. The radiative cooling $\overline{\left<{\partial F_{z}}/{\partial z}\right>}$ dominates the advective cooling $\overline{\left<{\partial ev_{z}}/{\partial z}\right>}$ at all heights. Actually, the advection picks up considerable amount of the dissipated energy near the midplane ($\overline{\left<{\partial ev_{z}}/{\partial z}\right>}$ is positive), but dumps all of it around $|z|/H\sim 1$ ($\overline{\left<{\partial ev_{z}}/{\partial z}\right>}$ is negative), which then the radiative diffusion takes over. When vertically-integrated, they are computed as \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 18.402pt\hbox{$\displaystyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-18.402pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 18.402pt}}{\hbox{\textstyle\kern 18.402pt\hbox{$\textstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-18.402pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 18.402pt}}{\hbox{\scriptstyle\kern 12.29442pt\hbox{$\scriptstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-12.29442pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 12.29442pt}}{\hbox{\scriptscriptstyle\kern 9.00713pt\hbox{$\scriptscriptstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-9.00713pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 9.00713pt}}}}{\partial F_{z}}/{\partial z}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 18.402pt\hbox{$\displaystyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-18.402pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 18.402pt}}{\hbox{\textstyle\kern 18.402pt\hbox{$\textstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-18.402pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 18.402pt}}{\hbox{\scriptstyle\kern 12.29442pt\hbox{$\scriptstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-12.29442pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 12.29442pt}}{\hbox{\scriptscriptstyle\kern 9.00713pt\hbox{$\scriptscriptstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}\kern-9.00713pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{{\partial F_{z}}/{\partial z}}\right.$}}\kern 9.00713pt}}}}}=1.82 and \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 19.19357pt\hbox{$\displaystyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-19.19357pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 19.19357pt}}{\hbox{\textstyle\kern 19.19357pt\hbox{$\textstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-19.19357pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 19.19357pt}}{\hbox{\scriptstyle\kern 12.84853pt\hbox{$\scriptstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-12.84853pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 12.84853pt}}{\hbox{\scriptscriptstyle\kern 9.40291pt\hbox{$\scriptscriptstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-9.40291pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 9.40291pt}}}}{\partial ev_{z}}/{\partial z}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 19.19357pt\hbox{$\displaystyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-19.19357pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 19.19357pt}}{\hbox{\textstyle\kern 19.19357pt\hbox{$\textstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-19.19357pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 19.19357pt}}{\hbox{\scriptstyle\kern 12.84853pt\hbox{$\scriptstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-12.84853pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 12.84853pt}}{\hbox{\scriptscriptstyle\kern 9.40291pt\hbox{$\scriptscriptstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}\kern-9.40291pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{{\partial ev_{z}}/{\partial z}}\right.$}}\kern 9.40291pt}}}}}=7\times 10^{-6}, in terms of $G\Sigma_{0}^{2}\Omega H$. The fact that the vertically-integrated advective cooling is negligible means that the dissipated energy is eventually ejected from the box via radiation only. In the lower panel, the profile of the hydrodynamical Brunt–Väisälä frequency $N$ squared (divided by $\Omega$ squared),

[TABLE]

is also shown, where $\Gamma_{1}\equiv d\ln p/d\ln\rho$ is the generalized adiabatic exponent.222$\Gamma_{1}$ is also precomputed as a function of $\rho$ and $e/\rho$ as one of the EOS tables. Since it is consistently positive near the midplane (note that the axis for $N^{2}/\Omega^{2}$ is on the right), the advection is not associated with thermal convection, but is driven by the upward flows created by the collisions of the turbulent density waves. We remind readers that this upward flows also created the dynamical pressure in the hydrostatic balance (section 3.2.2).

3.2.5 Turbulent velocity and sound velocity

Knowing the detailed turbulence properties are crucial to study the dust dynamics in self-gravitating protoplanetary discs (Booth & Clarke, 2016; Shi & Stone, 2016). Figure 8 shows time-averaged vertical profiles of the turbulent velocity defined as $\delta\bm{v}\equiv(\bm{v}-\bm{v}_{\text{K}})-\left<\bm{v}-\bm{v}_{\text{K}}\right>$ as well as the sound velocity $c_{\text{s}}\equiv\sqrt{\Gamma_{1}p/\rho}$. The figure reveals that the gravito-turbulence here is supersonic at all heights including the midplane.333The profiles of the turbulent velocities are not symmetric due to downflows near the boundaries that were consistent for some periods in this particular simulation. This is consistent with Shi & Chiang (2014), in which the turbulent flows are supersonic when the cooling time is rather short $\beta\leq 4$ (c.f. $\overline{\beta}_{\text{eff}}=3.17$ in our simulation). The turbulent velocity varies with disc height by a factor of a few, and the $x$-component of the turbulent velocity contributes most near the midplane, which are also consistent with Shi & Chiang (2014).

As discussed in Gammie (2001), the ratio of the sound velocity to the Keplerian velocity needs to be small enough for the local model to be applicable. The Keplerian velocity in our study is computed as $v_{\text{K}}={a\Omega}=({a}/{H}){\Omega H}=15.7{\Omega H}$. Therefore the ratio of the sound velocity, which reads $\sim 1.5\Omega H$ in figure 8, to the Keplerian velocity is about $0.096$. This might not be small enough since Gammie (2001) derived that the ratio needs to be much less than 0.12 based on Fourier analysis of the surface density in his simulations. We will discuss this issue from different point of views in section 4.3.

3.3 Dependence on the surface density $\Sigma$

In this section, we examine dependence on the surface density $\Sigma$. Here, we change $\Sigma$ from $30$ to $300$ (specifically 30, 60, 80, 100, 150, 200, 250, and 300) g cm*-2*, with all other parameters being fixed.

When $\Sigma\leq 60$ g cm*-2*, the flow is found to be laminar. In the range of $80\leq\Sigma\leq 250$ g cm*-2*, the gravito-turbulence is sustained. When $\Sigma=300$ g cm*-2*, a gravitationally-bounded clump was formed during the initial transient, which increases its mass and eventually underwent a rapid collapse around $t=2.4$ orbits. The collapse, which is actually the same physical process as the first core collapse in the star formation, couldn’t be resolved with our fixed numerical grid, the simulation was stopped there (see details for Appendix C). In figure 9, we compare snapshots of density and gas temperature for $\Sigma=60$, $100$, $200$, and $300$ g cm*-3*. The snapshot for $\Sigma=300$ g cm*-3* was taken just before the simulation was stopped.

Figure 10 shows how various quantities depend on $\Sigma$. The first panel shows the dependence of the time-averaged, vertically-integrated stress W_{xy}\equiv\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 20.43832pt\hbox{$\displaystyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\textstyle\kern 20.43832pt\hbox{$\textstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\scriptstyle\kern 13.73808pt\hbox{$\scriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-13.73808pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 13.73808pt}}{\hbox{\scriptscriptstyle\kern 10.29901pt\hbox{$\scriptscriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-10.29901pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 10.29901pt}}}}g_{x}g_{y}/4\pi G\mathclose{{\mathchoice{\hbox{\displaystyle\kern 20.43832pt\hbox{$\displaystyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\textstyle\kern 20.43832pt\hbox{$\textstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-20.43832pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 20.43832pt}}{\hbox{\scriptstyle\kern 13.73808pt\hbox{$\scriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-13.73808pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 13.73808pt}}{\hbox{\scriptscriptstyle\kern 10.29901pt\hbox{$\scriptscriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}\kern-10.29901pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{g_{x}g_{y}/4\pi G}\right.$}}\kern 10.29901pt}}}}}+\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 14.23245pt\hbox{$\displaystyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\textstyle\kern 14.23245pt\hbox{$\textstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\scriptstyle\kern 9.99673pt\hbox{$\scriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-9.99673pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 9.99673pt}}{\hbox{\scriptscriptstyle\kern 7.62663pt\hbox{$\scriptscriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-7.62663pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 7.62663pt}}}}\rho v_{x}\delta v_{y}\mathclose{{\mathchoice{\hbox{\displaystyle\kern 14.23245pt\hbox{$\displaystyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\textstyle\kern 14.23245pt\hbox{$\textstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-14.23245pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 14.23245pt}}{\hbox{\scriptstyle\kern 9.99673pt\hbox{$\scriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-9.99673pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 9.99673pt}}{\hbox{\scriptscriptstyle\kern 7.62663pt\hbox{$\scriptscriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}\kern-7.62663pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{\rho v_{x}\delta v_{y}}\right.$}}\kern 7.62663pt}}}}} in terms of $\dot{M}$ via equation (26). When the flow is laminar in $\Sigma\leq 60$ g cm*-2*, the values of $\dot{M}$ is below the displayed range. When the gravito-turbulence is sustained, $\dot{M}$ is strongly correlated with $\Sigma$; specifically, $\dot{M}(\Sigma)\sim\Sigma^{7/3}$ for $100\leq\Sigma<300$ g cm*-2*, where $\dot{M}$ changes from $10^{-5}$ to $10^{-4}$ $M_{\odot}$ yr*-1*. In that range of $\Sigma$, the fraction of the Reynolds stress in the total stress gradually increases with $\Sigma$, from $0.45$ to $0.56$.

The Shakura & Sunyaev’s $\alpha$ is almost constant at $\sim 0.25$ in the gravito-turbulence regime as shown in the second panel in figure 10. This means that the equilibrium thermal pressure also rises as the surface density increases, in a similar way to the stress $W_{xy}$ shown in the first panel. We note that $\alpha$ here is different from the conventional $\alpha$ by a factor of $2/3\gamma$ (see equation 25); the value of $\sim 0.25$ here therefore corresponds to $\sim 0.1$ in terms of the conventional $\alpha$ with $\gamma=5/3$. Although this is rather close to the maximum $\alpha$ sustainable in the gravito-turbulence found by Rice et al. (2005), we note that the value of $\alpha$ may change with the radius (see discussion in section 5).

The third panel in figure 10 shows how the time-averaged cooling time near the midplane $\overline{\beta}_{\text{mid}}$, defined in equation (22), changes with $\Sigma$. In the laminar flow range of $\Sigma$ ($\leq 60$ g cm*-2*), since there is little dissipation (and thus little cooling rate) while the thermal energy is kept finite due to the reprocessed irradiation, $\overline{\beta}_{\text{mid}}$ become very large and are not displayed. When the gravito-turbulence is sustained, the lowest $\Sigma$ case exhibits a larger $\overline{\beta}_{\text{mid}}$ for a similar reason while in other cases $\overline{\beta}_{\text{mid}}$ tends to be constant at $\sim 4$.

The bottom panel in figure 10 shows the time-averaged Toomre parameter near the midplane, $\overline{Q}_{\text{mid}}$, as a function of $\Sigma$. The initial ${Q}_{\text{mid}}$ (small open circles) is proportional to $\Sigma^{-1}$ since the initial sound velocity is the same in all cases. For $\Sigma\geq 100$ g cm*-2*, where the gravito-turbulence is sustained, $\overline{Q}_{\text{mid}}$ is larger than the initial value because the sound velocity is raised by dissipation of the turbulence. More importantly, $\overline{Q}_{\text{mid}}$ decreases (as $\sim\Sigma^{-1/3}$) down to $\sim 0.7$ at $\Sigma=250$ g cm*-2*, beyond which fragmentation occurred. This indicates that the critical value of $\Sigma$ corresponds to the minimum value of the Toomre parameter that can be realized in the gravito-turbulence. On the other hand, the critical value of $\Sigma$ seems not related to a specific value of the cooling time $\overline{\beta}_{\text{mid}}$ since it stays almost constant in the gravito-turbulence regime.

Figure 11 shows dependence on $\Sigma$ of the vertical profile of gas temperature.444We did not plot the radiation temperature for clarity. As in the fiducial run (figure 4), the radiation temperature always follows the gas temperature near the midplane, and continues to decrease monotonically in the upper layers. When the flow is laminar ($\Sigma\leq 60$ g cm*-2*), the temperature of the disc interior is flat at the base temperature ($T\sim 10$ K) that is determined by the reprocessed irradiation from the hot upper layers ($T\sim 60$ K) (Chiang & Goldreich, 1997).555The location of the boundaries between the upper layers and the interior are lower in the case of larger $\Sigma$ since downflows in the upper layers due to hydrostatic imbalance are stronger and thus the disc is more compressed. On the other hand, when the gravito-turbulence is sustained for $\Sigma\geq 80$ g cm*-2*, the interior is also heated by dissipation of the turbulence and there appears a peak at the midplane. As $\Sigma$ is increased, the midplane temperature increases, and thus the wave length of the axisymmetric mode of GI becomes longer as seen in the snapshots (the second and third panels) in figure 9. The midplane temperature takes a maximum value of $\sim 100$ K at $\Sigma=250$ g cm*-2*. Therefore, at least at this radius of 50 AU, heating by the gravito-turbulence will not activate MRI.

3.4 Dependence on the grazing angle $\theta$

In this section, we will see dependence on the grazing angle $\theta$, ranging from $0.005$ to $0.64$, for the surface density $\Sigma=100$ g cm*-2*, the same as in the fiducial run. Note that the amount of irradiation energy injected into the simulation box per unit time is also changed as $\sim\sin\theta$ (see Appendix A). The initial condition of these simulations was taken from a snapshot in the steady state of the fiducial run. (Therefore, the unit length $H$ is unchanged.)

The top panel in figure 12 shows dependence of the total stress $W_{xy}$ and $\alpha$. As $\theta$ is increased, both stay almost constant until $\theta\sim 0.1$, beyond which they decrease gradually, and then drop suddenly at $\theta=0.32$. At $\theta=0.64$, the gravito-turbulence is very weak and the total stress $W_{xy}$ is about $1/20$ the constant value at the lower $\theta$s. Figure 13 shows snapshots of density and gas temperature for $\theta=0.64$ and $0.24$ as well as for the fiducial run ($\theta=0.02$).

The middle and bottom panels in figure 12 shows, respectively, how vertical profile of gas temperature and irradiation heating changes with $\theta$. When $\theta\leq 0.24$, the temperature profiles are similar, having a small peak at the midplane, which means that the main heating source there is the dissipation of the turbulence. Temperatures at $|z|/H\sim 4$ rise as $\theta$ increases since the incoming irradiation energy increases as $\sim\sin\theta$. When $\theta\geq 0.32$, the temperatures near the midplane are greatly affected by the reprocessed irradiation. At $\theta=0.32$, the peak temperature at the midplane is reduced because the gravito-turbulence is weakened by the heat of the reprocessed irradiation. At $\theta=0.64$, the gravito-turbulence is almost shut off and temperatures near the midplane are mainly determined by the reprocessed irradiation.

3.5 Cases with the simple cooling function

So far we have shown results using our radiative transfer solver with the FLD approximation. For comparison, we show results of simulations using the simple cooling function commonly used in gravito-turbulence simulations (e.g. Gammie, 2001). In those simulations, we simply replaced the energy equation (5) with the following equation and solve it using the predictor-corrector method:

[TABLE]

where $\beta$ is a constant cooling time. At the same time, we dropped off the radiation energy equation (6) and the radiation force in the momentum equation (4). We performed four cases of $\beta=30$, $10$, $3$, and $1$, where a snapshot of the fiducial run was used as the initial condition. After an initial transient that lasted for several orbits, they reached a steady state except for the $\beta=1$ case. In the $\beta=1$ case, about 10% of the total mass was lost via the vertical boundaries during the first 10 orbits, and thus we stopped the calculation there.666The time averaging analysis was done for this first 10 orbits.

Snapshots of density and gas temperature are compared in figure 14. Comparing with the FLD case (shown in the bottom for reference), the flow structures look quite different even in the case of $\beta=3$, which is most close to the FLD case in terms of $\beta$ (see equation 21). Especially, fine structures are seen in the cases with the simple cell-by-cell cooling function while structures look more diffusive in the FLD case, presumably due to radiative diffusion. These indicates that the simple cooling function does not approximate the realistic radiative transfer.

The top panel of figure 15 shows profiles of time-averaged total stress for the four cases as well as the FLD case. The shape of the profiles is similar to that of the FLD case except for the $\beta=1$ case. The profile of the $\beta=3$ case is quantitatively similar to the FLD case, in which $\overline{\beta}_{\text{eff}}=3.17$. In the case of $\beta=1$, the profile is irregular due to fragmentation as seen in the snapshot (figure 14), indicating that the fragmentation criterion is similar to the Gammie (2001)’s condition (equation 1) when the irradiation is included.

The middle panel of figure 15 shows profiles of time-averaged gas temperature. Unlike the stress in the above, the profiles of temperature of the four cases are completely different from the profile of the FLD case. As $\beta$ decreases, temperatures in the upper layers decreases because the cooling time is shorter, but the midplane temperature increases according to the increase of stresses. More importantly, temperatures except the midplane are consistently much higher than the FLD case, indicating that the simple cooling function with a constant $\beta$ value does not represent the realistic cooling regardless of the $\beta$ value. We remind readers that the cell-by-cell cooling time in the FLD run is actually quite short ($\beta\ll 1$) except the midplane (figure 3). In the cases of $\beta\geq 3$, it is seen that the gas temperatures in the upper layers are saturated just above $10^{3}$ K. This is because increase in the internal energy is used to dissociate H2 molecules, rather than to increase the translational energy of the molecules, at that temperature (see $\Gamma_{1}(\rho,T)$ for the upper panel in figure 20).

The bottom panel of figure 15 shows profiles of time-averaged density. In the simple cooling function cases, corresponding to the high temperatures at higher altitudes shown in the middle panel, the density scale height there is generally larger than that in the FLD case. On the other hand, the density profiles near the midplane are similar (except the case of $\beta=1$). This is the reason why the stress profiles are similar as shown in the top panel because the profiles of gravitational potential are mostly determined by the density profiles near the midplane.

4 Discussion

4.1 Numerical convergence of vertical profiles of stress work and gas temperature

In figure 16, the upper panel shows dependence of the vertical profile of the stress work on the horizontal box size $L_{x}(=L_{y})$, the vertical box size $L_{z}$, and the resolution. The solid curves correspond to the cases that have the same box size as the fiducial ($L_{x}/H=24$), but are different in the vertical box size (blue: $L_{z}/H=18$, red: $L_{z}/H=9$) or in the resolution (grey: two times finer than the fiducial). Since they show similar profiles, we see that the vertical box size or the resolution does not affect much the result. On the other hand, the black curves compare the cases that have different horizontal box sizes. The smallest box run (dash-dotted: $L_{x}/H=12$) shows consistently smaller stress work than others. The other cases ($L_{x}/H\geq 18$) show similar profiles near the midplane, but the stress work in the upper layers ($z/H\geq 2$) increases as $L_{x}/H$ is increased.

The different amounts of the stress work in the upper layers do not affect temperatures there since they are mainly determined by the irradiation heating. This can be seen by comparing the black curves in the lower panel, where the vertical profiles of gas temperature are compared using the same notation. On the other hand, comparing the solid curves in the lower panel, it is seen that the temperature profile in the upper layers is affected by the vertical box size. This is because we assume that the irradiation always begins to be absorbed at the vertical boundaries, and thus the temperatures there are virtually fixed at the same value ($\sim 60$ K). In summary, at least for the gravito-turbulence in the main body of the disc, which we are interested in, we may conclude that our results are firm.

4.2 Dependence on the surface density

In section 3.3, we explored dependence of some key quantities on the surface density and found that i) $W_{xy}(\Sigma)\sim\Sigma^{7/3}$, ii) $\alpha$ and $\overline{\beta}_{\text{mid}}$ stays almost constant, and iii) $\overline{Q}_{\text{mid}}(\Sigma)\sim\Sigma^{-1/3}$, in the regime of gravito-turbulence. These properties seem robust since we confirm that they hold also for a different grazing angle ($\theta=0.32$) and for a different box size ($L_{x}/H=12$), as shown in figure 17. (The turbulence is weaker when $\Sigma\leq 100$ g cm*-2* for $\theta=0.32$ since the irradiation affects the midplane temperature more when $\Sigma$ is smaller.) Although the smaller box size ($L_{x}/H=12$) case shows smaller values in $W_{xy}$ (and thus smaller $\alpha$s and larger $\overline{\beta}_{\text{mid}}$s), which is a box size effect as discussed in the previous section, the slope is the same as others, and $\overline{Q}_{\text{mid}}(\Sigma)$ is even quantitatively similar to others.

Here, we examine what determines those scalings on $\Sigma$ found in our simulations. First, we note the weak scaling of the Toomre parameter on $\Sigma$. This is contrasted with the assumption of $Q\sim 1$ in the literature that discuss scalings in the gravito-turbulence (c.f. Levin, 2007; Clarke, 2009; Rafikov, 2009; Paardekooper, 2012). Given the weak scaling of the Toomre parameter on $\Sigma$, we have \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}\sim\Sigma^{4/3} from equation (2). Then, the cooling time $\overline{\beta}_{\text{mid}}$, approximated as the thermal energy content divided by the cooling rate by radiative diffusion,777In our simulations, the primary cooling mechanism is still radiative diffusion although advection also contributes to cooling as discussed in section 3.2.4. can be written as a function of $\Sigma$ as

[TABLE]

where we assumed \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}\kappa\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}}_{\rho}\sim\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}^{b}. We see that the rather constant $\overline{\beta}_{\text{mid}}$ ($\sim\Sigma^{0}$) found in our simulations corresponds to $b=3/2$. This is roughly confirmed in figure 18, where the solutions of \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}\gtrsim 30 K correspond to those following the scalings shown in figure 17. Thus, the constant $\overline{\beta}_{\text{mid}}$ can be explained as a consequence of $\overline{Q}_{\text{mid}}(\Sigma)\sim\Sigma^{-1/3}$ and \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}\kappa\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}}_{\rho}\sim\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}^{3/2}. We note that the same constant $\beta\sim\Sigma^{0}$ is also derived with the assumption of $Q\sim 1$ and $\kappa\sim T^{2}$ at given radius (Paardekooper, 2012).

The scaling \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}\kappa\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.32285pt\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\textstyle\kern 6.32285pt\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}\kern-6.32285pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{\kappa}\right.$}}\kern 6.32285pt}}{\hbox{\scriptstyle\kern 4.426pt\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}\kern-4.426pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{\kappa}\right.$}}\kern 4.426pt}}{\hbox{\scriptscriptstyle\kern 3.1614pt\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}\kern-3.1614pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{\kappa}\right.$}}\kern 3.1614pt}}}}}_{\rho}\sim\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}^{3/2} for \overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left<\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left<\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{T}\right.$}}\kern 3.44643pt}}}}T\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.89288pt\hbox{$\displaystyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\textstyle\kern 6.89288pt\hbox{$\textstyle\left>\vphantom{T}\right.$}\kern-6.89288pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{T}\right.$}}\kern 6.89288pt}}{\hbox{\scriptstyle\kern 4.82501pt\hbox{$\scriptstyle\left>\vphantom{T}\right.$}\kern-4.82501pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{T}\right.$}}\kern 4.82501pt}}{\hbox{\scriptscriptstyle\kern 3.44643pt\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}\kern-3.44643pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{T}\right.$}}\kern 3.44643pt}}}}}_{\rho}\leq 100 K found in our simulations is a result of mixing $\kappa_{\text{R}}(\rho,T)\sim T^{2}$ for $T\leq 100$ K and $\kappa_{\text{R}}(\rho,T)\sim T^{0}$ for $T\geq 100$ K (e.g. Johnson & Gammie, 2003, see also the inset in the lower panel in figure 20). This is because the temperature $T$ can be locally larger than $100$ K even when its horizontal average $\left<T\right>$ is less than $100$ K, as seen in the second and third panels in the right column in figure 9.

As discussed in Gammie (2001), a thermal balance condition that the cooling rate equals to the stress work 3/2\Omega W_{xy}=3/2\Omega\alpha\overline{\mathopen{{\mathchoice{\hbox{\displaystyle\kern 6.03984pt\hbox{$\displaystyle\left<\vphantom{p}\right.$}\kern-6.03984pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left<\vphantom{p}\right.$}}\kern 6.03984pt}}{\hbox{\textstyle\kern 6.03984pt\hbox{$\textstyle\left<\vphantom{p}\right.$}\kern-6.03984pt\hbox to0.0pt{\hss\hbox{$\textstyle\left<\vphantom{p}\right.$}}\kern 6.03984pt}}{\hbox{\scriptstyle\kern 4.22789pt\hbox{$\scriptstyle\left<\vphantom{p}\right.$}\kern-4.22789pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left<\vphantom{p}\right.$}}\kern 4.22789pt}}{\hbox{\scriptscriptstyle\kern 3.01991pt\hbox{$\scriptscriptstyle\left<\vphantom{p}\right.$}\kern-3.01991pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left<\vphantom{p}\right.$}}\kern 3.01991pt}}}}p\mathclose{{\mathchoice{\hbox{\displaystyle\kern 6.03984pt\hbox{$\displaystyle\left>\vphantom{p}\right.$}\kern-6.03984pt\hbox to0.0pt{\hss\hbox{$\displaystyle\left>\vphantom{p}\right.$}}\kern 6.03984pt}}{\hbox{\textstyle\kern 6.03984pt\hbox{$\textstyle\left>\vphantom{p}\right.$}\kern-6.03984pt\hbox to0.0pt{\hss\hbox{$\textstyle\left>\vphantom{p}\right.$}}\kern 6.03984pt}}{\hbox{\scriptstyle\kern 4.22789pt\hbox{$\scriptstyle\left>\vphantom{p}\right.$}\kern-4.22789pt\hbox to0.0pt{\hss\hbox{$\scriptstyle\left>\vphantom{p}\right.$}}\kern 4.22789pt}}{\hbox{\scriptscriptstyle\kern 3.01991pt\hbox{$\scriptscriptstyle\left>\vphantom{p}\right.$}\kern-3.01991pt\hbox to0.0pt{\hss\hbox{$\scriptscriptstyle\left>\vphantom{p}\right.$}}\kern 3.01991pt}}}}} relates the cooling time $\beta$ and $\alpha$ as

[TABLE]

where $\Gamma_{1}=5/3$ is substituted. Therefore, $\overline{\beta}_{\text{mid}}\sim\Sigma^{0}$ readily leads to $\alpha\sim\Sigma^{0}$. Finally, the scaling of stress work (or mass accretion rate) can be derived as

[TABLE]

Note that the scalings shown in the above are obtained from our simulations at a single radius of $a=50$ AU. To see whether our scalings are universal, we need to perform simulations at different $\Omega$s (i.e. different radii) (Hirose & Shi, in preparation).

4.3 Locality of angular momentum transport and shearing box

Here, we examine the locality of angular momentum transport in the gravito-turbulence in our simulations. As pointed out by Balbus & Papaloizou (1999), the shearing boundary conditions automatically exclude the global angular momentum transport by self-gravity. Therefore, to be self-consistent, the locality of angular momentum transport must be satisfied.

Physically, the locality of angular momentum transport is determined by the disc mass $M_{\text{disc}}$ in terms of the stellar mass $M_{*}$, which was studied by global simulations using the $\beta$ cooling prescription or more realistic cooling. Generally, the local description of angular momentum transport is valid when $M_{\text{disc}}/M_{*}\lesssim 0.25$ (Lodato & Rice, 2004; Boley et al., 2006; Cossins et al., 2009) while the non-local transport cannot be negligible when $M_{\text{disc}}/M_{*}\gtrsim 0.1$ (Cai et al., 2008; Harsono et al., 2011) and becomes significant when $M_{\text{disc}}/M_{*}\gtrsim 0.5$ (Lodato & Rice, 2005; Forgan et al., 2010). The disc mass does not appear in local shearing box simulations since they are independent from the global disc model. Here, we evaluate corresponding disc mass from the surface density $\Sigma$ simply as $M_{\text{disk}}\approx\pi a^{2}\Sigma$, as often done in the literature (e.g. equation (2) in Kratter & Lodato, 2016). Then, the ratio of the disc to stellar mass is $0.088$ in the fiducial case ($\Sigma=100$ gcm*-2*) and 0.27 in the maximum surface density case ($\Sigma=300$ gcm*-2*). Given the ambiguity in evaluating $M_{\text{disc}}$ from $\Sigma$, the non-local angular momentum transport might not be negligible in most of our cases and could be significant in the largest $\Sigma$ cases.

Next, we numerically examine the locality of angular momentum transport by computing a 3-dimensional auto-correlation function of density fluctuation (e.g. Gammie, 2001),

[TABLE]

where $\delta\rho(x,y,z)=\rho(x,y,z)-\left<\rho\right>(z)$. Figure 19 shows the midplane slice $\xi(x,y,z=0)$ for the fiducial run ($L_{x}/H=24$) as well as one-dimensional cuts along the long and short axes (such as shown in the plot of the midplane slice) for $L_{x}/H=12$, 18, 24, and 32. If we look at the cuts along the short axis, the locality looks fairly good except $L_{x}/H=12$. On the other hand, along the long axis, convergence cannot be seen although the locality improves as $L_{x}$ increases (in the sense that the minimum value decreases). Apparently, we need even a larger box than $L_{x}/H=32$ to assure the locality of angular momentum transport.

However, we also have another limitation on the box size, which is the validity of using the local shearing box. In the fiducial run, the box, whose half size is $L_{x}/2=12H=38$ AU, is located at the distance, $a=50$ AU, from the central star. Therefore, increasing the box size even larger than the fiducial run ($L_{x}/H=24$) makes no sense in terms of the validity of using the local shearing box.

We can see this problem more quantitatively by taking the ratio of the wavelength of the fastest growing mode of axisymmetric GI $\lambda=2\pi Qc_{\text{s}}/\Omega$ and the distance from the central star $a$, which scales as

[TABLE]

We would like to have this ratio as small as possible so that $\lambda\ll L_{x}\ll a$ is satisfied. Then, we can take a box whose size $L_{x}$ is much larger than $\lambda$ (so that the locality of angular momentum transport is satisfactory), but still smaller than the distance from the star $a$ (so that the locality of the shearing box is guaranteed). To do that, we need to make both $T$ and $a$ small. However, it would be difficult to lower the ratio $\lambda/a$ less than order of $0.1$ because the ratio only weakly depends on the two parameters (the exponents are both $\frac{1}{2}$) and they are usually anti-correlated. The bottom line is that the parameters that we chose for the fiducial run were the best we could do for both the locality of angular momentum transport and the validity of the shearing box. This problem was not explicitly mentioned in the previous studies of gravito-turbulence since they usually used dimensionless quantities or ignoring the irradiation.

5 Summary

We explored the gravito-turbulence in irradiated protoplanetary discs using radiation hydrodynamics simulations. We used a stratified shearing box located at the radius of $a=50$ AU, which is irradiated by the central star of $T_{*}=4000$ K, $M_{*}=1M_{\odot}$, and $R_{*}=1R_{\odot}$. Under these conditions, we found that:

Gravito-turbulence is sustained for a finite range of the surface density, $80\leq\Sigma\leq 250$ gcm*-2*.The flow is laminar below the range and fragmentation occurs above the range. 2. 2.

In the regime of gravito-turbulence, the Toomre parameter $\overline{Q}_{\text{mid}}$ decreases monotonically from $\sim 1$ to $\sim 0.7$ as the surface density increases while an effective cooling time is rather constant at $\overline{\beta}_{\text{mid}}\sim$4. 3. 3.

In the gravito-turbulence, successive collisions of turbulent density waves contribute to hydrostatic balance via the dynamic pressure as well as thermal balance via the advection cooling. The turbulence dissipates through both shock heating and compressional heating. 4. 4.

The irradiation heating does not affect much the gravito-turbulence of the main body of the disc unless the grazing angle is as large as $\theta\geq 0.32$, where the turbulence is orders of magnitude smaller in terms of $\alpha$. 5. 5.

The simple cooling function used in the previous studies of gravito-turbulence does not approximate the realistic radiative transfer in the irradiated discs, regardless of the cooling time $\beta$.

The point (ii) indicates, for a fixed radius in protoplanetary discs, that there is a minimum Toomre parameter that can sustain the gravito-turbulence and that fragmentation is determined by the Toomre parameter (c.f. Takahashi et al., 2016). As discussed in Kratter & Lodato (2016), fragmentation may be driven by either cooling (by decreasing the temperature) or accretion (by increasing the surface density). The criterion found here, which is based on the surface density dependence, may apply to the latter type of fragmentation. On the other hand, the cooling time $\beta$ is expected to strongly depend on the radius (Clarke, 2009; Clarke & Lodato, 2009; Paardekooper, 2012). Therefore, studies on the radial dependence will be required to address the criterion for the cooling-driven fragmentation criterion.

We also discussed locality of angular momentum transport as well as locality of the simulation box in our simulations. The fact that both are not quite satisfactory indicates limitations on using shearing box simulations for gravito-turbulence in realistic protoplanetary discs.

Acknowledgments

We thank the anonymous referee for his/her useful comments to improve the manuscript. We thank Kengo Tomida, who kindly provided us with the EOS tables that were used in our simulations. Numerical calculations were carried out partly on Cray XC30 at CfCA, National Astronomical Observatory of Japan, and partly on Cray XC40 at YITP in Kyoto University. SH was supported by Japan JSPS KAKENH 15K05040 and the joint research project of ILE, Osaka University. JS was supported in part by the National Science Foundation under grant PHY-1144374, "A Max-Planck/Princeton Research Center for Plasma Physics" and grant PHY-0821899, "Center for Magnetic Self-Organization".

Appendix A Irradiation heating rate

In this section, we describe how we compute the irradiation heating rate $q_{\text{irr}}$ by a simple ray tracing method. We inject a single ray into every single cell located at the top and bottom surfaces of the simulation box along the direction $\bm{\Omega}=(\cos\theta,0,\mp\sin\theta)$. A single ray represents the irradiation flux entering the box through the surface of a single cell. Then the initial energy flux [erg s*-1*] of each ray is $F_{\text{irr}}\sin\theta\Delta S$, where $F_{\text{irr}}=\sigma_{\text{B}}T_{*}^{4}({R_{*}}/{a})^{2}$ and $\Delta S$ is the surface area of the cell.

As a ray travels through the box, it is attenuated by absorption, which in turn heats the gas.888We assume that a ray that escapes from the box through a radial boundary re-enters the box through the opposite side of the box. When it enters $n$-th cell, its energy flux [erg s*-1*] is written as

[TABLE]

where $\Delta l^{n^{\prime}}$ is the path length within the $n^{\prime}$-th cell, and \kappa_{\text{P{}{*}}}^{n^{\prime}}=\kappa_{\text{P{}{}}}^{n^{\prime}}(\rho,T) is the star-temperature Planck-mean opacity in the cell. Here, \rho\kappa_{\text{P{}_{}}} is fixed for simplicity. The rate at which the gas in the cell is heated by the ray is then computed as

[TABLE]

where $\Delta z$ is the thickness of the cell. The (total) irradiation heating rate $q_{\text{irra}}$ in the cell is computed as a sum of contributions by the all rays considered.

Appendix B Equation of self-gravitational energy

The self-gravitational energy equation is written as (e.g. Balbus & Papaloizou, 1999)

[TABLE]

or equivalently

[TABLE]

In a steady state, the time and horizontally-averaged version of the equation is written as

[TABLE]

Applying the shearing periodic boundary conditions, the first term in RHS is rewritten as (see Balbus & Papaloizou, 1999)

[TABLE]

where, $|_{x^{\pm}}$ denotes averaging in the $y$ direction at each height on the plane of $x=x_{\pm}$. The last equality would be validated when the box size is large enough compared with the typical length of the turbulence.

Thus, the equation 41 is written as

[TABLE]

Appendix C Evolution of the gravitationally bounded clump

In this section, we describe the rapid collapse of a gravitationally-bounded clump observed in the case of $\Sigma=300$ g cm*-2*. In figure 20, we plotted the evolutionary track of $(T(t),\rho(t))$ of the cell at the centre of the clump, on the colour contour of the Rosseland-mean opacity $\kappa_{\text{R}}(T,\rho)$ as well as on that of the adiabatic exponent $\Gamma_{1}(T,\rho)$, from $t=2.07$ to $t=2.39$ orbits. The track evolves from the lower left (low $T$ and low $\rho$) to the upper right (high $T$ and high $\rho$). As seen from intervals between the marks, which are placed for every 0.01 orbits on the track, the density increase is accelerated when the adiabatic exponent $\Gamma_{1}$ decreases from $4/3$ to $\sim 1.1$ due to dissociation of H2 molecules. This is exactly the same physical process as the first core collapse in the star formation (e.g. Tomida et al., 2013).

We can learn more about the collapse by comparing the two extra tracks plotted in white in the same figure; one is $(T(t),p(t)\times 10^{-12})$ and the other is $(T(t),\rho(t)^{\Gamma_{1}(t)})$. We see that the pressure increase is actually much smaller than that expected in the adiabatic evolution ($\sim\rho^{\Gamma_{1}}$) for $3.1\leq\log T\leq 3.25$. This means that the cooling is very effective for that temperature range. Since the beginning of the range $\log T\sim 3.1$ exactly corresponds to the dust sublimation temperature, the effective cooling should come from large reduction of Rosseland-mean opacities due to the dust sublimation. Therefore, the rapid collapse of the clump was first triggered by the dust sublimation, and then followed up by the dissociation of H2 molecules.