© ESO, 2016

1. Introduction

When the components in close binary systems are rotationally and/or tidally distorted, the fluxes are not uniformly distributed over their surfaces. The same applies, mutatis mutandis, to fast, isolated rotators. In 1924, von Zeipel (1924) quantified the flux distribution in distorted stars and found that ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\propto }{\mathit{g}}^{{\mathit{\beta }}_{\mathrm{1}}}$, where g is the local gravity, T_eff is the effective temperature, and β₁ is the gravity-darkening exponent (GDE). For the case of configurations in radiative equilibrium, β₁ = 1.0. For some decades, this was the only available value of theoretical GDE to be compared with observations. The next step in the theoretical derivation of GDE was taken by Lucy (1967). In this paper Lucy, using the properties of an adiabatic in convective envelopes, derived an average value of β₁ equal to 0.32; this average value was derived for a small effective temperature range: 5800–6500 K. Later, Webbink (1976) and Sarna (1989) confirmed Lucy’s results, although Anderson & Shu 1977 found that there was no relation between local gravity and fluxes for late-type stars, that is, β₁ = 0.0 for contact binaries.

In 1998 Claret introduced a numerical method to compute the GDE that embraces stars with radiative and convective envelopes as a function of the stellar mass, effective temperature, radius, chemical composition, and time. For hotter configurations, the von Zeipel value is recovered while for the colder configurations the results reported by Lucy were in general confirmed. Also, a smoother transition zone was achieved between the radiative and convective energy transport mechanisms. More recent theoretical investigations on the temperature distribution over distorted stellar surfaces were carried out, for instance by Claret (2000, 2012), Espinosa Lara & Rieutord (2011, 2012), McGill et al. (2013), and Claret (2015).

With respect to observations, semi-empirical values of the GDE were obtained only for a few cases in eclipsing binaries. The situation changed with the pioneering work by Rafert & Twigg (1980). In that paper, the authors inferred the semi-empirical GDE for several detached and semi-detached systems, covering a wide range of effective temperatures. Later, some authors extended this study, in particular, the investigations carried out by Pantazis & Niarchos (1998), Niarchos (2000), and Djurasevic et al. (2003, 2006). These semi-empirical GDE are shown in Fig. 2 by Claret (2015). In that paper, the author derived an equation for the GDE with a perturbation theory. When compared with observations, there is a scattering around the von Zeipel theoretical prediction for the hotter systems. Such a scattering can be explained by the mentioned equation assuming a moderate degree of differential rotation. For those hotter systems with β₁ lower than 1.0, a possible explanation could be the presence of an accretion disk that is rotating faster than the stellar surface. With respect to colder systems, the small values of β₁ were interpreted as due to the role of convection in the stellar upper layers and also to the change of the source of thermonuclear energy from CNO cycle to proton-proton chain (see Sect. 4.2 for a more detailed discussion). The resulting equation for β₁ is also capable of predicting the drop-off/transition zone around log T_eff ≈ 3.9, which is in good agreement with the observational data provided by Rafert & Twigg (1980) and Djurasevic et al. (2003, 2006).

The case of isolated, fast rotators is somewhat different. Reliable observations of fast rotating stars was only possible after the advent of long-baseline optical/infrared interferometry, which provides high resolution data. In 2005 McAlister and collaborators detected for the first time the GDE for a single star, α Leo (an early-type star), whose β₁ was in good agreement with the classical value, that is, 1.0. The ratio ω/ω_crit for this star is around 0.9, where ω is the angular velocity and ω_crit is the break-up velocity (Table 1). Some years later, Che et al. (2011) using the CHARA array and the instrument MIRC at the band H, derived ${\mathit{\beta }}_{\mathrm{1}}\mathrm{=}\mathrm{0.75}{\mathrm{2}}_{-0.116}^{\mathrm{+}\mathrm{0.048}}$ for the same star, which is significantly different from von Zeipel’s value. In addition to the deviation of the GDE for α Leo, a small value was also reported of ${\mathit{\beta }}_{\mathrm{1}}\mathrm{=}\mathrm{0.58}{\mathrm{4}}_{-0.028}^{\mathrm{+}\mathrm{0.052}}$ for the case of β Cas. This small value was expected for β Cas since its equatorial and polar temperatures locate the system near the drop-off/transition zone predicted by Claret (1998, 2000). For details, see Fig. 11 in Che et al. (2011). However, the accuracy of the value of β₁ for α Leo (radiative envelope) revealed that the deviation from the value of von Zeipel was real for this fast rotating star. From an observational point of view, the study by Che et al. represents a turning point in the empirical study of the temperature distribution of rapidly rotating stars. After the work by McAlister (2005), gradually more and more data for fast, isolated systems were obtained, as for example those by van Belle et al. (2006), Aufdenberg et al. (2006), Peterson et al. (2006), Monnier et al. (2007), Zhao et al. (2009), Che et al. (2011), and Domiciano de Souza et al. (2012, 2014). The more recent papers dealing with hotter systems and listed in Table 1 confirm the tendency reported by Che et al. (2011).

