© ESO 2021

1. Introduction

Gravity-darkening is an important piece of the stellar structure and evolution that has been studied for almost a hundred years (von Zeipel 1924). Gravity-darkening exponents (GDE) are key tools for analysing light curves of eclipsing binaries or in isolated rotating stars through long-baseline optical and infrared interferometry. Considering stars in strict radiative equilibrium (pseudo-barotrope), in 1924 von Zeipel showed that the variation of brightness over the surface is proportional to the effective gravity, that is,

$$\begin{array}{c}\hfill F=-\frac{4ac{T}^3}{3\kappa \rho }\frac{\mathrm{d}T}{\mathrm{d}\psi }{g}^{{\beta }_1},\end{array}$$(1)

or equivalently

$$\begin{array}{c}\hfill {T_{\mathrm{eff}}}^4\propto g^{{\beta }_1},\end{array}$$(2)

where ψ is the potential, g the local gravity, T the local temperature, κ the opacity, ρ the local density, a the radiation pressure constant, c the velocity of light in vacuum, T_eff the effective temperature, and β₁ = 1.0 is the GDE, which is a bolometric quantity. Although we are still far from fully understanding the gravity-darkening phenomenon, in the last 22 years several theoretical papers related to the GDE have appeared in the literature that have shed some light on the scenario, such as Claret (1998, 2000, 2012), Espinosa Lara & Rieutord (2011, 2012), and McGill et al. (2013). For a historical summary of the theoretical research on the GDE see Claret (1998, 2016). From now on we also designate the GDE as β₁.

Another important and complementary ingredient in the analysis of systems with non-spherical configurations are the so-called gravity-darkening coefficients (GDCs). As we can only observe determined band-limited stellar flux (not the bolometric), the introduction of these coefficients is necessary in order to model a distorted star. Such a coefficient can be written as (Claret & Bloemen 2011):

$$\begin{array}{cc}& y(\lambda ,T_{\mathrm{eff}},log[\mathrm{A}/\mathrm{H}],logg,V_{\xi })\hfill \\ \hfill & =\left(\frac{\mathrm{d}ln{T}_{\mathrm{eff}}}{\mathrm{d}lng}\right){\left(\frac{\partial ln{I}_o(\lambda )}{\partial ln{T}_{\mathrm{eff}}}\right)}_g+{\left(\frac{\partial ln{I}_o(\lambda )}{\partial lng}\right)}_{T_{\mathrm{eff}}},\hfill \end{array}$$(3)

where λ is the wavelength, I_o(λ) the intensity at a given wavelength at the centre of the stellar disc, and V_ξ is the microturbulent velocity. We note that the expression $\left(\frac{\mathrm{d}ln{T}_{\mathrm{eff}}}{\mathrm{d}lng}\right)$ can be written as β₁/4. In order to progress in our investigation of stellar configurations distorted by rotation and/or tides, we recently studied the GDCs for the case of compact stars (Claret et al. 2020). In that study we computed GDC for DA, DB, and DBA white dwarf models, covering the transmission curves of the Sloan, UBVRI, Kepler, TESS, Gaia, and HiPERCAM photometric systems. These computations are available for log [H/He] = −10.0 (DB), −2.0 (DBA), and He/H = 0 (DA). The covered range of log g was 5.0–9.5, while for the effective temperatures the respective range was 3750 K–100 000 K.

From an observational point of view, discoveries of binary systems whose components are tidally and/or rotationally distorted white dwarfs are ongoing (see e.g. Burdge et al. 2019a,b, 2020; Kupfer et al. 2020 and references therein). However, as far as we know, there are no specific GDE calculations for white dwarf sequences or even for an individual model. In this short paper we present, for the first time, GDEs for three white dwarf cooling sequences. We also introduce some improvements to our methods for calculating the GDE: β₁ is computed as a function of the optical depth, that is, β₁ = β₁(τ). In addition, we introduce an extra condition in our calculations, namely the relationship between the local gravities and the corresponding optical depths for a given equipotential surface. In the following, sections we discuss these points in more detail.

