text:
# Thermo-Hydrodynamics of Circumstellar Disks with High-mass Planets1footnote 11footnote 1To appear in The Astrophysical Journal (v599 n1 December 10, 2003 issue). Also available as ApJ preprint doi: 10.1086/379224.

## 1 Introduction

Since the very first attempt of modeling a protoplanet interacting with its primitive nebula by Miki (1982), two decades have passed. Seventeen years had to go by for tackling this problem with the strength of new and more powerful computational tools (Bryden et al., 1999; Kley, 1999; Lubow, Seibert, & Artymowicz, 1999; Miyoshi et al., 1999). As a consequence of the mounting interest in extrasolar planets, boosted by eight years of uninterrupted new detections333Data archives on extra-solar planets, with the latest information, can be found at the Extrasolar Planet Encyclopedia (http://www.obspm.fr/planets) and at the California & Carnegie Planet Search (http://exoplanets.org/)., the scientific community has made an extraordinary effort, over the past four years, to analyze disk-planet interactions by means of numerical calculations.

Yet, regardless of the processes accounted for in the implemented numerical models, limitations and restrictions remain an issue to deal with. As a matter of fact, numerical calculations still lag behind the present knowledge of the possible physical effects which occur in a circumstellar environment containing massive bodies. Thereupon, simplifying assumptions are always demanded. Nonetheless, steps forward have been made in many directions to reduce their number.

Three-dimensional (3D) disks containing Jupiter-size objects have been investigated by Kley, D’Angelo, & Henning (2001) and, including also magnetic fields, by Papaloizou & Nelson (2003), Nelson & Papaloizou (2003), and Winters, Balbus, & Hawley (2003). D’Angelo, Henning, & Kley (2002, hereafter Paper I) carried out global two-dimensional (2D) high-resolution simulations of Jupiter-mass as well as Earth-mass planets, resolving also the flow structure within the planetary’s Roche lobe by applying a resolution enhancement strategy. This study has been lately extended to the 3D space by D’Angelo, Kley, & Henning (2003, hereafter Paper II) and by Bate et al. (2003), who used instead a single grid with variable step size. Nelson & Benz (2003a, b) performed 2D computations, considering the effects due to disk self-gravity. Morohoshi & Tanaka (2003) have investigated the interactions between a planet and an optically thin disk with local simulations in the shearing sheet approximation, in order to determine the one-sided gravitational torque acting on the planet.

With the present work we intend to relax the widely adopted local-isothermal approximation, i.e., that of treating a disk as a system having a fixed temperature distribution which depends neither on time nor on any other system variable. Yet, we will preserve the global viewpoint which we believe to be crucial because disk-planet interactions can be highly non-local.

Many of the sophisticated accretion disk studies adopt 2D or 1+1D schemes due to the difficulty of dealing with all the necessary ingredients in a full 3D space. Until very recently, in fact, radiation (e.g., Dullemond et al., 2002), convection (e.g., Agol et al., 2001), magnetic fields (e.g., Matt et al., 2002) , mixtures of gas and dust (e.g., Suttner & Yorke, 2001), and chemical evolution (e.g., Markwick et al., 2002) have been treated via two-dimensional models. Furthermore, numerical simulations generally involve parameter studies, because many physical quantities are poorly known. Thereby, with a restricted number of dimensions, one is at least able to acquire a spectrum of possible approximate solutions.

In order to overcome the local-isothermal hypothesis and to maintain a global approach at the same time, we introduce an energy transport in 2D ($`r`$$`\phi `$) models. We include all the supposedly major causes responsible for generation, transfer, and loss of energy in low-temperature circumstellar disk environments. We take advantage of the small aspect ratio of the disk, and assume that all the radiation transport is effective only in the vertical direction. Such approximation works reasonably well in accretion disks, away from the boundary and surface layers (see, e.g., Pringle, 1981). Though such an approach relies on the “slimness” assumption, it is as well appropriate in the local environment around a protoplanet. In fact, around there, the disk becomes even thinner when the planet’s gravitational action is accounted for. Hence, the amount of energy transported by radiation in the vertical direction still overwhelms that transported horizontally. Although this represents a simplistic kind of thermal description, it allows to investigate processes which have not been considered thus far, i.e., the joint thermo-hydrodynamical evolution of disk-planet interactions. Furthermore, since these simulations benefit of a nested-grid numerical technique (Paper I, Paper II), though relying on the global approach we are capable of achieving resolutions large enough to accurately describe the system within the planetary’s Roche lobe.

The outline of the paper is the following. In the next section we introduce the physical formulation of the problem, focusing on the equation for the thermal energy density and related issues. Section 3 concerns the numerical method utilized to solve the energy equation along with a series of tests. In § 4 we specify the adopted parameters and other various details. Sections 5 and 6 describe, respectively, global and local properties of models with different masses and viscosity regimes. In § 7 we report the results for mass accretion and orbital migration. Our conclusions are presented in § 8.

## 2 Physical Formalism

As mentioned before, we are going to tackle the problem of energy transfer in an accretion disk with an embedded planet. Since solving the full set of equations in three dimensions and simultaneously achieving sufficient resolution around the protoplanet is computationally difficult at the moment, we restrict to two-dimensional disk models.

The reason for this choice is two-folded. On the one hand, with affordable computational times and sufficiently high resolutions, it allows to investigate the coupling between hydrodynamics and thermodynamics both in the disk and nearby the embedded body. On the other, it permits to adopt a strategic assumption in writing the energy equation, i.e., that the horizontal energy transport can be neglected compared to the vertical one. In other words, in a 2D disk one can assume that radiation transport is only effective along the direction perpendicular to the equatorial plane of the system.

Hence, let us describe the disk material through the Navier-Stokes equations for the surface density $`\mathrm{\Sigma }`$, the linear momentum (density) $`\mathrm{\Sigma }u_r`$, and the angular momentum (density) $`r\mathrm{\Sigma }u_\phi `$, where we choose as working coordinate frame a cylindrical one $`\{O;r,\phi ,z\}`$. The origin is located at the center of mass of the star and the planet while the disk lies in the plane $`z=0`$. The hydrodynamics equations have been already stated, by using the same notations, in Paper I.

As in most of the computations performed so far, we avoid dealing with the Poisson equation and rather assume that the gravitational field is determined only by the star and the planet, neglecting the self-gravity of the disk material. This is allowed by the circumstance that the Toomre $`Q`$-parameter is well beyond $`1`$ in all of our models. Nelson & Benz (2003a, b) found that some effects indeed arise when disk self-gravity is taken into account. However, we note that they consider disk masses ten times as large as those considered here.

Thereby, we rely on a smoothed point-mass gravitational potential, which well reproduces more realistic protoplanetary gravitational potentials (Paper II),

$$\mathrm{\Phi }=\frac{GM_{}}{\left|𝒓𝒓_{}\right|}\frac{GM_\mathrm{p}}{\sqrt{\left|𝒓𝒓_\mathrm{p}\right|^2+\delta ^2}},$$
(1)
where $`𝒓_{}`$ and $`𝒓_\mathrm{p}`$ are the radius vectors indicating the positions of the star and the planet, respectively. As done in Paper II, we set the parameter $`\delta `$ as a constant length (see § 4 for details).

The thermal structure of the fluid has been usually fixed through a power law of the distance $`r`$. Thereupon, the set of equations has been closed by adding an equation of state which connects the local gas (two-dimensional) pressure to the surface density by means of a local isothermal sound speed $`c_s^{\mathrm{iso}}`$

$$c_s^{\mathrm{iso}}=h\sqrt{\frac{GM_{}}{r}}.$$
(2)
Thus, the Mach number of the flow is considered constant throughout the system’s evolution, constrained solely by the disk aspect ratio $`h=H/r`$ (e.g., Lubow, Seibert, & Artymowicz, 1999; Miyoshi et al., 1999; Nelson et al., 2000; Tanaka, Takeuchi, & Ward, 2002; Masset, 2002; Nelson & Papaloizou, 2003; Bate et al., 2003; Nelson & Benz, 2003b). Alternatively, Kley (1999) and Tanigawa & Watanabe (2002) did some simulations using a polytropic equation of state of the type $`P=\overline{K}\mathrm{\Sigma }^\gamma `$. However, the evolution of the thermal properties of the system have not been generally considered.

In the present calculations we use an ideal equation of state which directly ties the gas pressure $`P`$ to the thermal energy density (energy per unit area) $`E`$

$$P=\left(\gamma 1\right)E,$$
(3)
where the adiabatic index $`\gamma `$ is a constant. Furthermore, if we suppose that disk material behaves as a perfect gas, then the temperature is

$$T=\left(\frac{\mu m_\mathrm{H}}{k}\right)\frac{P}{\mathrm{\Sigma }},$$
(4)
where $`\mu `$ is the mean molecular weight, $`m_\mathrm{H}`$ is the hydrogen mass, and $`k`$ the Boltzmann constant. The adiabatic sound speed is given by

$$c_s=\sqrt{\gamma \left(\frac{kT}{\mu m_\mathrm{H}}\right)}.$$
(5)
Equation (3) implicitly assumes that radiation pressure $`P_{\mathrm{rad}}`$ is negligibly small compared to gas pressure. Such hypothesis is connected to the relatively low temperature regimes we deal with (in fact $`P_{\mathrm{rad}}T^4`$) and its validity was checked afterward. In all of the models presented here, the ratio $`P_{\mathrm{rad}}/P_{\mathrm{gas}}`$ never exceeds $`10^4`$.

### 2.1 Energy Equation

Equations for energy transport differ according to the processes that have to be included for a correct description of the energy budget of fluid elements. In our case, we suppose that a parcel of disk fluid can gain or lose thermal energy only because of flow advection, compressional work, viscous dissipation, and dissipative effects due to radiation transport. In this sense, rather than strictly treating the transfer of radiation in the disk, we will account only for the cooling effects that radiation causes. Then, the energy equation takes the following form

$$\frac{E}{t}+\left(E𝒖\right)=P𝒖+\mathrm{{\rm Y}}\mathrm{\Lambda }.$$
(6)
In equation (6) we indicated with $`\mathrm{{\rm Y}}`$ the dissipation function and with $`\mathrm{\Lambda }`$ the radiated energy. In two dimensions, either of these quantities is an energy flux. Both functions are always positive, thereby the former acts as a heating source, whereas the latter as a cooling term. In this study we do not consider convection because we do not expect it to be a major energy source (D’Alessio et al., 1998). Furthermore, we neglect irradiation from the central star which is presumably an effective heating mechanism only in the upper layers of circumstellar disks (see, e.g., D’Alessio et al., 1998), whose effects however are reduced by small disk scale heights. More importantly, as a consequence of the gap formation and the protoplanet’s gravity, circumplanetary material is shielded from stellar radiation, as it will be clarified later.

The function $`\mathrm{{\rm Y}}`$ can be directly computed from the components of the viscous stress tensor $`S_{rr}`$, $`S_{\phi \phi }`$, and $`S_{r\phi }`$ (Mihalas & Weibel Mihalas, 1999; Collins, Helfer, & van Horn, 1998)

$$\mathrm{{\rm Y}}=\frac{1}{2\nu \mathrm{\Sigma }}\left(S_{rr}^2+S_{\phi \phi }^2+S_{r\phi }^2\right)+\frac{2\nu \mathrm{\Sigma }}{9}\left(𝒖\right)^2.$$
(7)
The divergence term $`𝒖`$ in the equation arises because a three-dimensional definition of stress tensor S is adopted. In fact, the reduction from three to two dimensions is merely achieved by setting the vertical component of the flow field $`u_z`$ and the derivative operator $`/z`$ to zero.

In a three-dimensional disk, the last term on the right-hand side of equation (6) is equal to $`𝑭`$, where $`𝑭`$ is the frequency-integrated radiation flux

$$𝑭=\frac{16\sigma _\mathrm{R}}{3\kappa \rho }T^3T.$$
(8)
In the above equation $`\sigma _\mathrm{R}`$ is the Stefan-Boltzmann constant, $`\kappa `$ is a frequency-integrated opacity coefficient and $`\rho `$ is the mass density. As mentioned before, we suppose that the amount of energy transported by radiation in the vertical direction is much larger than that transported horizontally, i.e., $`\left|F_r\right|`$ and $`\left|F_\phi \right|\left|F_z\right|`$. The validity of this statement holds as long as the vertical extent of the disk remains very small compared to the disk characteristic size in the other directions. Thus, in a two-dimensional cylindrical disk we have

$$\mathrm{\Lambda }=_{\mathrm{}}^+\mathrm{}𝑭dz_{\mathrm{}}^+\mathrm{}\frac{F_z}{z}dz.$$
(9)
Since the vertical disk structure is not meant to be accounted for, it is possible to equate the pressure scale height to the photospheric scale height and assume that all the radiation is liberated at $`z=\pm H`$. Thereupon, equation (9) becomes

$$\mathrm{\Lambda }=_H^{+H}\frac{F_z}{z}dz=F\left(H\right)F\left(H\right)=2F\left(H\right).$$
(10)
In 2D disks it is useful to measure the emitted flux by means of an effective temperature: $`F\left(H\right)=\sigma _\mathrm{R}T_{\mathrm{eff}}^4`$, therefore

$$\mathrm{\Lambda }=2\sigma _\mathrm{R}T_{\mathrm{eff}}^4.$$
(11)
The factor $`2`$ indicates that in a disk radiation escapes from both of its sides.

A simple relation between the midplane temperature and the emergent radiation flux can be found by writing

$`\mathrm{\Lambda }`$ $`=`$ $`2{\displaystyle _0^H}{\displaystyle \frac{F_z}{z}}dz`$ (12)
$`=`$ $`2\left[F\left(H\right)F\left(0\right)\right]`$
$``$ $`{\displaystyle \frac{8}{3}}{\displaystyle \frac{\sigma _\mathrm{R}}{\kappa \rho H}}\left[T^4\left(0\right)T^4\left(H\right)\right].`$
In the inner parts of a circumstellar disk the inequality $`T^4\left(0\right)T^4\left(H\right)`$ generally holds (e.g., Bell et al., 1997; D’Alessio et al., 1998), hence equation (12) yields

$$\mathrm{\Lambda }\frac{2\sigma _\mathrm{R}T^4}{\left(3/4\right)\tau },$$
(13)
in which we introduced the optical thickness $`\tau =\kappa \rho H=\frac{1}{2}\kappa \mathrm{\Sigma }`$ (we generally adopt a Rosseland mean opacity). We notice that with the previous definition the total disk optical thickness is $`2\tau `$. From now on, the quantity $`T`$ refers to the disk midplane temperature.

Equation (13) represents a fairly good approximation when the medium is very optically thick, i.e., $`\tau 1`$. This is indeed the case in those regions of unperturbed accretion disks which we simulate, because $`\mathrm{\Sigma }`$ is large enough. Yet, because of the action of massive bodies, in our case deep density gaps form where material is very diluted. In addition there are zones, where the disk spirals interact with the circumplanetary disk spirals, which can have very low densities (see Paper I, Fig. 12 and 13). Such conditions make equation (13) not always applicable. Hubeny (1990) found a more rigorous relation between the effective and midplane accretion disk temperature, which represents a generalization of the gray model of classical stellar atmospheres in local thermodynamic equilibrium. Following Hubeny’s theory, in a circumstellar disk we can write

$$\mathrm{\Lambda }=2\sigma _\mathrm{R}T^4\left[\frac{3\tau }{8}+\frac{\sqrt{3}}{4}+\frac{\epsilon _\mathrm{H}}{4\tau }\right]^1.$$
(14)
According to equation (14), the emitted flux is inversely proportional to $`\tau `$ in optically thick disk portions, whereas it becomes proportional to $`\tau `$ in the opposite limit. The quantity $`\epsilon _\mathrm{H}`$ is equal to the ratio between the Rosseland and the Planck mean opacities. This quantity can be also roughly interpreted as the ratio of total extinction (i.e., absorption plus scattering coefficient) to the pure absorption coefficient: $`\kappa _{\mathrm{ext}}=\kappa _{\mathrm{abs}}+\kappa _{\mathrm{sc}}`$. In these computations, we set $`\epsilon _\mathrm{H}=1`$ because we checked that we can neglect radiation scattering (see § 3.1). Thus, no distinction is made between the extinction and the absorption coefficients.

Equation (14) mimics the radiative losses of both optically thick and optically thin media, and therefore is very suitable to our purposes. In fact, it has been successfully applied to many studies on accretion disks (Popham et al., 1993; Godon, 1996; Collins, Helfer, & van Horn, 1998; Huré et al., 2001).

### 2.2 Opacity Table

In order to calculate the disk optical semi-thickness $`\tau `$, we adopt the opacity formulas derived by Bell & Lin (1994) (see Fig. 1), as an improvement of those by Lin & Papaloizou (1985). Eight temperature regimes are identified according to the dominant processes active in each of them. Contributions from dust grains, molecules, atoms, and ions are accounted for. Since we simulate a distance range where disk material can be relatively cold, dust opacity is crucial for an accurate estimation of radiative losses (Lin, 1981). Because of this, Bell’s opacity includes grain absorption as tabulated by Alexander, Augason, & Johnson (1989).

However, for the sake of comparison we ran some test cases (see § 3.1) with a new opacity coefficient developed by Semenov and collaborators (2002, private communication) and based on the improved grain opacity tables provided by Henning & Stognienko (1996). They aimed at coupling gas and dust opacities, focusing on the temperature range proper to protostellar disk environments.

In either case, the midplane temperature $`T`$ and mass density $`\rho `$ have to be provided in cgs units. In turn $`\kappa =\kappa (T,\rho )`$, has units of square centimeters per gram. Within the infinitesimally thin disk approximation, dynamics variables are vertically integrated. Thus, the mass density $`\rho `$ is not directly available. It must be retrieved from the surface density and the disk semi-thickness $`H`$. For this purpose, we assume that

$$\rho =\frac{\mathrm{\Sigma }}{2H}.$$
(15)
In the next section we explain how the disk scale height $`H`$ is calculated in a fashion such to account also for the gravitational influence of the planet. This guarantees a more refined and consistent modeling of the system.

### 2.3 Disk Scale Height

Inside a circumstellar disk it is natural to assume that material is in hydrostatic equilibrium along the vertical direction (Pringle, 1981; Frank, King, & Raine, 1992). Let $`g_z`$ be the gravitational acceleration in the $`z`$-direction within the disk, then the hydrostatic equilibrium reads

$$\frac{P}{\mathrm{\Sigma }}=2_0^Hg_zdz.$$
(16)
The factor $`2`$ is necessary because the midplane pressure sustains both sides of the disk.

When the temperature distribution is not coupled to the rest of the dynamical variables, equation (16) is used to obtain the disk pressure from the density and the pressure scale height, assuming that $`g_z=GM_{}/\left(r^2+z^2\right)`$ and omitting the planet gravitational attraction.

In the present study, we compute the pressure directly from the thermodynamical processes occurring in the gas (eq. ), thus we can use equation (16) to estimate the disk scale height $`H`$ in a consistent fashion. This quantity, in fact, is needed to obtain the mass density from the surface density, as explained above. Since the left-hand side of equation (16) is equal to $`c_s^2/\gamma `$, we can integrate the right-hand side in order to get an implicit function of $`H`$. Furthermore, the effects due to the gravitational field of the protoplanet can be included in this evaluation. If we set $`q=M_{}/M_\mathrm{p}`$ and $`s=\left|𝒓𝒓_\mathrm{p}\right|`$, equation (16) generates the implicit function

$$\frac{H^2}{2r^2}\frac{q}{\sqrt{\left(s/r\right)^2+\left(H/r\right)^2}}=\frac{r}{2\gamma }\left(\frac{c_s^2}{GM_{}}\right)q\left(\frac{r}{s}\right).$$
(17)
Although there are no constraints on the ratio $`s/r`$, because with nested grid computations that ratio can be very small, the aspect disk ratio $`H/r`$ is smaller than one, by working hypothesis, otherwise the 2D geometry would not be consistent. Therefore, the second term on the left-hand side can be expanded in a binomial series in $`H/r`$ (up to the its second power), obtaining $`H`$ as an explicit function

$$H^2=\frac{r}{\gamma }\left(\frac{c_s^2}{GM_{}}\right)\left[\frac{1}{r^2}+q\left(\frac{r}{s^3}\right)\right]^1.$$
(18)
We note that, in the limit $`q0`$ or $`sr`$, equation (18) yields $`\sqrt{\gamma }H=c_s/\mathrm{\Omega }_\mathrm{K}`$, as in regular accretion disks. Close to the planet, i.e., in the circumplanetary disk where $`sr`$, the above relation reduces to $`\sqrt{\gamma }H=c_s/\sqrt{GM_\mathrm{p}/s^3}`$, which resembles the previous limit because the circumstellar Keplerian angular velocity is replaced by the circumplanetary one. Rearranging the terms in equation (18), a more compact expression can be written

$$H^2=\frac{1}{\gamma }\left(\frac{c_s}{\mathrm{\Omega }_\mathrm{K}}\right)^2\left[1+q\left(\frac{r}{s}\right)^3\right]^1.$$
(19)
Equation (19) has a singularity at $`s=0`$. This arises from the corresponding singularity in the gravitational potential $`\mathrm{\Phi }`$, which is overcome by introducing the smoothing length (eq. ). In our computations, for consistency reasons, the distance $`\sqrt{s^2+\delta ^2}`$ substitutes $`s`$ in equation (19).

### 2.4 Artificial Viscosity

Non-linear effects in wave propagation inevitably lead to shock formation. This is indeed the case in disk-planet simulations. In ideal fluids shocks are mathematical discontinuities. Therefore, in finite-differencing schemes, they must extend only over one or two grid points. In order to provide the correct jump conditions, ahead and behind a shock front, dissipative terms have to be present in the equations. This is usually done by introducing a non-linear viscous pressure otherwise known as artificial viscosity. Since the pressure gradient is no longer proportional to the density gradient, shocks stronger than those observed in local-isothermal computations (see, e.g., Paper I) may develop. Therefore, the physical viscosity might not be sufficient to provide the correct jump conditions across shock fronts.

The most rigorous treatment of shocks in a multidimensional space, with a generic metric tensor, requires the definition of an isotropic viscous stress tensor Q (Winkler & Norman, 1986), which can be written as (Mihalas & Weibel Mihalas, 1999; Stone & Norman, 1992)

$$\text{Q}=\mu _\mathrm{Q}\left[𝒖\frac{1}{3}\left(𝒖\right)\text{I}\right],$$
(20)
being I the unit tensor. Here we choose to discard the off-diagonal, i.e., shear tensor components because they may lead to artificial angular momentum transport (Stone & Norman, 1992). The coefficient of artificial viscosity is defined as

$$\mu _\mathrm{Q}=^2\mathrm{\Sigma }\mathrm{min}(𝒖,0),$$
(21)
where $``$ represents the length over which the shock is smeared out. This is usually fixed to the maximum grid spacing. The coefficient $`\mu _\mathrm{Q}`$ is positive only for compression and zero for expansion, so the artificial viscosity is large in shocks and negligibly small elsewhere.

Since the artificial viscous tensor acts as a pseudo-pressure, it affects both momentum and energy equations through the terms $`\text{Q}`$ and $`\text{Q}𝒖`$ (meant as a tensor product), respectively. The aforementioned equations are updated, by using the correct components $`Q_{ij}`$, as explained in Stone & Norman (1992). One severe side-effect is that Q can reduce the viscous time step by a factor which is proportional to $`\left[^2𝒖\right]^1`$.

## 3 Energy Equation Solver

The general numerical method employed to solve the hydrodynamical equations on a hierarchy of nested grids, applied to simulations of disk-planet interactions, has been explained in Paper I. In contrast to previous calculations, thermal energy density $`E`$ now appears as an independent variable. Here we outline how the source terms (right-hand side of eq. ) are dealt with. In the framework of the nested-grid scheme, whenever required, $`E`$ is interpolated from a finer to a coarser grid and vice versa, according to the procedures outlined in Paper I (see also Paper II).

The energy equation (eq. ) is solved by means of a multi-step operator splitting method. The first step takes care of the energy advection and this is done in the same fashion as the advection of the other dynamical quantities, i.e., through the van Leer’s algorithm. Viscous dissipation and radiative cooling are treated by means of a predictor-corrector scheme, which is second-order accurate in time, where heating and cooling terms are considered simultaneously for better accuracy and stability.

During the course of the first experiments, it was discovered that, due to very high pressures and large values of the flow divergence (in the wake of the circumplanetary spirals) the compressional work could lower the thermal energy by a large amount. Consequently, the predictor-corrector procedure applied to this equation term occasionally occurred to be unstable, producing negative energy values. Therefore we decided to take advantage of equation (3) and use an analytic solution for updating the energy.

In the framework of the operator-splitting approach, in order to correct the energy because of the gas compression/dilation, we have to solve numerically the equation

$$\frac{E}{t}=P𝒖.$$
(22)
With the aid of the state equation, this becomes

$$\frac{E}{t}=\left(\gamma 1\right)E𝒖.$$
(23)
Equation (23) has an analytic solution, given by $`E=E_0\mathrm{exp}\left[\left(\gamma 1\right)\mathrm{\Delta }t𝒖\right]`$, where $`\mathrm{\Delta }t`$ is a time interval over which the divergence $`𝒖`$ has not appreciably changed. Hence, thermal energy density can be updated as follows:

$$E_{n+1}=E_n\mathrm{exp}\left[\left(\gamma 1\right)\mathrm{\Delta }t_n𝒖_n\right].$$
(24)
We notice that, expanding the exponential of equation (24) in a Taylor series and keeping all the terms up to the second order in $`\mathrm{\Delta }t_n`$, one finds

$`{\displaystyle \frac{E_{n+1}E_n}{\mathrm{\Delta }t_n}}`$ $`=`$ $`P_n𝒖_n`$ (25)
$`+`$ $`{\displaystyle \frac{1}{2}}\left(\gamma 1\right)P_n\mathrm{\Delta }t_n\left(𝒖_n\right)^2,`$
where $`P_n=\left(\gamma 1\right)E_n`$. It is easy to prove that equation (25) is equivalent to a predictor-corrector scheme. The above equation also shows that the second-order correction always increase the thermal energy.

Although the procedure in equation (24) is not so general as the predictor-corrector algorithm, it is more accurate and unconditionally stable. Furthermore, contrary to what may seem at first glace, from the computational viewpoint it is faster than the other because it requires less mathematical operations.

### 3.1 Solver Test

We have tested the equation energy solver from both the numerical and the physical point of view. Here we present some of these tests.

The first test is intended to check whether the equation of energy furnishes physically consistent results by comparing the computed temperature to that derived from some analytical solutions of equation (6). If we deal with a stationary Keplerian disk then the energy budget simplifies enormously. In fact, in a pure Keplerian flow both energy advection and compressional work are negligibly small. Therefore equation (6) reduces to (e.g., see Pringle, 1981)

$$\mathrm{{\rm Y}}\mathrm{\Lambda }=0.$$
(26)
From the same hypotheses, it follows that the dissipation function can be written in a simple form

$$\mathrm{{\rm Y}}=\mathrm{\Sigma }\nu \left(r\frac{\mathrm{d}\mathrm{\Omega }_\mathrm{K}}{\mathrm{d}r}\right)^2=\frac{9}{4}\mathrm{\Sigma }\nu \mathrm{\Omega }_\mathrm{K}^2.$$
(27)
In fact, in a Keplerian flow, the terms $`\left(u_r/r\right)^2`$ and $`\left(u_r/r\right)^2`$ are both much less than $`\left(r\mathrm{\Omega }_\mathrm{K}/r\right)^2`$. If we additionally assume that the disk is optically thick then the emitted flux can be approximated to

$$\mathrm{\Lambda }=2\sigma _\mathrm{R}T_{\mathrm{eff}}^4\frac{2\sigma _\mathrm{R}T^4}{\left(3/8\right)\tau },$$
(28)
as implied by equation (14), when $`\tau 1`$. Recalling that $`\tau =\frac{1}{2}\kappa \mathrm{\Sigma }`$, then, equating equation (27) and equation (28) yields

$$T^4=\frac{27}{128}\frac{\kappa \nu }{\sigma _\mathrm{R}}\mathrm{\Sigma }^2\mathrm{\Omega }_\mathrm{K}^2.$$
(29)
Now we assume that disk material has an opacity which we can be cast in the form

$$\kappa =\kappa _0T^2\mathrm{cm}^2\mathrm{g}^1.$$
(30)
with $`\kappa _0=2\times 10^6`$ (Morfill, Tscharnuter, & Voelk, 1985). Placing equation (30) in equation (29), one finds

$$T=\sqrt{\frac{27}{128}\frac{\kappa _0\nu }{\sigma _\mathrm{R}}}\mathrm{\Sigma }\mathrm{\Omega }_\mathrm{K}.$$
(31)
Such expression allows a direct comparison with the temperature distribution obtained from simulations, whose setup involves equation (27), equation (28), and the opacity law (30).

In order to carry out a comparison of this kind, we simulated an unperturbed disk (i.e., without any embedded body), with borders at $`1`$ and $`20\text{AU}`$, surrounding a Solar-mass star. The disk mass is $`M_\mathrm{D}=0.03M_{}`$ and the kinematic viscosity is $`\nu =5\times 10^{16}\mathrm{cm}^2\mathrm{s}^1`$. Both the initial surface density and temperature are constant: $`\mathrm{\Sigma }\left(t=0\right)=197\mathrm{g}\mathrm{cm}^2`$ and $`T\left(t=0\right)=352\mathrm{K}`$. Once the system has reached a stationary state, the surface density is expected to decay as $`1/\sqrt{r}`$, because $`\nu `$ is constant (Lynden-Bell & Pringle, 1974). Indeed, this is what one can observe in Figure 2 (left panel), where the azimuthally averaged, computed surface density $`\mathrm{\Sigma }`$ (crosses) is fitted by the power-law

$$\mathrm{\Sigma }=300\sqrt{\frac{5\text{AU}}{r}}\mathrm{g}\mathrm{cm}^2.$$
(32)
On the other hand, equation (31) states that $`T\mathrm{\Sigma }\mathrm{\Omega }_\mathrm{K}1/r^2`$, or more precisely, $`T=104\left(5\text{AU}/r\right)^2\mathrm{K}`$. Figure 2 (right panel) shows how the calculated averaged temperature $`T`$ (crosses) fits to this analytic estimation.

For the sake of completeness, we repeated the same type of test in which, instead of the Morfill et al.’s relation (eq. ), we chose the Kramer’s law $`\kappa =6.6\times 10^{22}\rho T^{3.5}\mathrm{cm}^2\mathrm{g}^1`$, as in Frank, King, & Raine (1992), obtaining the same good agreement between analytical and numerical results.

Although not reported here, we also tested the complete algorithm against convergence stability, i.e., we checked that it always converges to the same solution whatever the initial conditions are. Two simulations were set up in which the full form of equation (6) is solved along with equation (7), equation (14), and the Bell’s opacity coefficient (§ 2.2). The disk radial domain is the same as before, but now the disk mass is smaller ($`M_\mathrm{D}=0.01M_{}`$) and so is the viscosity ($`\nu =10^{16}\mathrm{cm}^2\mathrm{s}^1`$). The initial $`\mathrm{\Sigma }`$ is imposed to be constant and equal to $`66\mathrm{g}\mathrm{cm}^2`$ for both models. As initial temperature, in one model $`T\left(t=0\right)=14.5\mathrm{K}`$, while in the other $`T\left(t=0\right)=580\mathrm{K}`$. The “cold” as well as the “hot” model soon evolve toward a stationary state. In each of the cases, the solutions exactly match, which indicates that initial values play no role in the equilibrium solution attained after the system has relaxed. The relaxation time somehow depends on the initial $`\mathrm{\Sigma }`$, whereas scarcely any dependency on $`T`$ at $`t=0`$ is observed.

Finally, we examine some implications related to the opacity tables. A simulation was performed with the same parameterization as the “hot” model mentioned above. But instead of choosing Bell’s opacity, this model is based on the new tables by Semenov and collaborators (2002, private communication). The resulting surface density distribution, at equilibrium, is hardly distinguishable from that obtained with the “hot” model. Furthermore, not much difference is measured between the stationary averaged temperature profiles, as demonstrated in Figure 3. The result is not surprising, since already Liu & Meyer-Hofmeister (1997) tested the influence of different opacity tables on the vertical structure of accretion disks and found that the change in disk structure, due to an improved opacity coefficient, is hardly perceivable when compared to the uncertainties connected with the general disk parameterization.

Hence, throughout these simulations we use the opacity coefficient by Bell & Lin (1994). It is worthy to note that this opacity table is well tested and has been extensively adopted in accretion disk studies (see, e.g., Godon, 1996; Bell et al., 1997; Klahr et al., 1999; Armitage, Clarke, & Tout, 1999; Papaloizou & Terquem, 1999; Nomura, 2002).

Equation (14) contains the quantity $`\epsilon _\mathrm{H}`$ which can be interpreted as the ratio of extinction to the absorption coefficient. However, the physical conditions of the protostellar environment are such that it is allowed to neglect radiation scattering and write $`\kappa _{\mathrm{ext}}\kappa _{\mathrm{abs}}`$, and hence $`\epsilon _\mathrm{H}1`$. In this phase, we checked such approximation by repeating the calculation addressed in Figure 3 and setting $`\epsilon _\mathrm{H}=\kappa _{\mathrm{ext}}/\kappa _{\mathrm{abs}}`$, in equation (14). The resulting temperature profile is hardly distinguishable from the thick line in Figure 3. Therefore, the hypothesis $`\epsilon _\mathrm{H}=1`$ is quite correct in our applications.

## 4 Model Parameters

Model parameterization deserves more attention in these computations than it did in the simulations where energy transfer by radiation is not taken into account. The reason for this resides in the nature of the opacity coefficient $`\kappa `$, that is always in cgs units. Hence, lengths and masses are to be fixed in physical units and consequently outcomes will not be scale-free, as they were in other circumstances (see, e.g., Paper I; Paper II).

As in previous simulations, we model a circumstellar disk, orbiting a one-solar-mass star ($`M_{}=1M_{}`$), whose radial borders are $`r_{\mathrm{min}}=2\text{AU}`$ and $`r_{\mathrm{max}}=13\text{AU}`$. The initial mass enclosed in this domain is $`M_\mathrm{D}=4.8\times 10^3M_{}`$, i.e., $`0.01M_{}`$ within $`20\text{AU}`$ (D’Alessio et al., 1998). Since in stationary Keplerian disks, the energy budget is regulated by the balance between viscous dissipation and radiative losses (eq. ), the magnitude of viscosity may play an important role. Hence, we will cover various viscosity regimes. The bulk of the simulations were run with a constant kinematic viscosity coefficient, by setting $`\nu =1.0\times 10^{15}`$, $`5.0\times 10^{15}`$, and $`1.0\times 10^{16}\mathrm{cm}^2\mathrm{s}^1`$ (Bell et al., 1997). Besides physical viscosity, we consider an artificial viscosity, which is active only in compression regions, according to the tensor formulation explained in § 2.4. The length $``$ is chosen to be equal to the maximum grid spacing, on each grid level (Stone & Norman, 1992; Ziegler, 1998). The mean molecular weight is $`\mu =2.39`$ (Morfill, Tscharnuter, & Voelk, 1985) and the adiabatic index $`\gamma =1.4`$.

Figure 4 shows characteristic averaged quantities of non-perturbed, stationary disks with the aforementioned characteristics. To distinguish among the different viscosity (and temperature) disk regimes, models will be referred to as “Hot” (H), “Warm” (W), and “Cold” (C).

To further explore the parameter space, two models were executed applying the Shakura-Sunyaev viscosity prescription $`\nu =\alpha c_sH`$, with $`\alpha =10^2`$ and $`10^3`$ (models A2 and A3). Hence, in such models $`\nu `$ is a function of space and time.

An embedded planet is placed on a fixed circular orbit at $`r_\mathrm{p}=5.2\text{AU}`$. The azimuthal position of the planet is kept fixed ($`\phi _\mathrm{p}=\pi `$) by the employment of a rotating grid. In order to obtain reliable outcomes in the two-dimensional geometry, we consider a planetary mass range extending from $`M_\mathrm{p}=33M_{}`$ to $`M_\mathrm{p}=1M_\mathrm{J}`$ or, in terms of mass ratios $`q=M_\mathrm{p}/M_{}`$, we have models with $`q=1\times 10^4,2\times 10^4,5\times 10^4`$, and $`q=1\times 10^3`$. It has been recently demonstrated (Paper II) that, concerning gravitational torques and mass accretion, 3D models agree quite well with 2D ones in that mass range.

However, since the disk thickness is sufficiently small, even for low protoplanetary masses (see Fig. 4), we go as far as to simulate an additional C-model with $`M_\mathrm{p}=20M_{}`$. Figure 5 shows the disk semi-thickness in the planet’s neighborhood for three selected masses. Indicating with $`R_\mathrm{H}=r_\mathrm{p}\left(q/3\right)^{1/3}`$ the planet’s Hill radius, one can clearly see that the condition $`HR_\mathrm{H}`$ locally holds even in the smallest mass case. The figure also shows that protoplanets dwell in cavities and thus radiation from the central star cannot directly reach the surrounding matter.

The smoothing length $`\delta `$ is set to $`2\times 10^2R_\mathrm{H}`$ (see eq. ). All of the models are executed in an accreting and a non-accreting mode. When accretion onto the planet is allowed, the procedure outlined in Paper I is utilized. The accretion region has a radius $`\kappa _{\mathrm{ac}}=0.1R_\mathrm{H}`$, which should be short enough to let the whole accretion procedure be almost independent of the removal characteristic time, as proved by Tanigawa & Watanabe (2002).

Figure 6 illustrates the influence of the extent of the domain’s radial extent for Jupiter C-models. Simulations provide quite consistent outcomes over the matching domain. Differences exist only close to the inner border, due to inflow boundary conditions.

### 4.1 Initial and Boundary Conditions

Initial conditions require the assignment of the three functions $`\mathrm{\Sigma }\left(t=0\right)`$, $`T\left(t=0\right)`$, and $`𝒖\left(t=0\right)`$ that we suppose to be axi-symmetric because we start from a stationary disk.

The theoretical surface density profile of a stationary accretion disk, with constant $`\nu `$, is a power-law of the radial distance $`r`$ (Lynden-Bell & Pringle, 1974), as also found numerically in § 3.1. Thus, as starting density distribution for our simulations, we set

$$\mathrm{\Sigma }\left(t=0\right)=\mathrm{\Sigma }_0\sqrt{\frac{r_\mathrm{p}}{r}}.$$
(33)
where $`\mathrm{\Sigma }_0`$ is determined by the disk mass value $`M_\mathrm{D}`$.

Calculations carried out for different viscosity magnitudes imply different temperature regimes (see Fig. 4). According to them, the initial temperature distribution is prescribed. The behavior of the profiles in the top-right panel of Figure 4 can be roughly reproduced by the power-law

$$T\left(t=0\right)=T_0\left(\nu \right)\left(\frac{r_\mathrm{p}}{r}\right)^{1.8}.$$
(34)
The initial temperature at $`r=r_\mathrm{p}`$ is given in Table 1 for three values of $`\nu `$. Finally, the initial circumstellar flow $`𝒖\left(t=0\right)`$ is a Keplerian one, corrected for the grid rotation.

As for the boundary conditions, all models are run with a partially open inner radial border and a reflective outer one. Thus, matter is free to flow out of the computational domain, at $`r=r_{\mathrm{min}}`$, but the opposite is not allowed. This is the same expedient invoked in Paper I and Paper II to artificially mimic the mass accretion toward the central star and avoid spurious wave reflection at the inner domain edge, which is the closer to the planet. The flow field is Keplerian both at $`r_{\mathrm{min}}`$ and $`r_{\mathrm{max}}`$: $`𝒖=[0,r\left(\mathrm{\Omega }_\mathrm{K}\mathrm{\Omega }_\mathrm{p}\right)]`$.

### 4.2 Numerical Specifications

The basic computational mesh is made of $`N_r\times N_\phi =143\times 423`$ cells, while sub-grid patches are $`64\times 64`$. All of the computations are based on a five-grid hierarchy. Thus, the finest resolution ranges from $`1.3\times 10^2R_\mathrm{H}`$ (if $`M_\mathrm{p}=1M_\mathrm{J}`$) to $`2.9\times 10^2R_\mathrm{H}`$ (if $`M_\mathrm{p}=33M_{}`$). Further details on numerical issues have been presented in Paper I and Paper II.

Despite we have not implemented a sophisticated energy transport treatment, these simulations take between 35 and 40% longer than those discussed in Paper I, for equal size grid hierarchies.

## 5 Global Model Properties

We now examine the behavior of global quantities for some selected models because the large scale look of disks with planets is, observationally, of extreme importance. Figure 7 shows both the averaged surface density and temperature in Jupiter-mass computations with different kinematic viscosities. We note that matter (left panel) is more uniformly distributed in the H- (solid line) and W-model (short-dash line), with respect to the C-model (long-dash line). In these cases, the density gap is not very deep and actually looks like a depression, because of the larger viscosity coefficient and the depletion of the inner disk. On the one hand, in fact, a high viscosity prevents the formation of wide and deep gaps. On the other hand, a large viscosity facilitates the inner-disk depletion, since $`\dot{M}_\mathrm{D}\nu \mathrm{\Sigma }`$ (Lynden-Bell & Pringle, 1974). The disk depletion rate that we measure in these simulations ranges from few times $`10^8`$ to $`10^7M_{}\mathrm{yr}^1`$.

In contrast, in the low viscosity regime (C-model) a well-defined and deep gap is carved in, where the density is more than an order of magnitude lower than it is in the other two models. The density peak visible in the middle of the gap is the signature of the circumplanetary disk. Since material removed from gap regions is disposed at short distances, due to the small damping length of the waves launched by the planet, gap shoulders have high densities. It has been argued that these density enhancements might trigger planet formation. If so, low viscosity disks should be more efficient in doing that.

Following the reduction of thermal energy in the low density regions, a temperature gap accompanies the density gap (right panel, Fig. 7). In the C-model, temperatures drop below $`20\mathrm{K}`$ and it generally stays below $`40\mathrm{K}`$. As expectable, the temperature gap has a counterpart in the pressure scale height, as seen in Figure 8. Such a deep trough represents an additional evidence that, as long as an inner disk exists, protoplanets are shaded against stellar radiation.

Global views of both density and temperature 2D-distributions are displayed in Figure 9. Top panels illustrate the surface density, while the bottom ones show the temperature structure. From the density maps one can see that disk spirals get closer to each other toward the outer border of the domain. In fact, circumstellar spirals propagate at a velocity equal to the sound speed, which is proportional to $`\sqrt{T}`$ (see eq. ). Thus, the propagation velocity reduces as $`r`$ increases and therefore faster waves overtake slower ones. Past the gap region, temperature decays more rapidly in the C-model (Fig. 7, right panel), hence this phenomenon is more evident. From the bottom panels one can also deduce that spiral perturbations in the thermal energy and mass densities are quite phase coherent and thus disturbances in the thermal structure are rather weak.

Shifting toward smaller masses, the gap is gradually refilled. This is visible in Figure 10, where C-models of different mass are compared. At $`M_\mathrm{p}=0.1M_\mathrm{J}`$ (i.e., $`M_\mathrm{p}=33M_{}`$), the gap is one order of magnitude as shallow as it is when $`M_\mathrm{p}=0.5M_\mathrm{J}`$. We will further address the issue of the gap structure in § 5.2, since gaps are considered promising features for detecting embedded planets. As regards the circumstellar disk spirals, only feeble traces are left, when $`M_\mathrm{p}=0.1M_\mathrm{J}`$, both in the temperature and density distributions. The absence of global disk features produced by an embedded body invalidates direct imaging observations as a detection tool in the intermediate-mass range, i.e., around a few times of Uranus’ mass.

Concerning the emitted energy, we can notice that the optical thickness of the models is generally larger than one, throughout the computational domain (Fig. 11, left panel). The averaged $`\tau `$ stays above $`1`$ also at the outer disk edge. This is in agreement with what was found for accretion disks around T Tauri stars (D’Alessio et al., 1998). The density and temperature gap reflects in a similar feature in the distribution of $`\tau `$, as shown in Figure 11 (left panel). Nonetheless, even within the gap region, material is, on the average, optically thick ($`\tau 1`$). The exception is represented by the C-model, where gap optical thickness reaches values around $`0.01`$ at $`\phi =\phi _\mathrm{p}+\pi `$ (see Fig. 11, right panel). The most diluted portions of the gap, situated where the circumplanetary disk merges into the gap, have values of $`\tau `$ on the order of $`10^2`$$`10^3`$.

Finally, we examine the emitted flux $`\mathrm{\Lambda }`$, computed according to equation (14) (see Fig. 12). From the observational point of view, disk gaps could represent a probe for protoplanet detection. In fact, this is by far the most extensive imprint that a planet leaves on a circumstellar disks. Prospective studies on observability of gaps due to Jupiter-mass protoplanets, have already been presented by Wolf et al. (2002); Steinacker & Henning (2003), and speculations on a developing gap around the T Tauri star TW Hya have been reported by Calvet et al. (2002). Here we do not intend to address the issue of whether gaps are really observable or to what extent they are. However, we have to notice that a necessary but not sufficient condition is that they must be wide, deep, and with sharp edges in order for the flux emitted from this region to have the strongest contrast with respect to the surrounding environment (Rice et al., 2003). From the bottom panels of Figure 12, it appears clear that low kinematic viscosities would favor this kind of investigation. Disk spirals are probably too elusive to be captured by present-day ground-based instruments, as shown by the flux maps. We also argue that it would seem rather unlikely to detect planetary masses smaller than Jupiter’s, by means of these measurements. Right panels of Figure 12 display the emitted flux in case of a $`M_\mathrm{p}=0.1M_\mathrm{J}`$ model. The quantity $`\mathrm{\Lambda }`$ furnished by high-viscosity model (top) looks quite smooth and only some inhomogeneities appear in C-model (bottom).

### 5.1 Models with $`\alpha `$-viscosity

So far we have considered disks with a constant $`\nu `$. However, we have also investigated the effects of a different viscosity prescription, namely an $`\alpha `$-type viscosity (Shakura & Sunyaev, 1973). Two Jupiter-mass accreting models were run in this mode, one with $`\alpha =10^2`$ (model A2) and the other with $`10^3`$ (model A3). Such values are representative of those recently estimated by Papaloizou & Nelson (2003) and Winters, Balbus, & Hawley (2003) who performed MHD disk simulations, with embedded Jupiter-size bodies, and could therefore evaluate the stress parameter $`\alpha `$ in a self-consistent fashion. Indeed, referring to the constant viscosity models discussed above, by making the azimuthal mean $`2\pi /\alpha =\left(1/\nu \right)c_sHd\phi `$, we get a value of $`\alpha `$ ranging from $`3\times 10^2`$ in the inner disk to $`4\times 10^3`$ in the outer disk. Specifically, C-models yield an overall trend of $`\alpha `$ which is the closest to $`10^2`$. And in fact outcomes from $`\alpha `$-viscosity model A2 are not much different from those obtained from the Jupiter-mass C-model. The quantity $`\mathrm{\Sigma }`$ is very similar to that illustrated in left panel of Figure 7 (long-dashed line) whereas $`T`$ is only slightly different (near the inner border) from the one in the right panel of the same Figure.

Some more important deviations from what has been seen thus far are instead encountered in the A3-model, as displayed in Figure 13. The averaged surface density shows a significant hump (top-left panel, solid line), inside the gap region, not observed in the other viscosity regimes. This is not caused by the mass accumulation around the planet but it is produced by material lingering around Lagrangian points L4 and L5, as it appears from the ring-like traces in the middle of the density gap in the top-right panel Figure 13. This also appear to be a persistent feature, in fact no change is measured over the last 50 orbital periods. L4 and L5 points are equilibrium locations for the restricted three-body problem. Hence, one can expect that the thermal energy is larger near such points than anywhere else in the gap. In turn, viscous torques, which are enhanced by the locally increased kinematic viscosity are more efficient in contrasting gravitational torques. Thereby, material is more likely forced to evolve on the much longer viscous time scale $`t_{\mathrm{vis}}r_\mathrm{p}^2/\nu =r_\mathrm{p}/\left(\alpha c_sh\right)`$.

As for the temperature (Fig. 13, bottom panels), the largest discrepancies between models A2 and A3 occur toward the inner border.

### 5.2 Gaps as Planet Signatures

We have seen in the previous two sections that the gap structure drastically depends on the viscosity regime. The standard criterion to open a gap in a disk requires two conditions (Lin & Papaloizou, 1985, 1993). The first is a thermal condition given by:

$$\left(\frac{M_\mathrm{p}}{M_{}}\right)>3\left(\frac{H}{r}\right)^3=3h^3,$$
(35)
i.e., the Roche lobe must be larger than the disk scale height. Because of planet’s gravity (see eq. ) and the low temperatures established in gap regions, this condition is always fulfilled when $`M_\mathrm{p}0.1M_\mathrm{J}`$ (see Fig. 8). The second condition concerns the requirement that tidal torques exceed viscous torques and reads

$$\left(\frac{M_\mathrm{p}}{M_{}}\right)>\left(\frac{40}{}\right),$$
(36)
where $`=r^2\mathrm{\Omega }_\mathrm{K}/\nu `$ is the Reynolds number. The right hand side of equation (36) is proportional to $`\nu `$ and in our case goes from $`4\times 10^3`$ (H-models) to $`4\times 10^4`$ (C-models). Thereby, we should not obtain gaps in any of the H-models. Figure 7 (solid line) decently agrees with this prediction, since only a trough is dug in H-model. However, condition (36) is partly violated in the low-viscosity regime (Fig. 10), because a clear gap is visible for the $`M_\mathrm{p}=0.2M_\mathrm{J}`$ case and a trough appears for $`M_\mathrm{p}=0.1M_\mathrm{J}`$. In Table 2 we report the gap occurrence and depth for some of the models.

One reason why the above criterion could not properly apply to our calculations lies in the fact that it was derived for simple polytropic disks, while we simulate also the thermal evolution of the system. As implied by the profiles in the top-left panel of Figure 13, it is even harder to predict the gap presence when the kinematic viscosity depends on the local fluid temperature.

## 6 Circumplanetary Disks

The first distinction about the circumplanetary flow has to be made according to whether the planet is accreting or not. Therefore, we shall carry out separate discussions, based on H- and C-models, as well as $`\alpha `$-viscosity models A2 and A3.

### 6.1 Accreting Models

Figure 14 indicates that, inside of the Roche lobe, the general aspect of the flow circulation around $`1M_\mathrm{J}`$ and $`0.1M_\mathrm{J}`$ protoplanets resembles that obtained with local isothermal models in two dimensions (Paper I). However, specific characteristics of the flow do differ. A circumplanetary disk, extending over a fair fraction of Roche lobe, can be identified for both masses and viscosity regimes. The mass of such disks around one Jupiter-mass protoplanets, within the 80% of the Hill radius, is $`1.03\times 10^5M_\mathrm{J}`$ in the H-model and $`5.36\times 10^6M_\mathrm{J}`$ in the C-model. As regards calculations with $`M_\mathrm{p}=0.1M_\mathrm{J}`$ planets, the mass measured in the circumplanetary disk (or sub-disk, for brevity) of the C-model is $`4.03\times 10^6M_\mathrm{J}`$, while the amount is nearly two times less in the H-model. While the flow dynamics in the Jupiter-mass model with $`\alpha =10^2`$ (left panel in Fig. 15) is similar that of the C-model, the model with $`\alpha =10^3`$ presents some peculiarities. As shown in right panel of Figure 15, the size of the circumplanetary disk is smaller than the Roche lobe. Moreover, the streams of matter associated with the spiral waves are narrower. The sub-disk mass measured in these cases is $`3.9\times 10^6M_\mathrm{J}`$ (A2) and $`1.5\times 10^6M_\mathrm{J}`$ (A3).

The azimuthal average of the surface density distribution around the planet, inside the Roche lobe, can be fitted by a relation of the type (see Fig. 16):

$$\mathrm{\Sigma }\mathrm{\Sigma }_0\mathrm{exp}\left[a_1\left(\frac{s}{R_\mathrm{H}}\right)\right],$$
(37)
where $`s/R_\mathrm{H}1`$ and larger than a threshold length (either $`0.1`$ or $`0.2R_\mathrm{H}`$). The values of $`\mathrm{\Sigma }_0`$ and $`a_1`$ are reported in Table 3.

The presence of spiral features in the density distribution is an indication that the circumplanetary flow is Keplerian-like. Indeed, decomposing the velocity field $`𝒖`$ into the in-fall velocity (toward the planet) $`w_{\mathrm{in}}`$ and the rotational velocity (counter-clockwise around the planet) $`w_{\mathrm{rot}}`$, it turns out that the former is more than an order of magnitude less than the latter. In Jupiter-mass models, the rotational velocity within $`0.5R_\mathrm{H}`$ of the planet drops approximately as $`s^{0.6}`$, nearly independently of the viscosity regime investigated in these computations. For $`0.1M_\mathrm{J}`$ models, this ratio declines at a slightly higher rate. Consequently, for $`s[0.1,0.5]R_\mathrm{H}`$, the ratio of $`w_{\mathrm{rot}}`$ to the local Keplerian velocity is always larger than $`0.7`$.

Comparing the right panels (C-models) of Figure 14 to Figure 13 (top and middle right panels) in Paper I, one can clearly see that spiral perturbations are less intense, as we show below, and more open. The latter circumstance is related to the lower values of the Mach number in the circumplanetary flow, which governs the inclination of the spiral wave with respect to the direction of the rotational motion. In fact, in local isothermal models $`_{\mathrm{iso}}\sqrt{q\left(r_\mathrm{p}/s\right)}/h`$. If $`M_\mathrm{p}=1M_\mathrm{J}`$, $`_{\mathrm{iso}}`$ drops from $`8`$, at $`s=0.1R_\mathrm{H}`$, to $`2`$, at $`s=1R_\mathrm{H}`$. As comparison, in the circumplanetary disk displayed in the upper-left panel of Figure 14, $``$ lies between $`2.4`$ and $`1.4`$, whereas the low viscosity Jupiter-model (upper-right panel) provides values between $`3`$ and $`1.5`$. Models A2 and A3 yield similar values. The main reason for this difference resides in the larger value of the sound speed. Hence, perturbations can travel faster and are less distorted by the background motion of the flow. In the Jupiter-mass cases, illustrated in Figure 14, the azimuthally averaged $``$ is approximately proportional to either $`s^{0.1}`$ (H-model) or to $`s^{0.2}`$ (C-model). In the limit of a nearly constant Mach number, wave perturbations assume the form of Archimedes’ spirals.

As pointed out in Paper II, in a three-dimensional space circumplanetary spirals are weakened (see also Bate et al., 2003). In fact, because of the flow circulation above the disk midplane (see Paper II, Figure 3, middle and bottom panels), less angular momentum is carried inside the Roche lobe by the midplane flow. As anticipated above, the spiral density features look smoother in the models presented here. This is quantitatively demonstrated in Figure 17. In the Figure the average specific angular momentum $`sw_{\mathrm{rot}}`$ obtained from the $`M_\mathrm{p}=1M_\mathrm{J}`$ C-model (solid line) is compared to that of a 2D Jupiter-mass local-isothermal model (short-dash line) and to that of a 3D Jupiter-mass model from Paper II (long-dash line), which also assumes a quasi-isothermal gas around the planet. In the 3D case, the quantity $`sw_{\mathrm{rot}}`$ is computed at $`z=0`$, i.e., in the disk midplane. The latter two models are quasi-isothermal in the sense that the temperature distribution is fixed and depends only on the distance distance from the star $`r`$ according to the relation $`T\left(\mu /\gamma \right)h^2M_{}/r`$. Thus, around the planet’s location $`r=r_\mathrm{p}`$, $`T90\mathrm{K}`$ if the disk aspect ratio $`h`$ is $`0.05`$, $`\mu =2.39`$, and $`\gamma =1.4`$. From Figure 17 one can see that in the outer part of the Roche lobe ($`s0.3R_\mathrm{H}`$), where spiral perturbations are strongest, the discrepancy between the two quasi-isothermal models lies between 20 and 30%, whereas the two-dimensional thermal model provides a specific angular momentum only 10% larger than that of the three-dimensional quasi-isothermal model. This occurrence can be attributed to the better pressure structure achieved by the calculations that account for the thermal energy budget.

Circumplanetary disks should not suffer from gravitational instabilities since the Toomre parameter $`Q`$ is orders of magnitude larger than unity, as proved by Figure 18. This arises from the small value of the ratio between the mass of the circumplanetary disk and that of the planet.

Two-dimensional temperature distributions are rather symmetric with respect to the planet’s position, though they are marked by weak spiral perturbations. Figure 19 (top panels) shows the azimuthal averaged temperature around protoplanet with $`M_\mathrm{p}=1M_\mathrm{J}`$ (left) and $`M_\mathrm{p}=0.1M_\mathrm{J}`$ (right). In the Jupiter-mass case, maximum temperatures range from roughly $`1500`$ to $`1000\mathrm{K}`$, whereas $`T`$ reaches the value of the ambient medium ($`50\mathrm{K}`$) toward the border of the Hill sphere. Comparing H- and C-models, one can realize that the average temperature in the inner part of the Roche lobe differ by less than a factor two, whereas the $`\alpha `$-viscosity model A3 shows significant lower temperatures. In the low-mass case ($`M_\mathrm{p}=0.1M_\mathrm{J}`$), the peak temperature is around $`1000\mathrm{K}`$ and between $`80`$ and $`40\mathrm{K}`$ at $`s=R_\mathrm{H}`$, respectively for both H- and C-model. The average temperature differs by less than $`30\mathrm{K}`$ (Fig. 19, upper-right panel).

Over the entire sub-disk domain, $`T`$ can be well fitted by the curve:

$$TT_{\overline{s}}\left(\frac{\overline{s}}{s}\right)^\xi ,$$
(38)
in which the length $`\overline{s}`$ is set to $`0.1R_\mathrm{H}`$, for convenience. Fitting parameters are reported in Table 4, together with the validity range of the fitting function. From entries in the Table, it turns out that the temperature generally falls off as $`1/s`$.

Bottom panels of Figure 19 display the average effective temperature $`T_{\mathrm{eff}}`$ (see eq. ) in the vicinity of protoplanets. Figure 20 shows that the magnitude of the optical thickness $`\tau `$ is between $`10`$ and $`100`$, hence equation (14) yields $`T_{\mathrm{eff}}T/\tau ^{0.25}`$ and thus there is, at most, a factor three between the two temperatures. It should be noted that only the model A3 presents optically thin regions toward the edges of the sub-disk.

Accreted matter also contributes to the energy budget of the circumplanetary disk via dissipation of gravitational energy into heat. This source of energy is not treated explicitly in these simulations. However, it is considered implicitly as a compression work term arising from the removal of matter in the accreting region. Tanigawa & Watanabe (2002) assumed the dissipation of gravitational energy as the only heating source and derived temperature profiles of sub-disks with different sound-speed regimes, performing local simulations with an isothermal equation of state. They found that the temperature would scale as $`1/\sqrt{s}`$, independently of the planetary mass.

We have to point out that, since radial transport of radiation is not treated here, one should expect large temperature gradients to be smoothed out. But we have seen that temperature profiles are not extremely steep, in the range of distances covered by these models, hence correction by radial transfer should not be dramatically relevant, as we are going to demonstrate.

#### 6.1.1 Radial Radiation Transfer in Circumplanetary Disks

Since the temperature gradient in the vicinity of a protoplanet may be quite steep, one of the hypothesis used in § 2.1 might locally break down. We refer to the assumption that the radiative flux in the vertical direction overwhelms those in the horizontal direction. With a semi-analytical analysis it is possible to make an a posteriori check of that hypothesis and shed some light over this matter. For simplicity we will suppose that the temperature distribution has a cylindrical symmetry around the planet, which is indeed what is observed in all of the models. Hence the total flux of energy transferred via radiation is

$$\mathrm{\Lambda }_{\mathrm{}}^+\mathrm{}\frac{F_z}{z}dz+_{\mathrm{}}^+\mathrm{}\frac{1}{s}\frac{}{s}\left[sF_s\right]dz=\mathrm{\Lambda }_z+\mathrm{\Lambda }_s,$$
(39)
where $`s`$ is the distance from the planet. The two terms on the right-hand side can be written as

$$\mathrm{\Lambda }_z\frac{128\sigma _\mathrm{R}}{3\kappa \rho }\frac{T^4}{H},$$
(40)
and

$$\mathrm{\Lambda }_s\xi ^2\frac{128\sigma _\mathrm{R}}{3\kappa \rho }H\frac{T^4}{s^2}.$$
(41)
Equation (40) was obtained by setting $`T/zT/H`$ in equation (39), whereas equation (41) was derived by adopting equation (38), which we have checked to hold as long as $`\overline{s}sR_\mathrm{H}`$. The validity of the relation (9), namely that $`\mathrm{\Lambda }\mathrm{\Lambda }_z`$, depends on the ratio of the left-hand sides of equations (40) and (41):

$$\frac{\left|\mathrm{\Lambda }_s\right|}{\left|\mathrm{\Lambda }_z\right|}=\xi ^2\left(\frac{H}{s}\right)^2.$$
(42)
According to equation (42), only when $`H<s`$ the radiative cooling (eq. ) included in the energy equation is also locally a good approximation. As shown in Figure 21, this seems indeed to be the case in models involving Jupiter-size planets. To cast the above ratio in a more explicit form, one can use the expression for $`H`$ (eq. ) in the limit $`sr`$, which yields

$$\left(\frac{H}{s}\right)^2=\frac{c_s^2}{\gamma }\left(\frac{s}{GM_\mathrm{p}}\right).$$
(43)
By using the form of the sound speed in equation (5), one gets

$$\left(\frac{H}{s}\right)^2=\left(\frac{kT}{\mu m_\mathrm{H}}\right)\left(\frac{s}{GM_\mathrm{p}}\right).$$
(44)
Therefore, because of the assumption made above on the temperature profile, equation (42) becomes

$$\frac{\left|\mathrm{\Lambda }_s\right|}{\left|\mathrm{\Lambda }_z\right|}=\xi ^2\left(\frac{kT_{\overline{s}}}{\mu m_\mathrm{H}}\right)\left(\frac{\overline{s}}{GM_\mathrm{p}}\right)\left(\frac{\overline{s}}{s}\right)^{\xi 1}.$$
(45)
Since the distance $`s`$ is constrained by the applicability range of equation (38), the ratio given by equation (45) is supposedly meaningful roughly for $`s[0.1R_\mathrm{H},R_\mathrm{H}]`$, for the models considered here. With the appropriate values, the ratio (45) becomes

$$\frac{\left|\mathrm{\Lambda }_s\right|}{\left|\mathrm{\Lambda }_z\right|}=0.145\xi ^2\left(\frac{T_{\overline{s}}}{1000\mathrm{K}}\right)\left(\frac{M_\mathrm{J}}{M_\mathrm{p}}\right)^{2/3}\left(\frac{\overline{s}}{s}\right)^{\xi 1}.$$
(46)
Equation (46) is an evaluation of equation (45) at a distance $`s`$ from the planet, whose length shortens as $`M_\mathrm{p}^{1/3}`$. Since the squared sub-disk aspect ratio grows as $`1/M_\mathrm{p}`$, the flux ratio goes as $`M_\mathrm{p}^{2/3}`$. From Table 4, we see that $`\left|\xi 1\right|`$ is between $`0.03`$ and $`0.12`$ in C-models and between $`0.04`$ and $`0.21`$ in H-models. Therefore, the right-most term in equation (46) can be considered as a unity term. Thus, around $`s\overline{s}=0.1R_\mathrm{H}`$, $`\left|\mathrm{\Lambda }_s\right|/\left|\mathrm{\Lambda }_z\right|=0.13`$ for the Jupiter-mass H-model and a factor two as small for the C-model. Similar small numbers are obtained from both the A2 and A3 model. For $`M_\mathrm{p}/M_\mathrm{J}=0.1`$, the ratio increases, respectively, to $`0.2`$ and $`0.3`$, which is consistent with the applicability limit of the two-dimensional approximation. Thereby, apart from regions closer than $`0.1R_\mathrm{H}`$ to the accreting planet which, probably, cannot be properly investigated by these computations, radiation transport in the radial direction does not play a major role in the energy budget of circumplanetary disk material. This is in agreement with the standard scenario of the late stages of protojovian disks (e.g., Coradini et al., 1989; Canup & Ward, 2002).

### 6.2 Non-accreting Models

Accreting models assume that the planetary accretion rate can adjust to the rate at which matter can be supplied by the surrounding environment. We now investigate the opposite situation, namely when no further accretion onto the protoplanet is allowed.

The overall flow structure and density distribution of non-accreting protoplanets are drastically different from those of accreting ones. This can be clearly seen in Figure 22, where we show the surface density and the flow field for $`M_\mathrm{p}=1M_\mathrm{J}`$ (upper panels) and $`0.1M_\mathrm{J}`$ (lower panels) objects. The standard scenario of a Keplerian-like sub-disk is completely altered. This is a consequence of the large pressure gradient, built up by the large density gradient, which allows matter to move on orbits not constrained by the centrifugal balance. Therefore, the two-arm spiral feature is replaced by a more complex system of multiple shock fronts across which material is first deflected toward the planet and then ejected outward, as also indicated by the sign of the in-fall velocity $`w_{\mathrm{in}}`$ within the region.

The flow strongly diverges in a zone that has a radius equal to $`0.4R_\mathrm{H}`$ in the upper-left panel (H-Jupiter model) and to $`0.2R_\mathrm{H}`$ in the upper-right panel (C-Jupiter model). The average surface density inside of $`R_\mathrm{H}/3`$ is between one and two and a half orders of magnitude larger than that measured in accreting models. The mass collected in $`s0.8R_\mathrm{H}`$ is $`1.97\times 10^4M_\mathrm{J}`$ and $`4.30\times 10^5M_\mathrm{J}`$ for the Jupiter-mass H-model and the C-model, respectively. The azimuthal average of the circumplanetary density, around Jupiter-size objects, can be approximated to $`\mathrm{\Sigma }1.2\times 10^4\left(0.1R_\mathrm{H}/s\right)^{2.8}\mathrm{g}\mathrm{cm}^2`$ (H-model) and to $`\mathrm{\Sigma }1.5\times 10^3\left(0.1R_\mathrm{H}/s\right)^{2.7}\mathrm{g}\mathrm{cm}^2`$ (C-model). The top panel of Figure 23 shows the averaged surface density profiles around $`M_\mathrm{p}=1M_\mathrm{J}`$ (left).

A major difference can be also observed in the sub-disk circulation of low-mass models (Fig. 22, lower panels). In fact, the fluid rotates in a clockwise direction. In general, the direction of rotation in sub-disks is determined by the balance among the Coriolis force, the pressure gradient, and the gravitational attraction by the planet. Basically, referring to the lower-right panel of Figure 22, the Coriolis force deflects fluid elements, orbiting at $`rr_\mathrm{p}`$, rightward. Therefore, it forces matter to cross the gap and to reconnect to the other side, while still moving in a position upstream of the perturber body. In contrast, the term $`\mathrm{d}P/\mathrm{d}s`$, supposedly positive, opposes to reconnection when fluid is upstream of the planet’s location but favors it when matter is downstream of the planet. The prograde rotation that has been encountered so far indicates that the Coriolis deflection overwhelms the pressure gradient. For this reason reconnection across the gap always occurs upstream of the planet. Evidently, in the low-mass models ($`M_\mathrm{p}=33M_{}`$) displayed in the bottom panels of Figure 22, the pressure gradient drives material, flowing from upstream, past the planet’s position and forces it to circulate clockwise around the perturber. Consequently, matter entering the Hill sphere has an angular momentum anti-parallel (in the planet’s frame) to that of the circumstellar disk. In the range from $`0.1`$ to $`1R_\mathrm{H}`$, the magnitude of the rotational velocity $`\left|w_{\mathrm{rot}}\right|`$ increases roughly as $`\sqrt{s}`$. The behavior of $`\mathrm{\Sigma }`$ resembles that of a power-law, as in Jupiter-mass cases, with powers $`0.9`$ and $`0.5`$, respectively for the H- and C-model (Fig. 23, right panel). Interestingly, no large differences are seen in the two density profiles.

The high pressure gradient is also caused by the large temperature gradient. Parameters for analytic approximations of the average temperature $`T`$, according to equation (38), are reported in Table 5.

Maximum temperatures, at $`s=0.1R_\mathrm{H}`$, reach $`1900\mathrm{K}`$ and $`1200\mathrm{K}`$, in $`M_\mathrm{p}=1M_\mathrm{J}`$ H- and C-models respectively (bottom-left panel in Fig. 23). At the limit of the Roche lobe the temperature is somewhat higher than that measured in accreting models and depends on the gap structure. They vary from $`100`$ to $`200\mathrm{K}`$. Low mass models furnish temperatures that are very similar in both viscosity regimes.

It is worthy to notice that in either planetary-mass and viscosity case, the resulting temperature profiles are steeper than the ones shown in Figure 19 only when $`s<0.1R_\mathrm{H}`$, but they are nearly as steep past this distance. Actually, from Table 5 it appears the toward low planetary masses temperature profiles are even flatter. Performing the same check as done in the previous section, it turns out that $`\left|\mathrm{\Lambda }_s\right|/\left|\mathrm{\Lambda }_z\right|`$ is at most $`0.3`$. In fact, around Jupiter-mass planets, the ratio $`H/s`$ varies between $`0.3`$ and $`0.5`$, for distances larger than $`0.1R_\mathrm{H}`$. Hence, even in these circumstances the type of energy equation implemented here yields a reasonable description of the thermal structure of sub-disks, down to $`0.1R_\mathrm{H}`$.

## 7 Accretion and Migration

The evaluation of the mass accretion rate of protoplanets is carried out by reducing the surface density around the planet, according to a relation of the type $`\mathrm{\Delta }\mathrm{\Sigma }/\mathrm{\Sigma }=\mathrm{\Delta }t/\tau _{\mathrm{ev}}`$ (see, Paper I for details). This reduction process has a time scale $`\tau _{\mathrm{ev}}`$ that is much smaller than the integration time step $`\mathrm{\Delta }t`$. In Figure 24 the planetary accretion rate $`\dot{M}_\mathrm{p}`$ evaluated in these computations is compared to that achieved with local-isothermal 2D (Paper I) and 3D (Paper II) models, where the disk mass was rescaled to match the values adopted here (see § 4).

In the context of the accretion history of protoplanets, it is worthwhile to stress that 2D calculations are only able to simulate the late stages of the accretion process, i.e., when it mainly occurs via mass transfer in the circumplanetary disk. This is consistent with the parameterization adopted here to describe the circumstellar environment.

From Figure 24, one can see that C-models (filled squares) provide accretion rates quite similar to those obtained from models with a fixed temperature distribution (crosses). This occurrence is related to the quasi-Keplerian, circumplanetary flow observed in both kinds of computations. Jupiter-mass models with an $`\alpha `$-type viscosity (open triangles) give rates which are similar to that of the C-model when $`\alpha =10^2`$ but somewhat lower ($`\dot{M}_\mathrm{p}=6\times 10^3M_{}\mathrm{yr}^1`$) when $`\alpha =10^3`$, consistently with the low sub-disk mass measured in this case. H-models provide values of $`\dot{M}_\mathrm{p}`$ (filled diamonds) that, in case of Jupiter-size planets, fall above the estimates from all of the other models. According to the accretion procedure, the larger the amount of matter contained in the accretion region $`s=\kappa _{\mathrm{ac}}`$, the larger the accreted mass is. Within $`s=0.1R_\mathrm{H}`$ of a Jupiter-mass planet, the mean density computed in the H-model is roughly three times larger than that in the C-model (Fig. 14). This is related to the background gap density that is much higher in H-models, as reported in Table 2 (see also Fig. 7), due to the larger fluid viscosity.

Figure 25 shows the total gravitational torque $`𝒯_\mathrm{D}`$ exerted by disk material onto accreting planets, as the inverse of the migration time scale ($`𝒯_\mathrm{D}\left|\dot{r}_\mathrm{p}\right|/r_\mathrm{p}1/\tau _\mathrm{M}`$). This quantity is computed from the gravitational force acting on the planet, as explained in Paper I. The length $`\beta `$ represents the radius of the circular region around the planet’s position whose contribution is not taken into account. In Figure 25, this length is set equal to $`R_\mathrm{H}`$. Outcomes from C-models, and generally those from H-models, compare well to the local-isothermal computations. Models with $`\alpha `$-viscosity provide migration rates on the same magnitude as well. C-models show a rather constant value of $`\tau _\mathrm{M}`$ as a function of the planet’s mass, whereas H-models present a minimum at $`M_\mathrm{p}=0.2M_\mathrm{J}`$. Migration time scales of H-models differ by less than a factor two from those of C-models. Only for $`M_\mathrm{p}=33M_{}`$, $`\tau _\mathrm{M}`$ is appreciably longer. One may attribute this effect to the absence of a real gap in this model (see Table 2), and thus to an early onset of the Type I drifting regime.

In Table 6, we separate the contributions of different disk regions to the total torque $`𝒯_\mathrm{D}`$, introducing the inner disk ($`𝒜`$), the outer disk ($`𝒟`$) and the gap region ($``$ and $`𝒞`$). Further details on the regions are reported in the Table’s note. All contributions are normalized to $`\left|𝒯_\mathrm{D}\right|`$. The width of the gap region is fixed to $`5R_\mathrm{H}`$ (see Fig. 7) and the zone behind the planet ($``$) is distinguished from the one ahead of it ($`𝒞`$). Torques exerted from inside the Hill sphere are not accounted for. A large fraction of the contributions arising from zones $``$ and $`𝒞`$ is actually generated in the coorbital regions, within a few Hill radii from the planet. From the entries in Table 6 one can argue that coorbital torques in H-models are larger in magnitude than they are in C-models. This is a consequence of the shallower gap of H-models. The situation is reversed in the $`\alpha `$-viscosity models because of the density humps at L4 and L5 Lagrangian points.

As for the influence of the parameter $`\beta `$ on the total torque, it turns out that H-models are generally more sensitive than C-models. By examining the function $`𝒯_\mathrm{D}=𝒯_\mathrm{D}\left(\beta \right)`$, some more insight can be gained. For $`M_\mathrm{p}=1M_\mathrm{J}`$ models (Fig. 26, top-left panel), no significant net torque arises from within the Hill sphere in the C-model (dashed line). The situation appears to be different for the H-model (solid line) because $`𝒯_\mathrm{D}`$ decreases to a minimum, around $`\beta =R_\mathrm{H}/2`$, and then increases and changes sign around $`\beta =0.1R_\mathrm{H}`$. Jupiter-mass $`\alpha `$-viscosity models (Fig. 26, top-right panel) display a smooth behavior since the total torque is rather constant down to $`\beta =R_\mathrm{H}/2`$. When $`\beta `$ diminishes, $`𝒯_\mathrm{D}`$ increases (i.e., it has smaller absolute values) and the model with $`\alpha =10^3`$ (dashed line) shows a change of sign at $`\beta =0.4R_\mathrm{H}`$. When $`M_\mathrm{p}=0.1M_\mathrm{J}`$ (Fig. 26, bottom-left panel), in either case, $`\mathrm{d}𝒯_\mathrm{D}/\mathrm{d}\beta <0`$, because of the positive torques arising within the sub-disk. The rate of change is more drastic in the H-model, where the total torque becomes positive and keeps increasing when $`\beta R_\mathrm{H}/2`$. Torques acting on the $`20M_{}`$ planet (Fig. 26, bottom-right panel) are strongly dependent on $`\beta `$, as positive torques are exerted down to $`0.4R_\mathrm{H}`$ and large negative torques arise between $`0.2R_\mathrm{H}`$ and $`0.4R_\mathrm{H}`$.

## 8 Conclusions

In this paper we presented two-dimensional simulations of disk-planet interactions in which we considered the joint thermo-hydrodynamics evolution of the system. In order to do that, we solved an energy equation that accounts for the major processes responsible for the energetic balance of fluid parcels. We restricted the treatment of the problem to a flat geometry in order to make some simplifying assumptions on the radiative part of the energy budget, i.e., to treat the radiation transfer as a cooling term. A nested-grid technique was used to resolve different length scales, therefore both the global and local features could be investigated.

The two-dimensional approximation has been proved to work reasonably well as long as the planetary mass is larger than a tenth of the mass of Jupiter (Paper II). Because of this, only planets more massive than $`0.1M_\mathrm{p}`$ have been simulated. In order to compromise with the lingering uncertainty on the viscosity magnitude and to inquire into different temperature regimes, we considered different constant viscosity and $`\alpha `$-viscosity models. Both accreting and non-accreting planets were simulated.

From the global point of view, and roughly speaking, we can conclude that circumstellar disk models with a fixed temperature distribution provide a reasonably good description of the system, though details of the gap structure and shape of the disk spirals differ. Around Jupiter-mass protoplanets, only with kinematic viscosities as low as $`10^{15}\mathrm{cm}^2\mathrm{s}^1`$, a wide and deep gap is carved in whereas it reduces to a shallow density trough when the viscosity is ten times as large. Mean temperatures in the gap can be as low as $`20`$$`50\mathrm{K}`$, depending on the viscosity regime. The most diluted and cold gap material is located where the circumplanetary disk merges into the gap. The aspect ratio appears smaller than $`0.04`$ and it is also marked by a gap-like feature, which shades the protoplanet’s environment against stellar radiation. The models with $`\alpha =10^3`$ show a ring of matter within the gap that evolves on a viscous time scale. Regarding smaller planetary masses ($`M_\mathrm{p}0.1M_\mathrm{J}`$), both gap and spiral perturbations become quite faint features. Likely, only Jupiter-size bodies are capable of digging gaps deep enough to be observed by present-day instrumentation.

From the local point of view, we obtained the distribution of density, temperature, and other thermodynamical quantities in circumplanetary disks. Except for the $`\alpha `$-viscosity model with $`\alpha =10^3`$, rather comparable outcomes are provided by models with different viscosity regimes. The density, averaged around the planet, shows an exponential radial decline. Sub-disks have masses on the order of $`10^5`$$`10^6M_\mathrm{J}`$, which makes them gravitationally stable. Temperatures range from several hundred Kelvin degrees, at distances $`s0.1R_\mathrm{H}`$ from the planet, to values between $`50`$ and $`100\mathrm{K}`$, at the border of the Hill sphere. Azimuthal averages yield profiles dropping nearly as the inverse of the distance from the planet. Circumplanetary disks have aspect ratios $`H/s`$ around a few tenths and they tend to be optically thick. All of these properties appear to be consistent with the “fast-inflow” analytical models of protojovian disks recently developed by Canup & Ward (2002). By means of a semi-analytical treatment, we showed that, also in the sub-disk, the neglect of radial radiation transfer should not crucially affect our results. Non-accreting models furnish a rather different scenario. Around Jupiter-mass objects, circulation in the Roche lobe is not Keplerian-like because of the large pressure built up mainly by the density gradient and partly by the temperature gradient. Around $`M_\mathrm{p}=0.1M_\mathrm{J}`$, clockwise rotation establishes as a result of the weakening of the Coriolis force with respect to the pressure gradient term.

Accretion and migration rates inferred by these simulations are comparable to those evaluated with previous purely hydrodynamical local-isothermal computations. An analysis of the zonal torques indicates that coorbital torques increase as viscosity increases. Moreover, negative coorbital torques, arising from material outside of the Roche lobe and lagging behind the planet, are generally larger than those from material ahead of the planet. We also found that, in the two-dimensional approximation, the material lying inside of the Roche lobe is able to exert strong gravitational torques on the planet. Since the density gap is usually filled in H-models (see Table 2), Type I migration regime might extend to larger planetary masses (see Fig. 25).

In order to address the issue of how the third dimension affects the thermo-hydrodynamics of the system, both globally and locally around protoplanets, we are extending this study to 3D simulations in which radiation transfer is treated in the flux-limited diffusion approximation (Levermore & Pomraning, 1981).

We thank Udo Ziegler for having made available to us an early FORTRAN Version of his code Nirvana. G. D. is deeply grateful to all the friends at the INAF-OAC (Naples) and RSMAS/MPO (University of Miami) for their warm hospitality. The paper benefited from the useful comments of an anonymous referee. This work was supported by the German Science Foundation and by the UK Astrophysical Fluids Facility (UKAFF). The numerical computations reported here were carried out at the Computer Center of the University of Jena and at the UKAFF.


format: text

text:
# ROTSE-III Observations of the Early Afterglow From GRB 030329

## 1. Introduction

The ROTSE program is dedicated to recording optical observations of Gamma-Ray Bursts (GRBs) within seconds of their detection at high energies. Although the study of late-time afterglows has revolutionized our understanding of GRBs, comparatively little is known about the nature and diversity in the prompt optical and early afterglow emission (Piran, 1999). To date the ROTSE-I observations of GRB 990123 remain the only optical detections of a GRB while the burst was still bright in gamma-rays (Akerlof et al., 1999). ROTSE-I and LOTIS observations of many bursts demonstrate that prompt emission is not, in general, well represented by extrapolation of late time decay curves (Kehoe et al., 2001; Park et al., 1999). More sensitive instruments are required to study prompt emission from typical bursts.

In preparation for HETE-2 (Vanderspek et al., 1999) and the Swift satellites (Gehrels, 2000) which can distribute rapid localizations of GRBs with arcminute accuracy, the ROTSE team developed a new set of 0.45-m telescopes. The ROTSE-III instruments use fast optics (f/1.9) to yield a $`1\mathrm{°}.85\times 1\mathrm{°}.85`$ field of view over a four megapixel CCD camera. Four of these robotically controlled, automatically scheduled instruments are being erected around the globe, both to increase the amount of sky available for instant viewing and to enable around-the-clock coverage of GRB afterglow behavior (Akerlof et al., 2003). At the time of this writing, the first two instruments, labeled ROTSE-IIIa and ROTSE-IIIb, are operating at Siding Spring Observatory in New South Wales, Australia, and McDonald Observatory at Ft. Davis, Texas, USA, respectively. The installation of ROTSE-IIIc is nearing completion in Namibia, and ROTSE-IIId is slated to arrive in Turkey in the Fall of 2003.

Within a week of the first two ROTSE-III instruments becoming fully operational, the network was tested by the discovery with HETE-2 of GRB 030329. Spectroscopic observations found the source to lie at a redshift of $`z=0.168`$ (Greiner et al., 2003), which makes it the closest GRB to the Earth yet to be measured, with the possible exception of the anomalously faint GRB 980425 (Galama et al., 1999). This event provided the first direct spectral evidence for a supernovae associated with a high redshift GRB (Stanek et al., 2003; Hjorth et al., 2003; Kawabata et al., 2003). We report here on the light curve of GRB 030329 for the first two days of observations by ROTSE-IIIa & b.

## 2. Observations and Analysis

GRB 030329 was first observed by the instruments on board the HETE-2 satellite at 11:37:14.67 UTC. Ground analysis of the SXC data produced a localization that was reported in a Gamma-Ray Burst Coordinate Network (GCN) Notice 73 minutes after the burst (Vanderspek et al., 2003). Although the burst’s celestial location was above the Australian horizon at the time, the ROTSE-IIIa connection to the GCN was temporarily experiencing sporadic outages. It was therefore impossible for the ROTSE automated scheduler daemon to respond to the alert. However, electronic mail was also distributed to ROTSE team members, and a manual response by the ROTSE-IIIa telescope was initiated over the Internet within 15 minutes of the alert’s distribution, 88 minutes after the burst itself.

Clouds and intermittent rain interfered with observations and delayed analysis of the images, but once Peterson & Price (2003) and Torii (2003) identified the optical counterpart, we were able to extract a light curve from the setting source. Our earliest useful images began 92 minutes after the burst event (Rykoff & Smith, 2003). At 4.8 h after the burst, the air mass to the source grew too large for further observation.

ROTSE-IIIb (in Texas) began observations 14.8 h after the burst, and followed the source for a further 8.5 h. Once the source had risen again in Australia, contemporaneous observation from both ROTSE-III instruments was possible for 1.8 h. The ROTSE-III automatic scheduler normally plans observations of fading GRB afterglows to occur at logarithmically increasing intervals, leading to gaps in coverage at the start of each night’s observing. As it became clear how bright this source was, we overrode the default schedule to instruct the instruments to look at nothing else as long as GRB 030329 was above the local horizon.

The ROTSE-III images were bias-subtracted and flat-fielded in the standard way. For images recorded later than 10 h after the burst, we co-added sets of ten frames before deriving a light curve. We applied SExtractor (Bertin & Arnouts, 1996) with aperture photometry to these co-added images to detect objects and derive their intensities. We compared these intensities to the USNO A2.0 catalog to derive an approximate R-band magnitude zero-point offset for the field. Zero-point offsets among the ROTSE images are determined through application of a relative photometry procedure. Absolute calibration to a V-magnitude equivalent (dubbed $`V_{\mathrm{ROTSE}}`$) is derived by comparison of several nearby stable stars to the $`UBVR_cI_c`$ photometry of Henden (2003).

The offset between unfiltered $`V_{\mathrm{ROTSE}}`$ observations and standard V magnitudes is a function of color. Without color information for the early-time afterglow, we have not applied any color correction to these data. However, from 6.0–14.4 h after the burst, during observations reported by Burenin et al. (2003), the afterglow color is rather blue, with $`VR0.3`$. This color remained stable through a break in the lightcurve. If the afterglow maintained this color throughout the ROTSE-III observations, it would be appropriate to further adjust the $`V_{\mathrm{ROTSE}}`$ magnitudes reported here to be approximately 0.5 magnitudes brighter.

GRB 030329 was also observed with a time resolution of 1-second using a CCD time-series photometer mounted at the prime focus position of the McDonald Observatory 2.1-m telescope. The photometer “Argos” (Nather & Mukadam, in preparation) operates in frame transfer mode, and thus allows rapid exposures as short as one second. The time-series photometry started 15.87 hours after the burst trigger, and continued for 4290 seconds in essentially photometric conditions. No filter was used, so the spectral response of the instrument was dominated by the combination of a typical CCD sensitivity and the Earth’s atmosphere: the transmission is about 70% in the wavelength range 560–640 nm, and greater than 50% in the wavelength range 450–800 nm.

Typical ADU count rates for the optical transient were 17,000 per second, while the sky contributed about 10,500 counts per second in the chosen software aperture. There were two comparison stars in the 2.8 arcminute field of view of the instrument, and sky-subtracted light curves for each object were obtained by using an optimized circular software aperture. A corrected light curve for the target object was then generated by dividing its light curve to the sum of the two comparison star light curves (which have been normalized to unity).

## 3. Results

The full light curve from both ROTSE-IIIa & b is shown on a logarithmic scale in the top panel of Figure 1. Figure 2 shows two segments of the light curve on expanded scales. In this section, we describe the behavior of the afterglow during each observing sequence.

The poor observing conditions during the first observing run (from 1.5 to 4.8 h after the burst) introduce large systematic errors into the photometry, but there is no evidence for any significant deviations from a power law decay with an index of 0.88$`\pm `$0.05, a value not atypical for early GRB afterglows. The earliest detections reveal an unfiltered $`V_{\mathrm{ROTSE}}`$ magnitude of 12.55, over 100 times brighter than any other afterglow to date to be observed at this epoch (Price et al., 2003), consistent with the very low redshift for the object. The latest points in this sequence, recorded after an hour during which ROTSE-IIIa automatically ceased operation due to a light rain, are also consistent with the same decay curve.

When ROTSE-IIIb commenced observations at 14.8 h, the intensity of GRB 030329 was decaying much more rapidly, consistent with a power law index of 1.85$`\pm `$0.07 and strongly inconsistent with the earlier curve observed with ROTSE-IIIa. If we assume a single, sharp transition from one decay mode to the other, these data predict a transition time of $``$12–14 h after the burst. Over the course of the 8.5 h observing sequence with ROTSE-IIIb on 30 Mar 2003, the decay slope steepened further, demanding the introduction of a third power law with decay index of 2.25$`\pm `$0.08. The best-fit transition time to this decay slope is at 19.0 h after the burst.

The interval during which GRB 030329 was observed by the 2.1-m coincided with the early part of this decay phase, before the break at 19 h. The corrected light curve at 1-s resolution shows the optical transient to be fading smoothly. The fitted power law decline has an index of 1.98$`\pm `$0.02. This power law fits well to both the beginning and ending 500-second time intervals with the same normalization to within 0.2%, thus indicating that there were no departures from the overall power law decline. The largest deviation from the best fit power law was one point that was low by only $`4.4\sigma `$. Analysis of binned light curves reveals no significant fluctuations from the power law on any longer time scales. The computed discrete Fourier transform did not reveal any power out to the Nyquist frequency of 0.5 Hz up to a limit of 2 millimag (i.e., 0.2% change in intensity). In all, we found no deviations from a simple power law decline.

During the interval when both ROTSE-III instruments could observe the source simultaneously, our calibration procedure yielded consistent results. However, ROTSE-IIIa subsequently observed the source intensity to level off and remain constant for at least three hours. When the source again became visible to ROTSE-IIIb in Texas, it was found to still have the same intensity as when last seen by ROTSE-IIIa. It shortly commenced a rapid decay phase, consistent with a single power law of index 2.34$`\pm `$0.08. This decay rate is consistent with the decay slope observed at the end of the previous night. The decay phase lasted for approximately six hours. The source then entered a second flattening phase that was still continuing when observations ceased for the night two hours later111The ROTSE-III light curve is stored in electronic form at http://www.rotse.net/transients/grb030329.

The decay curve for GRB 030329 over its first two days is clearly inconsistent with a single power-law function. We find that at least four power laws plus two intervals of constant intensity are necessary to adequately reproduce the behavior of the light curve. Three of these functions are plotted as broken lines in Figure 2. The decay indices as well as the intervals over which each power-law segment applies are given in Table 1.

## 4. Discussion

The world-wide distribution of the ROTSE-III network allows for unprecedented coverage of the early decay curve of GRB afterglows, and in the case of GRB 030329, the afterglow was above the horizon for 56% of the first 47 h after the event. Since the ROTSE-III network was only half complete at the time of GRB 030329, there are still large gaps in the derived light curve when the source was not visible to either telescope. However, at least two groups have published observations of the afterglow’s behavior during the time on the first day when neither ROTSE-IIIa nor ROTSE-IIIb could observe it.

Uemura et al. (2003) observed GRB 030329 with a variety of telescopes from 1.2 h to $`14`$ h after the event, and they report three breaks in the light curve during this time: from 0.74$`\pm `$0.02 to 0.95$`\pm `$0.01 at $`t=2.04`$ h, to 0.65$`\pm `$0.04 at $`t=3.9`$ h, and again to 1.16$`\pm `$0.01 at $`t=5.4`$ h. This characterization is consistent with the ROTSE-IIIa data, which have neither sufficient coverage nor small enough errors to distinguish between this and the best-fit single power law reported above in Section 3. Burenin et al. (2003) present observations taken with the 1.5 m Russian-Turkish Telescope shortly after the time of this last transition and give a slope of 1.19$`\pm `$0.01 from 6.2 h to 14.4 h after the burst, after which the slope changes to $`1.9`$, which is consistent with the ROTSE-IIIb observations at 14.8 h. The Japanese and Turkish programs, therefore, neatly fill in the gap between ROTSE-III observations, and the joint data set provides evidence for five changes in slope decay index within the first 20 hours after the burst.

No GRB afterglow has been observed nearly as intensely as GRB 030329, and for most other bursts only a single steepening in the light curve decay has been reported. These breaks have been taken as evidence for deceleration of a relativistic jet: as the jet decelerates, the beaming angle must increase, and when the beaming angle exceeds the jet angle, the intensity of visible light is expected to decrease more rapidly (Rhoads, 1997; Sari, Piran, & Halpern, 1999; Perna, Sari, & Frail, 2003). Multiple breaks within a single decay curve are more difficult to interpret in this paradigm, but Wei & Jin (2003), for example, have suggested that they may result from off-axis viewing geometry or perhaps non-uniform structure within the jet. A detailed analysis of what geometry five decay breaks might demand is beyond the scope of this paper, but the data clearly are not consistent with the simplest models and demand caution when trying to interpret $`t0.5`$ d as the jet break (Granot, Nakar, & Piran, 2003; Burenin et al., 2003; Price et al., 2003).

The other striking feature of the afterglow is the deviation from its rapid decay, beginning at $`23`$ h after the burst. In the ROTSE-III observations, the intensity levels off and remains constant for several hours. A second decay phase begins $`39`$ h after the burst. It is unlikely that this episode is related to the emergence of a supernova in the light curve, as the first spectral evidence for a supernova associated with GRB 030329 emerged only after $`7`$ d (Stanek et al., 2003). At one day after the event, the optical intensity should be dominated by the burst afterglow. Furthermore, the light curve flattens yet again at $`44`$ h. Wei & Jin (2003) argue that geometric factors could produce a flattening of the light curve at about this time, but it is hard to understand the repeating nature of the phenomenon within that paradigm. On the other hand, Granot, Nakar, & Piran (2003) have interpreted these episodes as refreshed shocks produced by further ejection events from the central engine catching up with the decelerating blast wave. In this interpretation, the central engine remains an active part of the GRB process well into the afterglow phase, further underscoring the need for early and thorough observations of the entire burst event if the complex physical interactions are ever to be teased out.

This work has been supported by NASA grants NAG 5-5281 and F006794, NSF grant AST 01-19685, the Michigan Space Grant Consortium, the Australian Research Council, the University of New South Wales, the University of Texas, and the University of Michigan. DAS is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST 01-05221. Work performed at LANL is supported by NASA SR&T through Department of Energy (DOE) contract W-7405-ENG-36 and through internal LDRD funding.