Another interesting application of the GDE is related to the behaviour of massive rotating stars that are losing mass. Because of the influence of gravity darkening the associated wind is anisotropic; at the pole the wind is fast and dense, while it is thinner and slower at the equator. This affects not only the total mass-loss rate but also the angular distribution of the wind. The loss of angular momentum controls the fate of the star; if the stellar wind presents spherical symmetry, the loss of angular momentum is larger than the anisotropic case and the star achieves a slow rotation rate and, as a consequence, the formation of a long-duration gamma ray burst (LGRB) is not achieved. In the case of anisotropic winds, the star retains enough angular momentum and potentially explodes as a LGRB. The degree of spherical symmetry of the wind depends on the GDE. Often the classical value by von Zeipel is adopted to evaluate the characteristics of the wind. However, as shown in this paper, there is theoretical evidence of deviations of the von Zeipel theorem (there is also observational evidence, as we have seen) and this may change the wind structure and the history of massive stars.

The mentioned deviations from the von Zeipel’s value found in hotter systems pose a challenge to the theory of the gravity-darkening exponents. In this paper, we present new insight into the problem of temperature distribution over distorted stellar surfaces. In this approach, we make use of the numerical method described in Claret (1998, 2000, 2012), but inserting the dependence of the theoretical GDE with the optical depth as a novelty. The paper is structured as follows: Sect. 2 is devoted to present theoretical evidence of deviations of the von Zeipel theorem for the case of stars with radiative envelopes. Section 3 summarizes the properties of the main sources of opacities as a function of the spectral type and their relation with the infrared filters. In Sect. 4 we discuss the influence of convection on the flux distribution and derive an equation relating the convective efficiency with the GDE. Finally, we introduce the numerical method to compute the theoretical GDE as a function of the optical depth and compare the resulting calculations with accurate semi-empirical data for six fast rotators.

2. Theoretical evidence of deviations of the von Zeipel theorem

The von Zeipel theorem is valid for conservative rotation laws and when the diffusion approximation is used as transfer energy equation. However, if a different transfer energy equation is adopted, for example a grey atmosphere, one could expect some deviation from the result by von Zeipel. Recently Claret (2012), using a variant of the numerical method mentioned in the Introduction, found significant deviations from the classical theorem by von Zeipel at the upper layers of a distorted hot star. In particular, the case of α Leo discussed before, was satisfactorily explained as shown in Fig. 1 in Claret (2012). In that paper, more elaborated transfer equations were used, including realistic stellar atmosphere models (ATLAS). In this section we confirm the deviations of the von Zeipel theorem at the upper layers via a simple analytical method.

Before proceeding, it is important to establish the range of validity of the present calculations. The divergence of radiation flux in an atmosphere of a distorted star in hydrostatic equilibrium is not vanishing throughout and this causes meridional circulation, depending on the viscosity. When we calculate the divergence of radiation flux in this atmosphere (assuming, for example, the Eddington approximation) one finds that it vanishes only in two situations: a) in deepest layers or b) for small distortions. We restrict our calculations to small/moderate distortions for which the meridional circulation can be neglected. It is hard to establish a precise limit of validity of our calculations in the framework of our numerical method (described in Sect. 4.2). In future works, we plan to investigate the role of the meridional currents and accurately derive the range of validity of the present calculations. Another characteristic of our method is that we are dealing with plane-parallel geometry for which the effects of sphericity are negligible.

Fig. 1 Convective efficiency and gravity-darkening exponent in the HR diagram. The models were computed with the GRANADA code and the masses range from 0.8 up to 3.16 M⊙. Dashed line indicates the locus of the f = 0.5 (Eq. (10)). This line roughly divides the HR diagram in envelopes with high convective efficiency (right) and low convective effciency envelopes (left). The critical log Teff at the ZAMS is around 3.9.