The paper is organised as follows: Sect. 2 is dedicated to describing our methods and some aspects of the distorted stellar configurations. In Sect. 3 we discuss both observational and theoretical evidence for deviations from von von Zeipel’s theorem (1924) and present our results and conclusions.

2. The numerical method

Our numerical method is based on the triangles strategy introduced by Kippenhahn et al. (1967). A complete description of our method can be found in Claret (1998), but for the sake of clarity, we summarise it below. To save computing time, Kippenhahn et al. (1967) introduced an ingenious method: when an evolutionary sequence is being calculated, if the external boundary conditions are unchanged at the fitting point (envelope-interior), the outer layer integrations must be the same as the previous ones. Three envelopes corresponding to three points in the Hertzsprung-Russell (HR) diagram around the current values of luminosity and effective temperature are computed. As the model evolves, its properties are checked to see if they remain within this triangle. If so, the boundary conditions are unchanged. It is important to highlight that this strategy is valid only if the triangle is sufficiently small. This warning is particularly valid for convective envelopes. To simulate a distorted star, we are interested in several triangles showing different effective temperature distributions over the surface. Therefore, we increase the number of triangles in the HR diagram, that is, we compute several envelopes with different temperature distributions but imposing the same physical conditions at a given interior point. Figure 2 in Claret (2000) shows a simplified scheme of our method where only three points are shown in the HR diagram. Once the triangles have been computed for each point of the evolutionary track we can derive β₁ by differentiating the neighbouring envelopes. Such a procedure is performed for the next evolutionary track point and so on.

Our method has some advantages: (a) it can be applied to convective and/or radiative envelopes; (b) one can investigate the influence of the optical depth in the GDE by changing the fitting point to impose the boundary conditions, without loss of generality; (c) the GDE can be computed as a function of initial mass, chemical composition, evolutionary stage (in the present paper up to the white dwarf phase), and other ingredients of the input physics, and (d) more realistic atmosphere models can easily be incorporated as external boundary conditions, as done in Claret (2012).

We introduce an extra condition in our procedure: the relationship between local gravity and optical depth over a given equipotential surface. In reality, solving the hydrostatic equilibrium equation for two points on an equipotential characterised by [g(μ),τ(μ, ψ)], we obtain g(μ)τ(μ, ψ) = g(μ_o)τ(μ_o, ψ). In this relationship, g is the local gravity, ψ is the total potential (rotational one included), τ is the optical depth, and μ is given by cos(ϕ) where ϕ is the angle between the radiation field and the z axis. The subscript o indicates the reference point, for example the pole. To guarantee that the triangle technique represents an equipotential, we have taken as reference the absolute dimensions for each point of the evolutionary models, that is, [g(μ_o),τ(μ_o, ψ)]. The equipotentials are then configurated introducing additional triangles from this point.

As outlined in Sect. 1, for stars in strict radiative equilibrium, von Zeipel’s theorem predicts an exponent β = 1.0. However, if we inspect a typical HR diagram in the log g × T_eff plane for stars with different initial masses, for example, 1.0 M_⊙ (convection predominates in the envelope) and 10.0 M_⊙ (envelope predominantly in radiative equilibrium), it can be verified that both models have different average slopes. If we make a simple analogy between these two slopes and the GDE given by Eq. (2), the respective β₁ would be different, with the one corresponding to the star with 1.0 M_⊙ being smaller than that of the model with 10.0 M_⊙; see Fig. 3 by Claret (1998) for a graphic example. This figure seems to indicate that stars with convective envelopes do not strictly obey von Zeipel’s theorem. Additional and more elaborate evidence has also been found to support these statements, such as that outlined in Claret (2012, 2016) for example. In these latter studies, significant deviations were found when the GDEs are computed for the upper layers of a distorted star.

On the other hand, Kopal (1959) derived the following equation for the stellar distortions:

$$\begin{array}{c}\hfill \frac{g-g_o}{g_0}={\displaystyle \sum _j(1-\frac{5}{{\mathrm{\Delta }}_j})(\frac{r}{a_1}-1),}\end{array}$$(4)

