© ESO, 2012

1. Introduction

In close binary stars, both components are distorted by rotation and by tides, with the tides tending to elongate the star, while rotation tends to flatten it at the poles. In a pioneer work, von Zeipel (1924) showed that the flux is not constant on the surface of a distorted star: , where g is the local gravity, T_eff the effective temperature, and β₁ the gravity-darkening exponent (equal to 1.0 for configurations in radiative equilibrium). Later, Lucy (1967) found a mean value 0.32 for β₁ for convective envelopes. Still discussing theoretical investigations, Claret (1998) introduced a numerical method to compute β₁, which embraces both mechanisms of transport of energy and β₁ were presented as continuous functions of the stellar masses, age, chemical composition, and mixing-length parameter.

On the other hand, many researchers working on double-lined eclipsing binaries were able to infer β₁ semi-empirically for several systems (Rafert & Twigg 1980; Pantazis & Niarchos 1998; Niarchos 2000; Djurasevic et al. 2003) to be compared with the theoretical predictions. By using interferometric techniques, some very rapid isolated rotating stars were also investigated, opening a new window on the stellar rotation (van Belle et al. 2001, 2006; Domiciano de Souza et al. 2003; Aufdenberg et al. 2006; Peterson et al. 2006; Monnier et al. 2007; Zhao et al. 2009). In Che et al. (2011) the predicted theoretical transition zone between radiative/convective envelopes (Claret 1998) was confirmed by using rapid rotators, some of them near the break-up velocities.

This amount of data on gravity-darkening permitted a robust comparison with the theoretical values of β₁. A general good inter-agreement between theory and observations was found (see for example, Fig. 3 in Claret 2003). However, some discrepancies are present too, as such Che et al. (2011), who studied the hot star α Leo and found β₁ = 0.75, significantly lower than the value predicted by von Zeipel’s theorem (β₁ = 1.0). On the other hand, the detected scattering of the semi-empirical values of β₁ around the theoretical values of double-lined eclipsing binaries also demand a theoretical explanation. To try to advance in this direction, we introduce in this paper two different methods to derive the theoretical gravity-darkening exponents. In the first one, we use a generalization of von Zeipel’s theorem to obtain an approximated expression for β₁ for slow rotating stars. In the second procedure, we explore the validity of the mentioned theorem in the upper layers of a distorted star by adapting grey and non-grey atmosphere models to our numerical method. We have found deviations in the von Zeipel theorem that could explain the low value of β₁ for α Leo, as well as the mentioned scattering in double-lined eclipsing binaries.

2. An analytical approach to the gravity-darkening exponent

We are interested in calculating the gravity-darkening exponent of a distorted star by using the perturbation theory. To do this, we derive a generalization of von Zeipel’s theorem, first obtained by Kippenhahn (1977). In addition to the generalization, we use it to derive the gravity-darkening exponent. By choosing the integrating factor Λ (to be determined) in such a way that Λg can be derived from a potential, it can be shown, with the equation of hydrostatic equilibrium, that the pressure is constant on equipotential surfaces. For constant mean molecular weight (ideal gas), the temperature will depend on Λ^-1 while the density varies as Λ.

Here we follow Kippenhahn’s (1977) method of deriving the integrating factor Λ. Adopting a general law for the stellar rotation of the form (1)where θ is the angle characterizing the colatitude, ω is the angular velocity and h(r) a generic function. If we assume that the effect of the centrifugal force on the stellar structure is small, including that the gravitational potential at the stellar surface is the Roche potential, the integrating factor for the potential equation is (2)where , χ = r³ω² / (GM) and M is the stellar mass.

On the other hand, the equation that governs the radiative transfer in stellar interior is (3)If we assume that the density and temperature are functions of the potential, that is, ρ = d(Ψ)Λ and T = t(Ψ) / Λ, then Eq. (3) can be written as (4)The above equation is a generalization of von Zeipel’s theorem (Kippenhahn 1977). In fact, if Λ is constant (a conservative centrifugal force), we recover the classical version of the mentioned theorem; i.e., the flux will be proportional to the gravity. With Eqs. (2) and (4) we can evaluate the gravity-darkening exponent β₁. To advance in the derivation of β₁, we calculate the gravities and the fluxes at the pole and at the equator for an equipotential. The potential is given by (5)and the radius of the respective equipotential r can be written as (6)where (7)and P₂(θ,φ) is the second surface harmonic, a the mean radius of the level surface, η₂ is the logarithmic derivative of the spherical harmonic defined through Radau’s equation (η₂ = 3 at the surface for the Roche model), and M_ψ the mass enclosed by an equipotential. By using Eqs. (5) and (6) we can calculate the gravity at the pole and at the equator: (8)Similarly, we can derive the fluxes at the pole and at the equator. At the pole the expression for the flux has a simple form. However, at the equator the calculation of the flux is a little bit more complex. By using Eq. (4) and expanding the value of Λ, we have (in units of ): (9)and therefore (10)Equation (10) is different from the similar one obtained by Kippenhahn (1977), probably due to misprints or to the adopted level of approximation. At the surface, the term is much larger than the variation of ω and α with r. Therefore, the term containing in Eq. (10) can be neglected. We finally obtain the value of β₁ by approximating the logarithmic derivatives of the flux and effective gravity through Eqs. (10) and (8), respectively. By incorporating the effects of the variable opacity we get (11)with and .