To evaluate the effects of an atmosphere on gravity-darkening calculations, keeping in mind the above discussion, we consider two points in a distorted configuration that are on the same equipotential. Using the hydrostatic differential equation for these points, we have g(μ)τ(μ,Ψ) = g(μ_o)τ(μ_o,Ψ) ≡ g_oτ_o, where g is the local gravity, Ψ is the total potential (rotational one included), τ is the optical depth, μ is given by cos(θ), and θ is the angle between the radiation field and the z axis. The subscript o indicates a reference point. In addition, we have B(τ(μ)) = B_o(τ_o(μ_o)), where B denotes the Planck function. We adopt for the calculation of the fluxes the quadrature formula ${\mathit{E}}_{\mathit{n}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\approx }{\sum }_{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}^{\mathit{n}\mathrm{-}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathit{t}\mathit{/}{\mathit{\mu }}_{\mathit{i}}}$ as well as two properties of this approximation, ${\sum }_{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1}\mathit{/}\mathrm{3}$ and ∑ _iα_iμ_i ≈ 1/2, where α_i are the weights (Reiz 1950; Unno 1962). For simplicity, we use the Eddington approximation for a stellar atmosphere model. We allow F and F_o to be the fluxes in a generic optical depth and in a reference point, respectively. Starting with B(t) = 3/4F_o(tg/g_o + 2/3) for the flux at a given optical depth τ, we obtain $\begin{array}{ccc}\mathit{F}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathrm{\int }}_{\mathit{\tau }}^{\mathrm{\infty }}\mathit{B}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{E}}_{\mathrm{2}}\mathrm{(}\mathit{t}\mathrm{-}\mathit{\tau }\mathrm{)}\mathrm{d}\mathit{t}\mathrm{-}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\tau }}\mathit{B}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{E}}_{\mathrm{2}}\mathrm{(}\mathit{\tau }\mathrm{-}\mathit{t}\mathrm{)}\mathrm{d}\mathit{t}& & \\ & & \hspace{0.17em} \mathrm{\approx }{\mathit{F}}_{\mathrm{o}}\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\mu }}_{\mathit{i}}}\mathrm{+}\mathrm{3}{\mathit{F}}_{\mathrm{o}}\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathrm{o}}}\right)\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}\\ & & \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{F}}_{\mathrm{o}}\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathrm{o}}}\right)\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathit{/}{\mathit{\mu }}_{\mathit{i}}}\mathit{.}\end{array}$(1)Therefore, for τ = 0 we have $\begin{array}{ccc}\mathit{F}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{o}}}{\mathrm{2}}\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathrm{o}}}\right)\mathrm{+}{\mathit{F}}_{\mathrm{o}}\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}\mathrm{\Rightarrow }{\mathit{\beta }}_{\mathrm{1}}\mathrm{\approx }\mathrm{0.5}\mathit{,}& & \end{array}$(2)which deviates significantly from the classical value (1.0). In other words, the von Zeipel theorem is no longer valid in the outermost layers of a distorted star. Similar results were already obtained by Osaki (1966) and numerically by Claret (2012, Fig. 1) for the case of the exact grey atmosphere and non-grey atmospheres. For the deepest layers Eq. (1) yields $\begin{array}{ccc}\underset{\mathit{\tau }\mathrm{\to }\mathrm{\infty }}{\lim }\mathit{F}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{3}{\mathit{F}}_{\mathrm{o}}\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathrm{o}}}\right)\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{i}}{\mathit{\mu }}_{\mathit{i}}^{\mathrm{2}}\mathrm{=}{\mathit{F}}_{\mathrm{o}}\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathrm{o}}}\right)\mathrm{\Rightarrow }{\mathit{\beta }}_{\mathrm{1}}\mathrm{\approx }\mathrm{1.0.}& & \end{array}$(3)The classical von Zeipel theorem is only recovered (β₁ = 1.0) for large optical depths where the diffusion approximation is valid. Again, this was previously verified through numerical methods by Claret (2012). If one uses, for example, the exact solution for the grey atmosphere, the final result, concerning β₁, is unchanged. An important feature emerges from Eqs. (1)–(3): the theoretical gravity-darkening exponent depends on the optical depth, that is, β₁ ≡ β₁(τ). In Sect. 4.2 we extend these results via our numerical method. Therefore, one would expect that, in observational terms, the GDE measurements ought to depend on the spectral type of target stars. In the next section, we briefly discuss the physical conditions in the stellar envelopes and their relationship with the infrared observations.

3. Some considerations on opacities and photometric systems

In previous years, a considerable effort was devoted to the investigation of very fast rotators. The observations were carried out using optical/infra-red interferometers mainly in the H and K-bands (see for example McAlister et al. 2005 and Domiciano de Souza et al. 2005). An important characteristic of the H-band, for example, is that it is located in a region where the transparency of the Earth’s atmosphere is high (Wing & Jorgensen 2003).