where g_o is the reference local gravity, Δ_j = 1 + 2k_j, and k_j is the apsidal motion constant of order j. Therefore, the radius of an equipotential r (order 2) can be written as

$$\begin{array}{c}\hfill r=a_1(1-f_2{P}_2(\theta ,\varphi )),\end{array}$$(5)

with

$$\begin{array}{c}\hfill f_2=\frac{5{\omega }^2{a}_1^3}{3G{M}_{\psi }(2+{\eta }_2)},\end{array}$$(6)

where ω is the angular velocity, P₂(θ, ϕ) is the second surface harmonic, a₁ is the mean radius of the level surface, η₂ is the logarithmic derivative of the amplitudes of the surface distortions defined through Radau’s differential equation, and M_ψ is the mass enclosed by an equipotential. Claret (2000) has shown that there is a close connection between the GDE and the shape of the distorted stellar configuration, its internal structure, and the details of the rotation law. In fact, for stellar masses around 1.5 M_⊙ there is a change in the predominant source of thermonuclear energy from the proton–proton chain to the CNO cycle. This change causes a readjustment of the mass concentration through the parameter η₂. We reiterate the fact that η₂ is connected to k₂ through a simple equation: k₂ = (3 − η₂(R))/(2(2 + η₂(R))), where η₂(R) is evaluated at the stellar surface. This readjustment will affect how a star reacts under distortions and will consequently effect the parameter β₁ (see Eq. (3) and also Fig. 1 by Claret 2000 where log k₂ and β₁ are shown for ZAMS models with masses varying from 0.075 to 40.0 M_⊙). In addition, convection also begins to contribute significantly to the total flux in this range of effective temperatures. For masses greater than 1.5 M_⊙ the mass concentration decreases almost linearly with the stellar mass.

3. Discussion of the results and final remarks

Some years ago, Claret (2016, see the corresponding Eqs. (8) and (9)) adopted an expression for the convective efficiency (ratio of the convective to the total flux, denoted by the symbol f) to investigate the GDEs. Here we generalise that result for a range of opacities (through parametrised formulae) as indicated by Eq. (7):

$$\begin{array}{c}\hfill f\approx A\gamma {\left(\frac{r}{R}\right)}^2{\left(\frac{3{\mathrm{\Gamma }}_1}{5{\mu }_1\beta }\right)}^{1/2}\left(\frac{2c_{\mathrm{P}}{\mu }_1\beta }{5}\right)T^{1/2}\left[\frac{g}{T_{\mathrm{eff}}^{(4n+4+|n+s|)}}\right].\end{array}$$(7)

In the above equation, A is a constant, R the star radius, r the radial coordinate, μ₁ the mean molecular weight, β the ratio of gas to total pressure, c_P the specific heat at constant pressure, and Γ₁ is the adiabatic exponent. The parameter γ is given by

$$\begin{array}{c}\hfill \gamma =\frac{1}{2}\frac{\overline{v}}{v_{\mathrm{s}}}\frac{\mathrm{\Delta }T}{T},\end{array}$$(8)

with $\overline{\mathit{v}}$ being the mean convective velocity, ΔT the excess of temperature of a rising element over the mean temperature of the surrounding environment, and v_s is the velocity of sound. Here, we assume that for the derivation of Eq. (7) the opacity can be written as κ ≈ κ′PⁿT^−s or κ ≈ κ″ρⁿT^{−n − s}, where κ′ and κ″ are constants, and n > 0 and s < 0, assuming a perfect gas equation. Although we will not use this expression directly because it is only an approximation, it is useful in order to analyse the correlation between f – through the role of the opacities – and β₁, as we see in the following paragraphs.

On the other hand, the evolutionary track from the pre-main sequence (PMS) up to the white dwarf stage was computed using the MESA module (Paxton et al. 2011, 2013, 2015), version 7385. The basic input physics of the MESA code is given in the above references. Another set of DA- and DB-type white dwarf models was computed with LPCODE (Althaus et al. 2001a,b). The models were followed from the zero age main sequence (ZAMS) up to the white dwarf stage (Renedo et al. 2010). This code considers modern input physics such as a detailed network for thermonuclear reactions, OPAL radiative opacities, full-spectrum turbulence theory of convection, detailed equations of state, and neutrino emission rates.