Fig. 1 Gravity-darkening exponent for a 4 M⊙ model. Continuous thick line refers to grey calculation (τg = 80), continuous thin (the same but for τg = 5), dashed (the same but for τg = 1). The dashed-dotted line indicates computations by using the ATLAS atmosphere model at τA = 35 and the dotted line, the same but for τA = 1.0. The error bars denote the location of α Leo.

We recover von Zeipel’s value for β₁ (1.0) for uniform rotation and no angle dependence or if κ_ρ = κ_T − 4 and/or N = 2α, in the general case. If we use the semi-empirical value of the gravity-darkening exponent for α Leo by Che et al. (2011) as a constraint, we can estimate the effects of the geometry and differential rotation in that star. The logarithmic derivatives κ_T and κ_ρ can be evaluated from appropriated tables of opacities, such as those by Ferguson et al. (2005). For a typical chemical composition (X = 0.70, Z = 0.02) and for log T = 4.0, log ρ =  −6.0 → κ_T = 9.20 and κ_ρ = 0.28. Considering no dependence on the angle θ we obtain α =  −0.1. Of course, such value will depend on the depth of the stellar envelope we are dealing with. We should only take such an estimation as an approximation. Equation (11) was derived by using a perturbation theory (slow rotation), but since α Leo is rotating near the break-up velocity, it is not sufficient to analyse such stars. However, this equation can be applied in those distorted binary systems or slow isolated rotators for which the rotation rates are far from the critical ones. A more adequate approach to β₁ for rapid rotators is analysed in the next section.

3. Deviations from von Zeipel’s theorem at the upper layers; final remarks

Our numerical method which was introduced to derive β₁ is based on the triangle strategy used in stellar evolution calculations (Claret 1998). This method does not depend on the details of the rotation laws. It is based fundamentally on the star’s distortions simulated by different temperature distribution corresponding to different log g’s and log T_eff’s. We now comment on the main characteristics of this method. As is usual in stellar evolution calculations, the outer-layer integration should be carried out after each time step. If the external boundary conditions of the interior are unchanged at the fitting point, the envelope integrations would be the same as before. Such calculations can be avoided thus saving computational time. In practice, three envelope computations are performed corresponding to three points in the HR Diagram. If the next point on the evolutionary track is still within this triangle, the same envelope integrations can be used without loss of accuracy (provided that the triangle is sufficiently small). Of course, new integrations are necessary if this condition is not fulfilled.