When one performs observations using a given photometric system, an old and frequent question arises: How deep in a star are we observing? Not only are the characteristics of a given photometric system important to obtain information on the target stars, but also their spectral characteristics are relevant to estimate the depth of the layers we are observing. This amalgam requires knowledge of the physical conditions of the stellar envelopes, mainly the equation of state and monochromatic opacity. Schematically, the main sources of opacity as a function of the spectral type can be summarized as follows:

• O-B ⇒ Electron scattering, HeI, HeII(b-f), and H (f-f)

• B-A ⇒ HI (b-f), HI(f-f), HeII(b-f), electron scattering

• A-F ⇒ similar contributions of HI (b-b), HI(f-f), H⁻ (b-f), H⁻ (f-f)

• G-K ⇒ HI(b-f), H⁻ (b-f), H⁻ (f-f), Rayleigh scattering

For a star of spectral type G, the opacity due to H⁻ (b-f) dominates in the optical region (maximum around 8500 Å), while H⁻ (f-f) is important in the infrared. However, there is a local minimum of the total opacity right where the H-band is located (≈16.000 Å). Therefore, the observations in the H-band allow us to access warmer and deeper layers in the stellar atmospheres of these kind of stars. In late A-type stars the physical conditions are different, although H⁻ still has some contribution to the total opacity in the optical and infrared regions. However, H (b-f) begins to dominate and the total opacity at the H-band is larger than in the previous discussed case. In the external layers of a B-type star the total opacity is dominated by H (b-f). In the interval where the H-band is located, the total opacity can achieve values that are several times larger than that for G-type stars. Summarizing and roughly speaking, the infrared observations only permit us observe the uppermost layers of hotter stars, while we can go deeper in late-type stars. In the following, we see how these properties are, in a way, correlated with the adopted numerical method to compute the theoretical β₁(τ).

4. Theoretical gravity darkening as a function of optical depth: Numerical method and comparison with observational data

Given that Eqs. (2) and (3) are consistent with the numerical method introduced in Claret (1998, 2000, 2012), we adopt this method to investigate how the GDE depends on the optical depth. We briefly describe our method, which is based on the triangle strategy, often used in stellar evolution calculations and introduced to save computational time (Kippenhahn et al. 1967). Consider a reference model with surface g_o and effective temperature T_effo in the Hertzsprung-Russell (HR) diagram. Three envelopes with the following parameters are then computed around this reference model: (g_o + Δg₁, T_effo + ΔT_eff1), (g_o - Δg₂, T_effo – ΔT_eff2), and (g_o + Δg₃, T_effo – ΔT_eff3). If the next point of the evolutionary track to be computed is still within this triangle, the previous envelope integrations can be used with enough accuracy (provided that the triangle is sufficiently small). We adopted the non-rotating and 1D version of the GRANADA code to generate the envelopes/atmospheres. To simulate distorted models, we take as our reference a model characterized by (g_o, T_effo) and apply a similar triangle strategy, but we increase the number of trial triangles and distribute them strategically around the reference model for which we are interested in computing the GDE. In this way we have several envelopes characterized by different gs and T_effs. To represent a distorted star with different flux distribution over the surface, we use these envelope models with different temperature distributions, but present the same pressure-temperature relationship at a given shell. This condition is achieved by adjusting the values for each T_effo ± ΔT_eff adopted in the modified triangle strategy until these envelopes have the same pressure temperature of the reference model at a given point. In general, the outlet effective temperatures are different from the inlet effective temperatures. A simplified sketch of the method is shown in Fig. 2 (Claret 2000). Finally, differentiating the properties of the envelopes, ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\propto }{\mathit{g}}^{{\mathit{\beta }}_{\mathrm{1}}}$, we obtain the GDE for each point of the evolutionary track. In this paper mentioned above, we adopted large τ for all models, i.e. we imposed the boundary conditions where the diffusion equation is valid. The novelty of the present method is that, as discussed in Sect. 1, we extended our calculations to the uppermost layers and, therefore, the GDE is now derived as a function of the optical depth. We treat this point in more detail in Sect. 4.2.

4.1. Gravity-darkening exponents and convective efficiency