We first discuss the cases of ‘pure’ DA- and DB-type white dwarfs, without considering the previous evolutionary phases. We analyse models with 1.0 M_⊙ although the results are similar if we consider models with different initial masses. Figure 1 illustrates the resulting GDE computed at optical depth τ = 100.0 (thick continuous line) using the modified triangles strategy method. The log g during the evolution of this cooling sequence vary approximately between 8.3 and 8.7. As expected, for high effective temperatures where radiative transfer predominates and at large optical depth, the results are compatible with the equation of diffusion, that is, with those resulting from the von Zeipel (1924) approach. However, as the model evolves, the influence of convection begins to appear which translates into a drop of β₁. The drop-off threshold is located at log T_eff ≈ 4.12 and we have found large deviations from the von Zeipel theorem for effective temperatures smaller than 10 000 K. We note that this transition temperature is slightly different from that typical of main sequence stars (log T_eff ≈ 3.90). This is due to differences in chemical compositions and degree of compactness (log g). The dashed (log g = 8.5) and thin continuous (log g = 8.0) lines are useful to establish a more direct connection between the contribution of the convection to the flux and β₁. Such a contribution can be set in terms of the ratio F_c(τ)/F_r(τ) where F_c(τ) is the convective flux and F_r(τ) is the radiative one for a given optical depth (these data were kindly provided by Cukanovaite in prep.; see also Cukanovaite et al. 2019). The F_c(τ)/F_r(τ) ratio is slightly different from that given by Eq. (7) but shows the same appearance, taking into account the different scales: the F_c(τ)/F_T(τ) ratio varies from 0 to 1.0, F_T(τ) being the total flux.

Fig. 1. DA-type white dwarf models (1.0 M⊙). The thick solid line represents the GDE as a function of effective temperature. The continuous thin line indicates the ratio Fc(τ)/Fr(τ) for log g = 8.5 and the dashed one denotes the same but for log g = 8.0. All calculations were performed at τ = 100.0. We note that the Teff range of the atmosphere models by Cukanovaite (in prep.) is limited to Teff ≤ 60 000 K. The arrow indicates the direction of time evolution.

We can use Eq. (7) to help us visualise the behaviour of β₁ and its relation to opacity. The connection between β₁ and the convective contribution to the flux is evident; it is especially clear near the maxima and minima. The behaviour of the GDE shown in Fig. 1 (also that of the following figures) is the result of the combination of the different sources of opacity but we can gain some insight into its influence on the GDE by adopting some opacity laws in their parameterised form. Firstly, we analyse the effects of opacity in the colder models. For this range of effective temperatures one of the predominant sources of opacity is the negative hydrogen ion whose dependence on ρ and T is given by κ ≈ κ₁ρ^1/2T^7.7, where κ₁ is a constant. We note the high dependence of negative hydrogen ion opacity with temperature. The expression $\frac{g}{T_{\mathrm{eff}}^{(4n+4+|n+s|)}}$ in Eq. (7) drives the behaviour of the gravity-darkening exponent, and, introducing the relevant values of n and s for such models at a given convective efficiency, we get β₁ ≈ 0.30. In addition, the opacities ff, bf, and bb, because of electronic transitions, are given by Kramers law: κ ≈ κ₂ρT^−7/2, where κ₂ is a constant. The approximate corresponding value of β₁ in this case is ≈0.35. Considering that Eq. (7) is just a rough approximation, these results are surprising because they predict deviations from the von Zeipel’s theorem and, in addition, they are in reasonable agreement with the values of β₁ found in the literature for late-type stars (semi-empirical and theoretical ones; see below). On the other hand, one of the main sources of continuum opacity in hot star atmospheres is the so-called Thomson scattering. As known, this opacity is ‘grey’ because there is no dependence on temperature and density and it can be written as follows for the case of complete ionisation: κ ≈ 0.2(1.0 + X), where X is the hydrogen content. Using the same treatment as for the negative hydrogen ion case and Kramers law, we obtain β₁ ≈ 1.00, which is in good agreement with the predictions of the von Zeipel’s theorem for hot stars. Again it is gratifying that Eq. (7) is capable of predicting the typical values of β₁, even considering its limitations.

