The effects of gravitational darkening on the thermal structure of Be star disks of differing densities are systematically examined. Gravitational darkening is the decrease of the effective temperature near the equator and the corresponding increase near the poles of a star caused by rapid rotation. We also include the rotational distortion of the star using the Roche Model. Increasing the disk density increases the optical depths in the equatorial plane, resulting in the formation of an inner cool region near the equatorial plane of the disk. High rotation rates result in disks that have temperatures similar to those of a denser disk, namely cooler overall. However the effect of increasing rotation produces additional heating in the upper disk due to the hotter stellar pole. Cool regions in the equatorial plane normally associated with high density are seen in low density models at high rotation rates. Gravitational darkening increases the amount of very cool and very hot material in the disk and decreases the amount of disk material at moderate temperatures. We also present models which study the effect of gravitational darkening on hydrostatically-converged disks, in which the temperature structure is consistent with vertical hydrostatic equilibrium. Because the equatorial regions become cooler, hydrostatically converged models that include gravity darkening have smaller vertical scale heights, and $H/R$ is smaller by as much as 56% near $v_{\rm crit}$. Finally we explore differences in disk temperatures when alternate formulations of gravitational darkening, which lower the temperature difference between the pole and the equator, are used.

The effect of gravitational darkening on models of the thermal structure of Be star disks is systematically studied for a wide range of Be star spectral types and rotation rates. Gravitational darkening causes a reduction of the stellar effective temperature towards the equator and a redirection of energy towards the poles. It is an important physical effect in these star-disk systems because the photoionizing radiation from the central B star is the main energy source for the disk. We have added gravitational darkening to the bedisk code to produce circumstellar disk models that include both the variation in the effective temperature with latitude and the non-spherical shape of the star in the calculation of the stellar photoionizing radiation field. The effect of gravitational darkening on global measures of disk temperature is generally significant for rotation rates above 80% of critical rotation. For example, a B0V model rotating at 95% of critical has a density-averaged disk temperature \approx 2500 K cooler than its nonrotating counterpart. However, detailed differences in the temperature structure of disks surrounding rotating and non-rotating stars reveal a complex pattern of heating and cooling. Spherical gravitational darkening, an approximation that ignores the changing shape of the star, gives good results for disk temperatures for rotation rates less than \approx 80% of critical. However for the highest rotation rates, the distortion of the stellar surface caused by rotation becomes important.

A popular model for the circumstellar disks of Be stars is that of a geometrically thin disk with a density in the equatorial plane that drops as a power law of distance from the star. It is usually assumed that the vertical structure of such a disk (in the direction parallel to the stellar rotation axis) is governed by the hydrostatic equilibrium set by the vertical component of the star's gravitational acceleration. Previous radiative equilibrium models for such disks have usually been computed assuming a fixed density structure. This introduces an inconsistency as the gas density is not allowed to respond to temperature changes and the resultant disk model is not in vertical, hydrostatic equilibrium. In this work, we modify the {\sc bedisk} code of \citet{sig07} so that it enforces a hydrostatic equilibrium consistent with the temperature solution. We compare the disk densities, temperatures, H$\alpha$ line profiles, and near-IR excesses predicted by such models with those computed from models with a fixed density structure. We find that the fixed models can differ substantially from the consistent hydrostatic models when the disk density is high enough that the circumstellar disk develops a cool ($T\lesssim10,000 $K) equatorial region close to the parent star. Based on these new hydrostatic disks, we also predict an approximate relation between the (global) density-averaged disk temperature and the $T_{\rm eff}$ of the central star, covering the full range of central Be star spectral types.