Before entering the calculations of β₁(τ) themselves we make some general considerations about our numerical method as well as the relationship between GDE and convection. First, we examine how the pressure-temperature slope depends on the depth, which is an important issue concerning our method. Consider an opacity law of the form κ ≈ κ_oPⁿT^s, where P is the pressure, T the temperature, κ_o is a constant, and n> 0 and s< 0. In the envelope, the enclosed mass and luminosity are constant and equal to their surface values. Therefore, the thermal differential equation for the envelope can be written as$\begin{array}{ccc}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{T}}\mathrm{=}{\mathit{C}}_{\mathrm{1}}\frac{{\mathit{T}}^{\mathrm{3}}}{\mathit{\kappa }}\mathit{,}& & \end{array}$(4)where C₁ is a constant. Integrating this equation from the effective temperature (T_eff, P_eff) up to a generic point (T,P) we obtain $\begin{array}{ccc}{\mathit{P}}^{\mathrm{1}\mathrm{+}\mathit{n}}\mathrm{-}{\mathit{P}}_{\mathrm{eff}}^{\mathrm{1}\mathrm{+}\mathit{n}}\mathrm{=}\frac{{\mathit{C}}_{\mathrm{1}}}{{\mathit{\kappa }}_{\mathrm{o}}}\frac{\mathrm{1}\mathrm{+}\mathit{n}}{\mathrm{4}\mathrm{-}\mathit{s}}\mathrm{(}{\mathit{T}}^{\mathrm{4}\mathrm{-}\mathit{s}}\mathrm{-}{{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}\mathrm{-}\mathit{s}}}^{\mathrm{)}}\mathit{.}& & \end{array}$(5)For positive values of 1 + n and 4−s (for example, the Kramer’s law) and considering very deep regions inside the star, the relationship pressure temperature is practically independent of the surface values, in particular of the effective temperature. This is why our original numerical method (large τ) recovers the von Zeipel value for all radiative envelopes and is independent of T_eff (see for example Fig. 1 in Claret 2000). This is also consistent with Eq. (3). However, near the surface, the differences between P and P_eff and T and T_eff are sufficiently small and the pressure temperature relationship depends mainly on T_eff and does not necessarily reproduce the von Zeipel value (Eq. (2)). This is the case of hotter stars we analyse below.

For convective layers Eq. (4) is no longer valid. In the deep layers of convective envelopes, where the density is relatively high, the slope of pressure temperature is given by the adiabatic gradient and the resulting GDE is smaller than 1.0 and presents a local maximum around log T_eff ≈ 3.7 (Fig. 1, Claret 2000). This behaviour of the GDE for convective envelopes is also supported by calculations via transfer equations that are more elaborated. In more external layers, the analysis is more complicated since convection is inefficient and stratification is overadiabatic; is it only possible to determine the slope of pressure temperature by adopting a suitable convection theory.

Fig. 2 Rosseland optical depth, where the boundary conditions are imposed, as a function of the effective temperature adopted in the calculation of the gravity-darkening exponents (ZAMS models). We note the drop-off for radiative envelopes (log Teff ≈ 3.9).