In addition to the analysis of the influence of convective flux on β₁, there are other ways to correlate β₁ with some physical magnitudes. Considering the additive properties of the specific entropy in the nonrelativistic case we have

$$\begin{array}{c}\hfill s=B+\frac{N_{o}k}{{\mu }_i}ln\frac{T^{3/2}}{\rho }+\frac{N_{o}k}{{\mu }_{\mathrm{e}}}[\frac{5}{3}\frac{F_{3/2}(\alpha )}{F_{1/2}(\alpha )}+\alpha ]+\frac{4a}{3}\frac{T^3}{\rho },\end{array}$$(9)

where the symbols have the following meaning: B is a constant, α the degeneracy parameter, N_o is Avogadro’s number, k the Boltzmann constant, μ_i the mean molecular weight per ion, and μ_e is the mean molecular weight per electron. The functions F_1/2(α) and F_3/2(α) are auxiliary functions that can be written as

$$\begin{array}{c}\hfill F_{1/2}(\alpha )={\displaystyle {\int }_0^{\infty }\frac{u^{1/2}\mathrm{d}u}{e^{\alpha +u}+1},}\end{array}$$(10)

and

$$\begin{array}{c}\hfill F_{3/2}(\alpha )={\displaystyle {\int }_0^{\infty }\frac{u^{3/2}\mathrm{d}u}{e^{\alpha +u}+1},}\end{array}$$(11)

where u = p²/(2mkT), with p being the particle momentum and m its mass. The three components in Eq. (10) above can be easily identified: the second term is the entropy due to ions, the third is connected to electrons, and the fourth to radiation. We can see in Fig. 2 how the GDE is related to the entropy for the same models and conditions shown in Fig. 1. As mentioned, the convection onset is in log T_eff ≈ 4.12 (log T_{eff onset}) where a sudden variation of the entropy is indicated by a vertical arrow. This effective temperature is in agreement with those by Tremblay (in prep.). There are three other inflexion points (also marked with vertical arrows). These four points shown in Fig. 2 indicate that the entropy varies with the effective temperature (also with log g), which in turn drives the behaviour of the GDE. The changes in β₁ with entropy are not surprising. Indeed, the differential of the entropy is given by dS = c_p(∇ − ∇_adia), where ∇ = dlnT/dlnP and ∇_adia = (dlnT/dlnP)_adia. On the other hand, it can be shown that, alternatively, the convective flux F_c(τ) ∝ (∇ − ∇_adia)^x, where x = 3 (small convective efficiency) or x = 3/2 for large convective efficiency. Thus, under determined physical conditions, the entropy can also be used as a convective stability criterion.

Fig. 2. Entropy for the same models and conditions shown in Fig. 1 (red lines). The vertical arrows indicate the points where there is a marked variation of the entropy with effective temperature. The solid line represents the entropy for τ = 100 while the dashed-dotted line denotes τ = 500; both for log g = 8.5. The horizontal arrow indicates the direction of time evolution.

Another point to note is that the entropy does not significantly depend on the optical depth for effective temperatures ≤T_{eff onset}. Indeed, the curves for τ = 100.0 and 500.0 coincide for T_eff ≤ T_{eff onset}. This implies that the behaviour of the resulting GDEs should not vary significantly, at least within the range of optical depths, log g, and effective temperatures explored here. We note that, for hot models, because the entropy hardly varies with T_eff, the corresponding GDE values are almost constant and equal to 1.0, which restores von Zeipel’s theorem. These results are more general than those obtained decades ago, whose GDE value was constant and approximately equal to 0.32 for envelopes in convective equilibrium.

Figure 3 shows a comparison between the GDE for DB and DA models at τ = 100.0. The profiles of β₁ are very similar, with the exception of the transition zone where there is a shift due to the difference in the chemical composition of the models and its influence on pressure and temperature.