To simulate a distorted stellar configuration we increased the number of points in the triangle strategy. Given a point on an evolutionary track (characterized by log T_eff and log g), a series of envelopes with log g′ = log g + Δlog g and trial values of effective temperatures are computed. Often, these new envelopes do not present the same physical conditions – pressure x temperature relationship – at the bottom as the original envelope. By interpolation in these envelopes, another envelope is found that presents the same physical conditions as the reference point (such an envelope has an effective temperature , which is different from the original). By repeating the process n times (by varying the values of Δlog g), a family of n interpolated envelopes (characterized by log g′ (n) and log (n) is found with the same physical conditions at the bottom. Finally, by numerically differentiating this family of envelopes, we obtain β₁ for each point on a given evolutionary track (). However, in calculating the corresponding envelopes, the equation of diffusion was used throughout. On the other hand, von Zeipel’s theorem was also derived by using this approximation. Intuitively, diffusion is only an appropriate approach for the deeper layers. In fact, the moment of the transfer equation can be written as (12)where F_ν(τ) is flux, B_ν the Planck’s function and τ the optical depth scale. At great depth (there are too many photon mean free paths from the surface), the second term on the right-hand side of Eq. (12) vanishes and the transfer equation reduces to a simple form, i.e., the diffusion approximation. As a consequence, Eq. (12) reduces to Eq. (3), which is the basis of the derivation of the classical von Zeipel’s theorem (β₁ = 1.0). This result is used as a control of our method at large optical depths as well as the solution for τ = 0.

It is interesting to explore the behaviour of β₁ at the atmosphere of a distorted star where the diffusion equation is no longer valid. Fortunately, our numerical method can be extended to the surface layers, and we are able to evaluate the influence of the optical depth on β₁. For each stellar model generated with our evolution stellar code (Claret 2004), characterized by a given mass, log T_eff, log g, and metallicity, a grey atmosphere is computed. By transferring our numerical method to these upper layers we can derive the behaviour of β₁ as a function of the effective temperature and on the optical depth. We restrict our calculations here to hotter models because we are interested in exploring the range of validity of von Zeipel’s theorem. The maximum value of Δlog g adopted was 0.25 (corresponding to the difference between the log g at the pole and at the equator of α Leonis). Figure 1 shows the results of these calculations for three optical depths using an evolutionary model with 4 M_⊙. As expected, at great optical depth, we restore the classical value of β₁ (1.0). However, when we consider more external layers, the corresponding values of β₁ decrease and there are also shifts in the transition zone (radiative/convective). This indicates that von Zeipel’s theorem is not strictly valid at the upper layers. The values of β₁ shown in Fig. 1 for colder models (based on grey atmospheres) are not fully reliable.

To check the reliability of the results described in the last paragraph by adopting more elaborated atmosphere models, we also use our numerical method but with a flag to use the structure of the ATLAS stellar atmosphere models, instead of the grey models. The ATLAS 9 (Kurucz, priv. comm.) models make use of a much more complex transfer equation, and the resulting calculations are more reliable than those by adopting the grey approximation. As in the grey case, the values of β₁ are systematically lower than 1.0, except at great depth, where the von Zeipel value is restored for hotter models. There are some differences if we compare β₁(ATLAS) with the β₁s derived by adopting grey models, but they are not as large for a given optical depth. We also note the numerical noise for ATLAS that is caused by the interpolation process in the atmosphere structure. (The neighbouring grid points in log g and effective temperatures differ 0.5 and up to 2500 K, respectively.)

As mentioned in the Introduction, the semi-empirical values of β₁ obtained from double-lined eclipsing binaries and rapid rotating stars compare well with the theoretical predictions based on our original method (see for example, Pantazis & Niarchos 1998; Niarchos 2000; Djurasevic et al. 2003; Che et al. 2011). There is some scattering around the theoretical predictions based on our numerical method using stellar envelopes, based on the diffusion equation. Within this context, we may interpret the detected deviations from von Zeipel’s theorem at lower optical depths as theoretical error bars, i.e., as lower limits for β₁. The recent derivation of the semi-empirical β₁ for the hot rapidly rotating star α Leo (Che et al. 2011), which is between 0.64 and 0.80, could be explained if we consider the deviations found here (see Fig. 1). One question arising from this discussion is where we should put such a lower limit, say, what is the optical depth we should select to compute the deviations of β₁? A measure of the drift of the radiation field can be evaluated by taking the ratio between the anisotropic/isotropic terms, which is given approximately by  [T_eff / T(τ)] ⁴. The deviations from the diffusion equation become important for τ_c of the order of 1. We expect that the mentioned critical optical depth depends on spectral type. We adopt τ_c = 1 provisionally, although it may also depend on the adopted stellar atmosphere models and perhaps on more accurate observational constraints. The mentioned scattering for close binaries could be explained with the help of Eq. (11), including presence of disks, or a combination of that equation and the effects described in this section.

The progress in deriving semi-empirical values of β₁ is notorious for isolated rapid rotating stars, including those near

the break-up velocities. On the other hand, the researchers in eclipsing binaries are continually bringing new data to be compared with theoretical predictions. That conjunction surely will contribute to refine a little more the theoretical calculation of the gravity-darkening exponents. Among the improvements that can be implemented, we plan to introduce a dependence of β₁ on the details of the rotation laws, including the degree of the distortions.

Acknowledgments

This paper is dedicated to the memory of Jens Viggo Clausen. I would like to thank V. Costa and an anonymous referee for useful suggestions. The Spanish MEC (AYA2006-06375, AYA2009-14000-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 Figures

	Fig. 1 Gravity-darkening exponent for a 4 M⊙ model. Continuous thick line refers to grey calculation (τg = 80), continuous thin (the same but for τg = 5), dashed (the same but for τg = 1). The dashed-dotted line indicates computations by using the ATLAS atmosphere model at τA = 35 and the dotted line, the same but for τA = 1.0. The error bars denote the location of α Leo.
In the text