Still, concerning the role of convection in the surface temperature distribution, it is interesting to try to relate the convective efficiency with the GDE. To do this, we write the convective flux as $\begin{array}{ccc}{\mathit{F}}_{\mathrm{c}}\mathrm{=}\mathit{\eta \rho }{\mathit{v}}_{\mathrm{s}}{\mathit{c}}_{\mathrm{P}}\mathit{T,}& & \end{array}$(6)where ρ is the density, v_s is the local velocity of sound, c_P is the specific heat at constant pressure, and η is given by $\begin{array}{ccc}\mathit{\eta }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\overline{\mathit{v}}}{{\mathit{v}}_{\mathrm{s}}}\frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}\mathit{,}& & \end{array}$(7)where $\overline{\mathit{v}}$ is the mean convective velocity along a mixing length and ΔT is the excess of temperature of a rising element over the mean temperature of the surrounding. Defining the convective efficiency f as the ratio of F_c to the total flux (Cox & Giuli 1968) now we obtain $\begin{array}{ccc}\mathit{f}\mathrm{=}\mathit{A\eta }{\left(\frac{\mathit{r}}{\mathit{R}}\right)}^{\mathrm{2}}{\left(\frac{\mathrm{3}{\mathrm{\Gamma }}_{\mathrm{1}}}{\mathrm{5}{\mathit{\mu }}_{\mathrm{1}}\mathit{\beta }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{\mathrm{2}{\mathit{c}}_{\mathrm{P}}{\mathit{\mu }}_{\mathrm{1}}\mathit{\beta }}{\mathrm{5}}\right){\mathit{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{{\mathit{g}}^{\mathrm{2}\mathit{/}\mathrm{3}}}{{\mathit{T}}_{\mathrm{eff}}^{\mathrm{27.4}\mathit{/}\mathrm{3}}}\right)\mathrm{·}& & \end{array}$(8)In the above equation, A is a constant, R the star radius, r the radial coordinate, μ₁ the mean molecular weight, β is the ratio of gas to total pressure, and Γ₁ is the adiabatic exponent given by (dln P/dln ρ)_ad. The exponent for the effective temperature in Eq. (8) depends on the adopted opacity law approximation. The exponents of g and T_eff were evaluated assuming that we are dealing with an ideal gas equation of state and that the negative hydrogen ion is the main source of opacity, κ ≈ κ₁ρ^0.5T^7.7, where κ₁ is a constant. We note that this is a slightly different approximation for the opacity from that adopted above. The large exponent for T_eff in Eq. (8) drives the behaviour of the convective efficiency and predicts a very sharp transition in the HR diagram between the inefficient-efficient regimes.

On the other hand, Eq. (8) is also important to evaluate the role of the convective efficiency in the surface temperature distribution. Consider the locus where f is constant. Inserting typical values for a convective envelope in the above equation at r/R ≈ 1.0, as for example, T = 10⁴ K, Γ₁ = 5/3, μ₁ = 1.0, β = 1.0, we have $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\propto }{\mathit{g}}^{\mathrm{4}\mathit{/}\mathrm{13.7}}\mathrm{=}{\mathit{A}}_{\mathrm{1}}{\mathit{g}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathit{,}& & \end{array}$(9)where A₁ is a constant. From Eq. (9), we obtain β₁ ≈ 0.30 that is in very good agreement with the results of Lucy and those by Claret (1998, 2000) for stars with convective envelopes. Therefore, the locus of constant convective efficiency in the HR diagram can be written in terms of the gravity-darkening exponents. Using Eq. (9) and assuming that the luminosity is approximately proportional to M⁴, we can derive an equation for constant convective efficiency in the HR diagram $\begin{array}{ccc}\mathit{L}\mathrm{=}{\mathit{A}}_{\mathrm{2}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{-}\frac{\mathrm{4}}{\mathrm{3}}\mathrm{(}\mathrm{4}\mathit{/}{\mathit{\beta }}_{\mathrm{1}}\mathrm{-}\mathrm{4}\mathrm{)}}\mathit{,}& & \end{array}$(10)where A₂ is a constant that depends on the efficiency we are dealing with. In Fig. 1 the line labelled f = 0.5 represents the locus in the HR diagram showing this efficiency in the envelopes for β₁ = 0.3. The dashed line in Fig. 1 roughly divides the HR diagram in envelopes with high convective efficiency (right) and inefficient convective efficiency (left); the further to the right the models are from this line the greater the efficiency, while models located on the left present smaller convective efficiency. As commented, as a result of the large exponent of T_eff in Eqs. (8) and (9), the transition efficient-inefficient convection is very sharp and resembles the behaviour of β₁ × T_eff found in our previous studies on GDE. We note that the critical value of the temperature at the ZAMS (zero-age main sequence) is log T_eff ≈ 3.9.

4.2. The numerical method and comparison between β₁(τ) and observations of fast rotators

Returning to our numerical method, this procedure for late-type stars for which convection is dominant is also consistent with the works by Lucy (1967), Webbink (1976), Sarna (1989), and Claret (1998, 2000). The mentioned numerical method enables us to compute the GDE for radiative and convective envelopes and it also allows us to go as deep as necessary into the interior of the model. Moreover, the theoretical GDE can be computed as a function of the evolutionary status (mass, effective temperature, radius, chemical composition, and time), the changes in the chemical profile of the envelopes can be taken into account, and we can adopt more realistic atmosphere models as external boundary conditions. In addition, this method incorporates the change of the predominant source of thermonuclear energy from CNO cycle to proton-proton chain. In the region of the HR diagram where the main source of thermonuclear energy changes from CNO cycle to pp-chain occurs, as a consequence, a readjustment of the mass distribution of the stars. This variation leads to a well-established change of the slope log T_eff × log g and in log T_eff × log k_j (see Fig. 1 in Claret 2000), where k_j are the apsidal-motion constants of order j. These constants are derived by integrating the Radau equation for each stellar model from the centre up to the atmosphere and measuring the degree of mass concentration of a star. On the other hand, to compute the GDE we must have a good description of the shape of a distorted stellar configuration which, in turn, depends on k_j. Therefore, the variation in the slope log T_eff × log g (or log T_eff × log k₂) as a result of the change of the main source of thermonuclear energy is related to the GDE determination, not directly through the nuclear sources, but by how the star reacts under rotational perturbations through the apsidal-motion constant k₂.

To consider the physical conditions of stellar envelopes needed in our numerical method, we computed ZAMS stellar models covering the range of the observed effective temperatures (Table 1) via the code GRANADA. For details of the code, see Claret (2004). The basic input physics are X = 0.70, Z = 0.02, α_MLT = 1.68 and α_ov = 0.2. At the external layers of these models, we adopt the grey atmosphere approximation to investigate the behaviour of β₁(τ). A more elaborated atmosphere model, for example the ATLAS grids, should be preferred. However, the neighbouring points in log g and effective temperatures differ 0.5 and up to 2500 K, respectively, and this introduces some numerical noise due to the interpolation process in the atmosphere structures. Another limitation is that the maximum Rosseland optical depth of these atmosphere models is τ_Ross = 100, which restricts the application of our method to hotter systems. For our present purposes, we consider the grey atmosphere as satisfactory, at least as a first approach.

The structural readjustments due to the change of thermonuclear sources and to the onset of convection also modify, for example, the ratio between the electron pressure P_e and the total pressure, P_total. The resulting ratio presents a sharp fall just at the onset of convection (log T_eff ≈ 3.9) with direct consequences for the theoretical β₁(τ) calculations. This ratio is similar to the theoretical log T_eff × β₁ (see for example Fig. 1 in Claret 2000). Considering the implications of Eqs. (2) and (3), we can tentatively relate the theoretical β₁(τ) with the properties of the external layers of fast rotating stars or, more specifically, with their effective temperatures. As discussed before, our numerical method requires that the envelopes/atmospheres with different temperature distribution, representing distorted configurations, present the same temperature-pressure relationship as the reference model at a given point. Often, authors who adopted similar numerical methods established a point in the region where the helium ionization is complete. The corresponding points are a function of the effective temperature and a minor scale of the gravity. In this paper, we are interested in investigating the behaviour of the GDE at the upper layers. Therefore, the boundary condition is established in the zones where the ionization of the hydrogen is complete and the corresponding optical depth is designed by τ_RossBC. This means that τ_RossBC is larger for cold stars and smaller for hotter systems. Using the physical properties of the stellar envelopes/atmospheres of the ZAMS models generated by the GRANADA code, to compute β₁(τ) we infer a relationship between τ_RossBC and T_eff. Such a relationship is shown in Fig. 2 and Table 2. We note that τ_RossBC is also a function of the gravity of the models. The shape of this curve is similar to that for log T_eff × log k₂, which also presents a drop-off at the onset of convection (Claret 2000). The dependency of τ_RossBC with effective temperature is almost constant for T_eff> 4.2. However, the shape of the log T_eff × τ_RossBC relationship in this interval, and for smallest effective temperatures as well, could be even more pronounced than that shown in Fig. 2. The Rosseland optical depth scale is only a reference indicator because the corresponding opacities depend strongly on the completeness and/or on the accuracy of the line data adopted and these scales differ from author to author. In addition, the shape of the log T_eff × τ_RossBC relationship depends, for example, on the input physics and/or on the evolutionary status of the models. We note that we are not arguing that the observations in the bands H and K for cooler systems allow us to achieve τ_Ross ≈ 235 as shown in Fig. 2; that depth only indicates the region of the star where we establish the boundary conditions inherent to our numerical method, i.e. equal pressure × temperature relationship. In summary, to compute the values of the theoretical GDE as a function of the effective temperature to be compared with observations, we impose the boundary condition at the corresponding τ_RossBC listed in Table 2 (or Fig. 2) and apply our modified triangle strategy as described above. The resulting GDE is tabulated in the fourth column of Table 2.

Fig. 3 Theoretical gravity-darkening exponent as a function of the Rosseland optical depth. Model with 4 M⊙, log g = 4.293, and log Teff = 4.165. See text for more details.

As an example of the application of our numerical method, Fig. 3 shows the variation of the GDE with the Rosseland optical depth for a ZAMS model with 4 M_⊙, log g = 4.293, and log T_eff = 4.165. To compare with the results from Eqs. (2) and (3), we artificially vary the point where the boundary condition is applied, i.e. at τ_RossBC. The tabulated value for this model is τ_RossBC = 0.71. As expected, in deeper layers the von Zeipel result is reproduced but near the surface, for τ_Ross of the order of the unity or less, there is a systematic decreasing of β₁(τ), which is almost linear for log τ_Ross< 0.25. This result is consistent with Eqs. (2) and (3) and also with Fig. 1 by Claret (2012).

In spite of the problems (mainly numerical noise) of using ATLAS instead of grey atmospheres, we have performed some tests adopting these atmosphere models. As the maximum value of τ_Ross for these models is around 100, we limited our calculations to log T_eff> 3.9. The shape of the resulting theoretical curves for β₁(τ) is very similar to that obtained using grey atmospheres. For a given optical depth and effective temperature, the ATLAS models predict slightly larger GDEs (less than 10%). This is because of the differences in accuracy and completeness of the line opacities used in the calculation of the optical depth in the ATLAS and in the GRANADA code. In general, the calculations of β₁(τ), adopting more sophisticated atmosphere models, support those obtained using the grey approximation.

Fig. 4 Comparison between the present theoretical gravity-darkening exponents (full line) and the observational exponents (see Table 1). The effects of the evolution and metallicity are not taken into account. The dashed line and asterisks represent the calculations of β1 for convective envelopes adopting our original method (Claret 2000, 2004).

The observational data of six fast rotators with an accurate determination of GDE to be compared with the theoretical predictions are listed in Table 1 with the corresponding references. The comparison between the semi-empicial GDEs and the theoretical GDEs is shown in Fig. 4. In this figure each system is identified by its name and the full line denotes the theoretical β₁(τ) (Table 2). A good agreement is achieved for all six systems simultaneously. The region for T_eff ≤ 8000 K needs special attention. Our original method for convective envelopes requires that the distorted structures present the same temperature-pressure slope very deep in the models where the diffusion approximation is valid. More robust calculations of β₁ in this region should combine the two methods (the results using the original method are shown in Fig. 4 as a dashed line) to obtain more consistent results. Therefore, we do not consider the calculations in Fig. 4 as definitive albeit the good agreement with semi-empirical data and we shall explore in more detail the behaviour of the GDE in this critical region in future works.

Although the aspect of Fig. 4 is encouraging, some improvements are still necessary to achieve a more consistent comparison. Some of the these improvements are of theoretical nature, as for example: 1) the inclusion of the meridional current is an important piece to understand the flux distribution of distorted stars. One of our main objectives for the future is to investigate the role of the meridional circulation on the GDE in the frame of our numerical method and to verify whether its inclusion restores the von Zeipel value or not. Important tools for comparison and/or estimation of the effects of the meridional circulation are the rotating evolutionary models generated by the code MESA, which includes the Eddington-Sweet approximation (A. Dotter 2015, priv. comm.). 2) Even though β₁(τ) computed by adopting the grey approach are consistent with those computed with more realistic stellar atmosphere models, for completeness, it would be very helpful to generate grids with smaller intervals in log T_eff and log g to reduce the numerical noise in the atmosphere structure interpolation process. 3) The evolutionary effects, including the input physics of the models, should also be considered for comparison with the semi-empirical data. 4) The effects of differential rotation or the hypothetical presence of a disk should also be incorporated, and