Fig. 3. Comparison between the GDE for a DB model (solid line, 1.0 M⊙) and DA (dashed line, 0.5 M⊙). All calculations were performed at τ = 100.0. The horizontal arrow indicates the direction of time evolution.

An interesting comparison that can be made is related to the masses of white dwarfs. In Fig. 4 we show the evolution of the GDE for DA models with masses 1 M_⊙ (continuous line) and 0.52 M_⊙ (dashed line). The calculations were also performed for τ = 100.0. Because of the dependence of the onset of convection on the local gravity, the GDEs for the model with 0.52 M_⊙ are shifted towards lower temperatures, in reasonable agreement with other studies (Cunningham et al. 2019, Tremblay, in prep.).

Fig. 4. Effect of log g on the onset of convection for DA models. The continuous line represents a sequence of DA model (1.0 M⊙) while the dashed one also denotes DA models but with 0.52 M⊙. The horizontal arrow indicates the direction of time evolution.

As indicated, another set of evolutionary models was generated with the MESA module with an initial mass of 2.0 M_⊙, X = 0.703, and Y = 0.277. For stars with convective envelopes, we employed the standard mixing-length formalism (Böhm-Vitense 1958) with α_MLT = 1.80. The opacity tables adopted are those given by Grevesse & Sauval (1998). The models were followed from the PMS up to the cooling white dwarf stage. Figure 5 shows the complete evolutionary track in the HR diagram. The final mass in the cooling stage is 0.56 M_⊙. The results concerning gravity-darkening are shown in Fig. 6 where we add the profile of β₁ shown in Fig. 1 for comparison. Within our present level of approximation, the differences are small, being more appreciable only in the interval around log T_eff ≈ 3.90.

Fig. 5. HR diagram for an initial mass of 2.0 M⊙ from the PMS to cooling white dwarf stage. Initial chemical composition X = 0.703, Y = 0.277, αMLT = 1.80.

As we can see in Figs. 1–4 and 6, there are clear indications of deviations from the approach by von Zeipel (1924). Deviations from von Zeipel’s theorem were previously found using suitable evolutionary models in stars evolving towards the branch of red giants and/or in low-mass main sequence stars (low effective temperatures) where the flux is predominantly convective in their envelopes (Claret 1998, 2000). Additional evidence for deviations from von Zeipel’s theorem comes from 3D simulations of cold stars (Ludwig et al. 1999). A more complete historical review about this subject can be found in above references.

Fig. 6. GDE for the white dwarf cooling sequence for the models shown Fig. 5 (continuous line). The DA model shown in Fig. 1 is represented by a dashed line. The horizontal arrow indicates the direction of time evolution.

On the other hand, there is an analytical approach by which we can explore the behaviour of GDEs in compact stars whose outermost layers are in radiative equilibrium. Such an analysis was carried out in Claret (2012) for main sequence stars. Here we adapt the physical conditions for the case of white dwarfs. Such an equation can be written as (see Appendix A)

$$\begin{array}{c}\hfill {\beta }_1\approx 1+[\frac{2+N-\chi (4+N)}{2+N}-\frac{2+N+\chi (-4+8{\alpha }_1-5N)}{2+N+\chi ({\kappa }_{\rho }-{\kappa }_T)(N-2{\alpha }_1)}]\frac{{\tau }_{\mathrm{p}}}{{\tau }_{\mathrm{e}}}.\end{array}$$(12)

Equation (12) opens up some possibilities on the investigation of the distribution of rotational velocities through the parameter α₁ and on geometry through N. We note that this equation is valid only for hot stars. We note that the ratio τ_p/τ_e ≤ 1.0. For example, we would restore the predictions of von Zeipel’s theorem for the cases α₁ = N/2 or for uniform rotation and no θ dependence. An interesting feature of Eq. (12) can be explored and is related to κ_ρ and κ_T. Depending on the behaviour of these two variables, we could have values of β₁ significantly smaller than 1.0¹. For example, for α = −1.0, N = 1, χ = 0.2, τ_p/τ_e = 0.71, and β₁ = 0.80, and the resulting condition would be

$$\begin{array}{c}\hfill -6.0⪅[{\left(\frac{\partial ln\kappa }{\partial ln\rho }\right)}_T-{\left(\frac{\partial ln\kappa }{\partial lnT}\right)}_{\rho }]⪅-5.0.\end{array}$$(13)

Such a condition is approximately fulfilled for envelopes with effective temperatures in the interval 6000 K ⪅ T_eff ⪅ 13 000 K. There are no semi-empirical data yet for white dwarfs in this effective temperature range for comparison, but it is interesting to note that values of β₁ smaller than 1.0 were detected using long-baseline optical/infrared interferometry in isolated fast rotators. Some of these systems show effective temperatures within the mentioned range. A summary of their properties can be found in Table 1 of Claret (2016).

Values of the GDE smaller than 1.0 were also observationally detected in main sequence and/or subgiants stars in binary systems. However, a direct connection between the results from long-baseline interferometry and from eclipsing binaries (mostly main sequence stars) and those provided by Eq. (12) (compact ones) is not straightforward. However, it can give us some clue for future research. A turning point concerning semi-empirical β₁ was the pioneering paper by Rafert & Twigg (1980) using eclipsing binaries. Such research was followed by others, such as Pantazis & Niarchos (1998), Niarchos (2000), and Djurasevic et al. (2003, 2006), for example, who explored a wide range of effective temperatures. The semi-empirical GDEs derived by these latter authors for hot stars are scattered around the classical von Zeipel value. Some of these systems show β₁ as low as 0.60. Another important result of these investigations is that the derived values for systems with cooler components also contradict the predictions of von Zeipel’s theorem. These semi-empirical GDE values are quite significant but do not yet constitute a critical test of the theory of temperature distribution on a distorted stellar surface. However, despite the fact that these observations are not conclusive, they seem to indicate a transition zone for the GDE. That zone approximately coincides with the change in the prevalence of the process of energy transport in the envelopes, that is, radiative to convective, as indicated by the Fig. 3 in Claret (2003). Such a transition zone is also predicted in the present paper for compact stars.

Finally, it would be very interesting and useful if observers were to focus their attention on close binary systems constituted by white dwarfs distorted by rotation and tides, so that the validity of the preliminary calculations presented here can be verified. To constrain the GDE values observationally, it would be necessary to investigate eclipsing binary white dwarfs that are double-lined and bright enough to obtain good radial-velocity semi-amplitudes for both components. We hope that such systems will be found in the not too distant future.

Acknowledgments

I would like to thank E. Cukanovaite for providing models of the structure of white dwarfs atmospheres and an anonymous referee for his/her helpful suggestions. The Spanish MEC (ESP2017-87676-C5-2-R, PID2019-107061GB-C64, and PID2019-109522GB-C52) is gratefully acknowledged for its support during the development of this work. A.C. also acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). 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

Appendix A: Brief description of the derivation of Eq. (12)

For completeness, we give here a brief summary of the derivation of Eq. (12) which is a generalisation of von Zeipel’s theorem. For more details, see Kippenhahn (1977) and Claret (2012).

Let us assume that the angular velocity has the following general form

$$\begin{array}{c}\hfill {\omega }^2=h(r){sin}^N\theta ,\end{array}$$(A.1)

where h(r) is a generic function and θ is the colatitude angle. In our approach we assume that the effect of the centrifugal force on the stellar structure is small. The total potential ψ can be written as

$$\begin{array}{c}\hfill \psi =\frac{-GM}{r}-\frac{r^2{\omega }^2{sin}^2\theta }{N+2}.\end{array}$$(A.2)

The integrating factor was chosen in such a way that Λg can be derived from a potential. As a consequence, it can be shown that the pressure is constant on equipotential surfaces. If the effects of rotation are small, we obtain for the integrating factor:

$$\begin{array}{c}\hfill \mathrm{\Lambda }=1+\frac{1}{N+2}\chi (N-2{\alpha }_1){sin}^2\theta .\end{array}$$(A.3)

In the above equation we define ${\alpha }_1=\frac{\partial ln\omega }{\partial lnr}$ and χ = r³ω²/(GM_ψ).

On the other hand, the equation of the radiative transport can be written as