the last one is mainly applicable to the case of close binary stars. 5) Instead of using the Rosseland optical depth to compute β₁(τ), it would be more consistent to adopt the optical depths computed for each specific pass-band, e.g. τ_H and τ_K. We plan to introduce these improvements in future works to evaluate their impact on the theoretical β₁(τ) calculations. On the other hand, concerning observations, it would be very interesting to have more data for fast rotators with higher effective temperatures than those listed in Table 1 and if possible, in the critical region T_eff ≤ 8000 K also to test our theoretical predictions of GDE for fast rotators more widely.

Acknowledgments

I acknowledge Dr. J. Puls and B. Rufino for careful reading of the manuscript, comments, and suggestions. The Spanish MEC (AYA2012-39727-C03-01) is gratefully acknowledged for its support during the development of this work. This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Abstract Service.

References

All Tables

All Figures

Fig. 1 Convective efficiency and gravity-darkening exponent in the HR diagram. The models were computed with the GRANADA code and the masses range from 0.8 up to 3.16 M⊙. Dashed line indicates the locus of the f = 0.5 (Eq. (10)). This line roughly divides the HR diagram in envelopes with high convective efficiency (right) and low convective effciency envelopes (left). The critical log Teff at the ZAMS is around 3.9.

In the text

	Fig. 2 Rosseland optical depth, where the boundary conditions are imposed, as a function of the effective temperature adopted in the calculation of the gravity-darkening exponents (ZAMS models). We note the drop-off for radiative envelopes (log Teff ≈ 3.9).
In the text

	Fig. 3 Theoretical gravity-darkening exponent as a function of the Rosseland optical depth. Model with 4 M⊙, log g = 4.293, and log Teff = 4.165. See text for more details.
In the text

	Fig. 4 Comparison between the present theoretical gravity-darkening exponents (full line) and the observational exponents (see Table 1). The effects of the evolution and metallicity are not taken into account. The dashed line and asterisks represent the calculations of β1 for convective envelopes adopting our original method (Claret 2000, 2004).
In the text