$$\begin{array}{c}\hfill \mathit{F}=\frac{4ac}{3\kappa \rho }{T}^3\mathbf{\nabla }T.\end{array}$$(A.4)

Assuming that the density and temperature are functions of the structural form of the potential we have ρ = d(Ψ)Λ and T = t(Ψ)/Λ and

$$\begin{array}{c}\hfill \mathit{F}=\frac{4ac{t}^3}{3\kappa d}(\frac{\mathrm{d}t}{\mathrm{d}\mathrm{\Psi }}\frac{\mathit{g}}{{\mathrm{\Lambda }}^4}+\frac{t}{{\mathrm{\Lambda }}^6}\mathbf{\nabla }\mathrm{\Lambda }).\end{array}$$(A.5)

As an extra condition we adopt the relationship g(μ_e)τ(μ_e, ψ) = g(μ_p)τ(μ_p, ψ), where the subscripts e and p denote values at the equator and at the poles. Calculating the values of local gravity, of Λ and of Eq. (A.5) at the pole and at the equator we finally obtain the GDE through the logarithmic derivatives of the flux and effective gravity

$$\begin{array}{c}\hfill {\beta }_1\approx \frac{\left(\frac{F_{\mathrm{e}}-F_{\mathrm{p}}}{F_{\mathrm{p}}}\right)}{\left(\frac{g_{\mathrm{e}}-g_{\mathrm{p}}}{g_{\mathrm{p}}}\right)}.\end{array}$$(A.6)

After some algebra and incorporating the effects of variable opacity we finally get

$$\begin{array}{cc}& {\beta }_1\approx 1+[\frac{2+N-\chi (4+N)}{2+N}-\frac{2+N+\chi (-4+8{\alpha }_1-5N)}{2+N+\chi ({\kappa }_{\rho }-{\kappa }_T)(N-2{\alpha }_1)}]\hfill \\ \hfill & \frac{{\tau }_{\mathrm{p}}}{{\tau }_{\mathrm{e}}},\hfill \end{array}$$(A.7)

where ${\kappa }_{\rho }\equiv {\left(\frac{\partial ln\kappa }{\partial ln\rho }\right)}_T$ and ${\kappa }_T\equiv {\left(\frac{\partial ln\kappa }{\partial lnT}\right)}_{\rho }$.

All Figures

Fig. 1. DA-type white dwarf models (1.0 M⊙). The thick solid line represents the GDE as a function of effective temperature. The continuous thin line indicates the ratio Fc(τ)/Fr(τ) for log g = 8.5 and the dashed one denotes the same but for log g = 8.0. All calculations were performed at τ = 100.0. We note that the Teff range of the atmosphere models by Cukanovaite (in prep.) is limited to Teff ≤ 60 000 K. The arrow indicates the direction of time evolution.

In the text

	Fig. 2. Entropy for the same models and conditions shown in Fig. 1 (red lines). The vertical arrows indicate the points where there is a marked variation of the entropy with effective temperature. The solid line represents the entropy for τ = 100 while the dashed-dotted line denotes τ = 500; both for log g = 8.5. The horizontal arrow indicates the direction of time evolution.
In the text

	Fig. 3. Comparison between the GDE for a DB model (solid line, 1.0 M⊙) and DA (dashed line, 0.5 M⊙). All calculations were performed at τ = 100.0. The horizontal arrow indicates the direction of time evolution.
In the text

	Fig. 4. Effect of log g on the onset of convection for DA models. The continuous line represents a sequence of DA model (1.0 M⊙) while the dashed one also denotes DA models but with 0.52 M⊙. The horizontal arrow indicates the direction of time evolution.
In the text

	Fig. 5. HR diagram for an initial mass of 2.0 M⊙ from the PMS to cooling white dwarf stage. Initial chemical composition X = 0.703, Y = 0.277, αMLT = 1.80.
In the text

	Fig. 6. GDE for the white dwarf cooling sequence for the models shown Fig. 5 (continuous line). The DA model shown in Fig. 1 is represented by a dashed line. The horizontal arrow indicates the direction of time evolution.
In the text