© ESO, 2010

1. Introduction

1.1. Review of observational approaches

One of the most enduring unknowns in stellar physics has been the inner distribution of the angular momentum in a star. In the past few decades, significant progress has been made in describing theoretically the evolution of rotating stars. This has required an understanding of numerous hydrodynamic and magnetic instabilities triggered by the rotation, as well as the mixing processes of chemical elements unleashed by these instabilities (Tassoul 1978, 2000; Zahn 1983, 1992; Maeder & Eenens 2004; Maeder 2009). However, apart from the Sun, reliable observational information about the internal rotation of stars remains scarce or non-existent.

Nevertheless, many attempts have been made to obtain information on the internal rotation from detailed studies of: a) the position of stars in the HR diagram; b) the evolution of the Vsini parameter during the main sequence (MS) phase; c) the shape of absorption lines, whose characteristics can depend upon the rotational law in layers close to the stellar surface; d) the global stellar geometry described with interferometric data.

We briefly review these efforts:

• a)

  The most numerous efforts among those just mentioned are the statistical analysis of photometric data on the rotational spread of the MS, which can be described by (Roxburgh et al. 1966; Maeder 1968; Collins & Sonneborn 1977; Collins et al. 1991) $\mathrm{\Delta }{\mathit{M}}_{\mathrm{V}}\mathrm{=}\mathit{k}\mathrm{\left(}\mathit{n}\mathrm{\right)}{\mathit{V}}^{\mathit{n}}\mathit{,}$(1)where ΔM_V is the deviation in absolute magnitude from the zero-rotation MS, V is the true equatorial rotational velocity, and k(n) is a constant whose value depends on the power n. When n = 2, k(n = 2) is on the order of k_o10^-5 mag/(km s^-1)² so that for k_o ≲ 1 the deviations may indicate that the internal rotation is uniform, while for k_o ≳ 1 the internal rotation can be differential (Cotton & Smith 1983). This type of analysis found that stars do not seem to rotate uniformly. However, owing to the measurement uncertainties and difficulties in defining the MS of zero rotation, the available data could not provide any firm evidence of a particular law of non-uniform rotation (Strittmatter & Sargent 1966; Golay 1968; Maeder 1968; Maeder & Peytremann 1970; Smith 1971; Smith & Worley 1974; Moss & Smith 1981). Furthermore, using detailed model atmospheres for differentially rotating stars, Collins & Smith (1985) concluded that photometry alone can place albeit rather weak constraints on the degree of differential rotation within the stars.

• b)

  Depending upon the internal angular momentum redistribution and evolutionary rearrangements of the inertial momentum, the surface equatorial rotational velocity of stars changes accordingly. Thus, the study of the variation in the true rotational velocity, V, as a function of time was studied by several authors using the ratio ${\mathit{R}}_{\mathrm{LC}}\mathrm{=}\frac{\mathrm{⟨}\mathit{V}\sin \mathit{i}{\mathrm{⟩}}_{\mathrm{LC}}}{\mathrm{⟨}\mathit{V}\sin \mathit{i}{\mathrm{⟩}}_{\mathrm{ZAMS}}}\mathrm{=}\frac{\mathrm{⟨}\mathit{V}{\mathrm{⟩}}_{\mathrm{LC}}}{\mathrm{⟨}\mathit{V}{\mathrm{⟩}}_{\mathrm{ZAMS}}}\mathit{,}$(2)where  ⟨ Vsini ⟩ _LC is the average of the Vsini parameters of stars with in principle the same mass and luminosity class (LC),  ⟨ Vsini ⟩ _ZAMS is the average of Vsini for stars with the same mass, but located near the zero-age-main sequence (ZAMS). These ratios were compared with similar ones predicted theoretically for stars evolving as rotators in two different and extreme ways. On the one hand, the stars were assumed to evolve all their way as uniform rotators, which implies that the angular momentum is entirely redistributed at each evolutionary step. On the other hand, it was assumed that each stellar layer conserved its initial specific angular momentum, i.e., the stars did not undergo any redistribution of its internal angular momentum. Since in many cases the observed ratios R_LC were found to be situated in-between the two extreme theoretical predictions, it was suggested that stars should be differential rotators. However, these studies could not provide any information about the characteristics of the internal rotational law (Sandage 1955; Danziger & Faber 1972; Zorec et al. 1987; Zorec 2004). Somewhat related to this category of inquiries is the study of the evolution of the total angular momentum of B and Be stars carried by Zorec et al. (1990), who concluded that these objects should undergo some internal angular momentum redistribution to explain the observed evolution of the Vsini parameters.

• c₁)

  The study of the absorption line profiles of MS B-type stars found some evidence of possible surface differential rotation. The angular velocity in the surface of stars was assumed to depend on the colatitude angle, θ, as $\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{o}}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{S}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{\left)}\mathrm{\right]}\mathit{,}$(3)where Ω_o is the equatorial angular velocity, R(θ) is the equation of the stellar surface, S is the parameter that testifies to the differential rotation. Using stellar models that are more or less gravity darkened, Stoeckley (1968a) and Stoeckley & Buscombe (1987) found that in most cases S < 0, which suggested that the angular velocity tends to increase from the equator to the pole. Nevertheless, a dependence of the surface angular velocity on the latitude Ω(θ) could be due either to an actual differential rotation present under the stellar surface, or simply to zonal atmospheric currents, which could appear in rapidly rotating early-type stars, as speculated by Cranmer & Collins (1993). In Sect. 2.4, we recall that an acceleration of the angular velocity towards the equator, i.e. S > 0, can be promoted by a temperature gradient induced by the gravity darkening effect.

• c₂)

  The possibility of detecting surface differential rotation by means of the Fourier analysis of spectral line profiles was discussed by Huang (1961), Gray (1977), Bruning (1981), Garcia-Alegre et al. (1982), and Reiners & Schmitt (2002). Evidence of surface differential rotation in late-type stars with Vsini < 50 km s^-1 were given by Reiners & Schmitt (2003b), Reiners & Schmitt (2004), and Reiners & Royer (2004), but no differential rotation for late-type stars with Vsini > 50 km s^-1 and A-type stars with Vsini > 150 km s^-1 were reported by Reiners & Schmitt (2003a) and Gray (1977), respectively. It is possible that modest differential rotation is difficult to detect with the Fourier transform technique in slowly rotating A-type stars, because the rotational broadening is not large compared with the broadening caused by other mechanisms such as thermal turbulence and pressure effects (Gray 1977). However, in those cases where there is some evidence of differential rotation, the parameter S cannot be differentiated from the unknown inclination angle factor sini. Nevertheless, its sign indicates acceleration of the angular velocity towards the equator. To our knowledge, the Fourier technique for differential rotation has not yet been applied to early-type stars.

• c₃)

  Ando (1980) proposed a method to probe the inner angular velocity of stars based on the use of the rotational splitting of non-radial oscillations. However, owing to uncertainties in the identification of pulsation modes and rough determinations of stellar fundamental parameters, particularly their evolutionary stage, this method has not yet been able to be applied with reliable success. Nevertheless, from the analysis of pulsation modes derived from photometric variations (Dupret et al. 2004) and COROT data (Degroote et al. 2009), constraints on the internal rotation of β Cep stars have been inferred.

• d)

  In the past few years, interferometric methods have helped provide remarkable insights into not only the rotational distortion of stars (Domiciano de Souza et al. 2003; van Belle et al. 2004), but also the induced gravity darkening effect by means of imaging techniques (McAlister et al. 2005; Aufdenberg et al. 2006; van Belle et al. 2006; Zhao et al. 2009). New instruments with higher spectral resolutions of up to the 10 000 attained by VLTI/AMBER in the J and K bands and an angular resolution of about 1 mas in the K band (Petrov et al. 2007), or spectral resolution reaching 30 000 and angular resolutions as high as 0.3 mas in the visible using the VEGA/CHARA interferometric array (Mourard et al. 2009), will not only probably enable us to determine with greater detail than in previous studies the global geometry of stars deformed and gravity darkened by the rotation, but also carry out differential interferometry.

A method based on differential interferometry that requires high spectral and spatial resolution was presented by Domiciano de Souza et al. (2004b,c); Domiciano de Souza et al. (2004a) to distinguish observationally the parameter controlled by the degree of the surface differential rotation from the inclination angle factor sini.

1.2. Aims of the present attempt

Most theoretical predictions about the evolution of rotating early-type stars and the mixing of chemical species triggered by the instabilities set up by the rotation, come from calculations performed in the framework of two significant assumptions: a) the global rotational energy stored by the stars in the ZAMS is lower than the limit allowed by the rigid rotation in the critical regime; b) the internal angular velocity undergoes an instantaneous “shellular” redistribution at each evolutionary step. However, Clement (1979) using a cylindrical (conservative) distribution of the angular velocity, and Maeder et al. (2008) basing their calculation on a “shellular” distribution, showed in a more detailed way that in rapidly rotating early-type stars the envelope layers beneath the surface, may have wider convective zones in radius than in non-rotating stars. In the Sun, only the layers unstable to convection rotate differentially with a non-shellular pattern. This motivates the inquiry of whether in massive and intermediate mass stars some coupling may also exist between convection and rotation beneath their surface. In this case, the characteristics of the rotational law in the external stellar layers should differ from those currently assumed in the above evoked stellar models.

As demonstrated by many authors, the global geometry of a star depends not only on the total amount of angular momentum stored by the star, but also on its internal distribution (Bodenheimer 1971; Zorec 1986; Smith & Collins 1992; Uryu & Eriguchi 1994, 1995). This geometry mostly relies on the stellar surface rotation, which acts as an imprint of its properties in the layers beneath the surface. In this case, we should not exclude the resulting mixing of chemical elements in the stellar atmosphere being more or less dependent on the characteristics of the external rotational law, upon which the description of the stellar structure, based on the abundance determination of chemical elements, should also rely. Therefore, to provide new information and/or constraints to test the global assumptions currently made to calculate models of stellar structure with rotation and thus help deepen our understanding of the properties of early-type fast rotators, we might ask: 1) what can be deduced, using first principles, about the properties the rotation laws can have beneath the surface as a consequence of the coupling between rotation and convection; 2) what are the parameters needed to characterize these stars that may be accessible to observations; 3) whether the combined interpretations of spectroscopic and interferometric data of rapidly rotating early-type stars enable us to determine these parameters.

In this attempt, the most interesting information might be the indication of some differential rotation in the stellar surfaces and the sign of its latitudinal gradient. Both pieces of data can be obtained, as much as possible, in a consistent way by taking into account the stellar geometrical deformation produced by this rotation and the concomitant gravitational darkening effect that responds to possible non-conservative rotation laws.

The present paper is organized as follows. In Sect. 2, we use first principles to infer possible rotation laws in the convective layers beneath the stellar surface of early-type rapid rotators. Section 3 presents the equation of the stellar surface of stars with non-conservative rotation laws. A discussion of the gravity darkening effect for non-conservative rotation laws is presented in Sect. 4. The discussion about the validity of the Roche approximation in representing the gravitational potential is presented in Sect. 5. This discussion is based on 2D models of rotating stars where the evolutionary stages are taken into account in a simplified way. We briefly comment on the determination of the rotational profile in the stellar envelope in Sect. 7. In Sect. 6, we summarize the attainable information on rapidly rotating early-type stars with external differential rotation from the combined analysis of spectroscopic and interferometric data. Our conclusions are presented in Sect. 8.

2. Rotational law in the stellar envelope

2.1. The angular velocity distribution beneath the stellar surface

The effects of rotation are generally introduced in the structure equations of rotating stars by replacing the spherical stratification of non-rotating star-models by a rotationally distorted stratification, which keeps the whole calculation problem to one dimension (Kippenhahn & Thomas 1970; Endal & Sofia 1976; Pinsonneault et al. 1990; Fliegner & Langer 1995; Chaboyer et al. 1995; Meynet & Maeder 1997). This procedure is justified if the internal differential rotation has a shellular distribution law because it comes from theoretical inferences made by Zahn (1983, 1992), which rely on the assumption that the horizontal turbulence is much stronger than the vertical one. The limiting case of a shellular rotational profile is rigid rotation. From the calculations by Maeder & Meynet (2005), it follows that magnetic fields created by the Pitts & Tayler instability (Pitts & Tayler 1985; Spruit 1999, 2002) can lock the stellar layers to each other and force the star to evolve as a rigid rotator. However, Zahn et al. (2007) concluded that the dynamo action can be less efficient as previously expected (Spruit 1999, 2002) and that the magnetic fields created contribute little or less to the transport of the angular momentum. In addition, if the magnetic field in early-type rapid rotators is created, it could perhaps has some effect on their rotational profile in the convective parts of the envelope as in the Sun (Balbus 2009; Balbus et al. 2009).

The differential rotation in the surface of the Sun is a direct consequence of the differential rotation in the convective layers beneath the surface layers. Therefore, it is important to stress that this differential rotation in depth is not shellular, in spite of the strong turbulent viscosity that according to the above-mentioned theoretical assumptions, might otherwise cause shellular-like rotation.

Maeder et al. (2008) demonstrated that rotation does not inhibit convection as could be thought from Solberg-Høiland’s stability criterion, but it changes the thermal gradient so as to enhance convection. Hence, in rapidly rotating massive stars the two external convection zones, associated with increased opacities due to He- and Fe-ionization, respectively, are both considerably enlarged in depth so that the entire convective zone covers a non-negligible external region, which ranges from some 1/8 of the stellar radius in the pole to nearly 1/4 in the equator.

However, the angular momentum distribution in the convective regions remains a puzzling question. To account for it in the stellar model calculations, two extreme approximations have been used: a) rigid rotation, angular velocity Ω =  constant, which is supposedly promoted by the turbulent viscosity (Maeder et al. 2008); b) constant specific angular momentum j (j = Ωϖ²; ϖ is the distance to the axis of rotation), possibly due to the redistribution produced by the convective plumes (Tayler 1973). Nevertheless, the solar convective regions, characterized by significant turbulence, are the only ones with significant differential rotation (Schou et al. 1998), even though the solar rotational profile in the convective region, revealed by the helioseismological data, does not fall between these two extreme possibilities (Deupree 2001). Global insights into the rotational law in the convective regions could be obtained by exploring the solutions to the thermal wind balance relation under imposed conditions between the entropy and rotation. This procedure successfully explains the Solar rotation law in the convective regions (Balbus 2009; Balbus et al. 2009). However, we pay attention here to our dealing with rapidly rotating stars. Detailed physical justifications of the assumptions made in the present paper are postponed to future contributions.

2.2. Rotation inferred from the baroclinic balance relation

The curl of the time-independent momentum equation of an inviscid, axisymmetric rotating star with negligible magnetic fields, yields the baroclinic balance relation. Using the cylindrical coordinates (ϖ,φ,z), this balance condition reads (Tassoul 1978) $\frac{\mathrm{1}}{{\mathit{\varpi }}^{\mathrm{3}}}\frac{\mathit{\partial }{\mathit{j}}^{\mathrm{2}}}{\mathit{\partial z}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{2}}}\mathrm{(}\mathrm{\nabla }\mathit{P}\mathrm{\times }\mathrm{\nabla }\mathit{\rho }\mathrm{)}\mathrm{ê}\mathit{\phi }\mathit{,}$(4)where j is the specific angular momentum (j = Ωϖ²), P and ρ are pressure and density, respectively, and ê_φ is the azimuthal unit vector. Taking into account the equation of hydrostatic equilibrium $\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{P}\mathrm{=}{g}\mathrm{eff}\mathit{,}$(5)where g_eff is the effective gravity, and using for the specific entropy S the expression $\mathit{S}\mathrm{=}\frac{\mathit{k}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\ln \mathit{P}{\mathit{\rho }}^{\mathrm{-}\mathit{\gamma }}\mathrm{+}\mathrm{const}\mathit{.}\mathit{,}$(6)where k is the Boltzmann constant and γ is the ratio of specific heats at constant pressure and constant volume per unit mass, the wind equation in Eq. (4) can be rewritten as $\frac{\mathrm{1}}{{\mathit{\varpi }}^{\mathrm{3}}}\frac{\mathit{\partial }{\mathit{j}}^{\mathrm{2}}}{\mathit{\partial z}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{C}}_{\mathrm{P}}}\mathrm{(}\mathrm{\nabla }\mathit{S}\mathrm{\times }{g}\mathrm{eff}\mathrm{)}\mathrm{ê}\mathit{\phi }\mathit{,}$(7)where C_P is the constant pressure specific heat.

We can attempt a discussion of Eq. (4) by seeking solutions for the stellar internal rotation under at least three different conditions where for the moment the effects carried by the meridional circulation are neglected: a) marginal stability imposed by the Solberg-Høiland criterion; b) state enforced by parallel surfaces of specific entropy and specific angular momentum; c) frame where the specific entropy parallels the local specific kinetic rotational energy.

• a)

  The Solberg-Høiland stability criterion states that “a baroclinic star in permanent rotation is dynamically stable against axisymmetric perturbations if two conditions are satisfied: (i) the specific entropy S increases outwards, and (ii) on each surface S = const., the specific angular momentum increases from the pole to the equator”. The second condition is written mathematically as (Tassoul 1978) $\mathrm{-}{\mathit{g}}_{\mathit{z}}\mathrm{[}\mathrm{\nabla }{\mathit{j}}^{\mathrm{2}}\mathrm{\times }\mathrm{\nabla }\mathit{S}\mathrm{]}\mathrm{\ge }\mathrm{0}\mathit{,}$(8)where g_z is the z − component of the effective gravity. Making the ansatz for strict marginal equilibrium, we have $\mathrm{-}{\mathit{g}}_{\mathit{z}}\mathrm{[}\mathrm{\nabla }{\mathit{j}}^{\mathrm{2}}\mathrm{\times }\mathrm{\nabla }\mathit{S}\mathrm{]}\mathrm{=}\mathrm{0.}$(9)Equation (9) suggests then that surfaces j² = const. and S = const. should be parallel, S = S(j²), i.e. a displaced fluid element in baroclinic turbulence retains both entropy and angular momentum. We do not consider here the conditions that might render possible the balance implied by Eq. (9).

• b)

  The assumption that the surfaces of constant specific entropy and constant angular velocity coincide, i.e. S = S(Ω²), brings another alternative solution to the baroclinic equilibrium Eq. (4). This coupling can be enforced by magnetic fields (Balbus 2009), although hydrodynamic constraints in the Sun can justify it entirely (Balbus et al. 2009).

• c)

  Only as an extrapolation to the conditions S = S(j²) and S = S(Ω²) may we also consider S = S(ϖ²Ω²), since two energy-related quantities are parallel: specific entropy S and specific rotational kinetic energy ϵ_Ω = ϖ²Ω².

We now use H to represent in turn j², Ω², and ϵ_Ω. We have then S = S(H) and ∇S = (dS/dH)∇H. Rewriting Eq. (7) in spherical coordinates (r,θ,φ) and knowing that in these coordinates the r- and θ-components of the effective gravity g_eff are $\begin{array}{ccc}{\mathit{g}}_{\mathit{r}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{G}}}{\mathit{\partial r}}\mathrm{+}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathit{r}{\sin }^{\mathrm{2}}\mathit{\theta ,}\\ {\mathit{g}}_{\mathit{\theta }}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{G}}}{\mathit{\partial \theta }}\mathrm{+}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathit{r}\sin \mathit{\theta }\cos \mathit{\theta ,}\end{array}$(10)we readily obtain $\frac{\mathit{\partial }H}{\mathit{\partial r}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{H}H\mathrm{)}\mathrm{-}\frac{\tan \mathit{\theta }}{\mathit{r}}\frac{\mathit{\partial }H}{\mathit{\partial \theta }}\left[\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{H}H\mathrm{+}{\mathit{\alpha }}_H\mathit{GM}F\right]\mathrm{=}\mathrm{0}\mathit{,}$(11)where we have used the notations ${\mathit{\alpha }}_H\mathrm{=}\frac{\mathrm{1}}{{\mathit{C}}_{\mathrm{P}}}\frac{\mathrm{d}\mathit{S}}{\mathrm{d}H}$(12)and $F\mathrm{=}\{\begin{array}{c}\\ \mathit{r}{\sin }^{\mathrm{2}}\mathit{\theta }& \mathrm{if}H\mathrm{=}{\mathit{j}}^{\mathrm{2}}\\ \frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{3}}{\sin }^{\mathrm{2}}\mathit{\theta }}& \mathrm{if}H\mathrm{=}{\mathrm{\Omega }}^{\mathrm{2}}\\ \frac{\mathrm{1}}{\mathit{r}}& \mathrm{if}H\mathrm{=}{\mathit{ϵ}}_{\mathrm{\Omega }}\mathit{.}\end{array}$(13)The characteristic equation of Eq. (11) is $\frac{\mathit{dr}}{\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{H}H}\mathrm{=}\mathrm{-}\frac{\mathit{d\theta }}{\frac{\tan \mathit{\theta }}{\mathit{r}}\left[\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{H}H\mathrm{)}\mathrm{+}{\mathit{\alpha }}_H\mathit{GM}F\right]}\mathrm{·}$(14)As for the assumption S = S(H) made here, H and S are constant along the characteristic curves of Eq. (11), and the differential equation Eq. (14) integrates immediately to give $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}{\mathrm{{\rm Y}}}_{\mathrm{J}}\mathrm{+}{\mathit{C}}_{\mathrm{J}}\mathrm{if}H\mathrm{=}{\mathit{j}}^{\mathrm{2}}\\ \mathit{r}{\sin }^{\mathrm{2}}\mathit{\theta }& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathit{r}}{\mathrm{{\rm Y}}}_{\mathrm{\Omega }}\mathrm{+}{\mathit{C}}_{\mathrm{\Omega }}\mathrm{if}H\mathrm{=}{\mathrm{\Omega }}^{\mathrm{2}}\\ \frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}& \mathrm{=}& \end{array}$(15)with the C_H as the integration constants (where H has different meanings).

Once the function giving the angular velocity Ω_s(θ) in the stellar surface is specified, the iso-rotation curves can be obtained easily everywhere inside the star using the relations in Eq. (15). Each iso-rotation curve depends on an integration constant C_H obtained as C_H = C_H(Υ_H,θ_s), where θ_s is the colatitude angle at which a given iso-rotation contour intersects the stellar surface. The constants Υ_H are defined as ${\mathrm{{\rm Y}}}_H\mathrm{=}\mathrm{-}\left(\frac{\mathrm{2}{\mathit{\alpha }}_H\mathit{GM}}{\mathrm{1}\mathrm{-}{\mathit{\alpha }}_H}\right)\mathrm{·}$(16)Since we consider convective regions where ∇S ≤ 0, we find that Υ_H ≥ 0. We note that depending on the meaning of H, the expressions Υ_H =  Υ_JR_e, ${\mathrm{{\rm Y}}}_{\mathrm{\Omega }}\mathit{/}{\mathit{R}}_{\mathrm{e}}^{\mathrm{2}}$ and Υ_ϵΩ are dimensionless constants of the order of unity (R_e is the stellar radius at the equator). Defining ${\mathit{A}}_H\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\gamma }}\left(\frac{\mathrm{d}\ln \mathit{P}{\mathit{\rho }}^{\mathrm{-}\mathit{\gamma }}}{\mathrm{d}\ln \mathit{r}}\right)\left(\frac{\mathrm{d}\ln \mathit{r}}{\mathrm{d}\ln H}\right)\mathit{,}$(17)which in convection regions is A_H ≥ 0, we have ${\mathrm{{\rm Y}}}_H\mathrm{=}\mathrm{2}\left(\frac{\mathit{GM}}{{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}{\mathrm{\Omega }}^{\mathrm{2}}}\right)\left(\frac{{\mathit{A}}_H}{\mathrm{1}\mathrm{-}{\mathit{A}}_H}\right)\mathrm{·}$(18)From the model of a rapidly rotating star with M = 20   M_⊙ given in Maeder et al. (2008), we derived the adiabatic gradient ∇ = dlnPρ^−γ/dlnr shown in Fig. 1. This quantity is shown as a function of stellar radius in the polar and equatorial directions. The value of Υ_H depends also on the adopted value for the gradient dlnr/dlnΩ² that represents the differential rotation. For the Sun, the term dlnr/dlnΩ² comes from seismology data and its value ranges from 1 to 10. However, until now in early-type stars, the factor dlnr/dlnΩ² has been a perfectly unknown quantity. For a rough estimate, we may assume that it can be at least of the same order of magnitude as in the Sun. Similar values are also obtained in the models calculated by Meynet & Maeder (2000). Using the values of the gradient ∇ shown in Fig. 1, it can be shown that the ratio A_H/(1 − A_H) is nearly constant over the entire convective region and that depending on the assumed value of dlnr/dlnΩ², it changes from 0 to 1. Hence, we have that $\mathrm{0}\lesssim {\mathrm{{\rm Y}}}_H\lesssim \mathrm{2}\mathrm{(}\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{)}\mathit{,}$(19)where $\mathrm{0.5}\lesssim \mathrm{(}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{/}\mathit{GM}\mathrm{)}\lesssim \mathrm{1}$ for rapid rotators. For the sake of comparison, we note that in the convective region of the Sun, it is |∇| ≃ 10^-6 and $\mathrm{(}\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{)}\mathrm{\simeq }{\mathrm{10}}^{\mathrm{5}}$ (Balbus 2009).

Examples of isorotation curves obtained from Eq. (15) are shown in Fig. 2, where we used stellar external contours whose rotational deformation is characterized by the ratio of centrifugal to gravitational acceleration at the equator ${\mathit{\eta }}_{\mathrm{o}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{e}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}\mathit{/}\mathit{GM}\mathrm{=}\mathrm{0.8}$ and a Maunder surface differential rotation law with parameter α = 0.3. The equation of the stellar surface with latitudinal differential rotation is discussed in some detail in Sect. 3. Solutions similar to those shown in Figs. 2a and 2b were obtained by Balbus (2009) for the Sun. This author assumed that ∂P/∂r ≫ ∂S/∂θ and ∂S/∂θ ≫ ∂P/∂r, which according to our expressions are for the limits α_Jj² → 0 and α_ΩΩ² → 0 that suit slowly rotating objects like the Sun.

2.3. Comments on the rotational profiles obtained

The rotational profiles Ω(r,θ) in the convective layers beneath the surface of early-type rapidly rotating stars were inferred here using the marginal condition of the Solberg-Høiland stability criterion against axisymmetric perturbations. The derived solutions are the consequence of a dominant thermal wind balance where the entropy is: a) a function of the specific angular momentum, S = S(j²); b) a function of the angular velocity, S = S(Ω²), and c) a function of the specific rotational kinetic energy, S = S(ϵ_Ω). The first case was widely studied in the literature (Tassoul 1978,references therein), and a number of independent numerical attempts to obtain baroclinic stars inexorably ended up with results close to the category shown in Figs. 2a₁, i.e. conservative-like (Uryu & Eriguchi 1994, 1995). Two-dimensional, implicit hydrodynamic numerical simulations by Deupree (1998) and Deupree (2001) showed that in convective regions Ω tends to adopt a cylindrical-like profile rather than a shellular or Ω =  constant form. Intermediate solutions to the classes shown in Figs. 2b₁ and b₂ are compatible with the seismological data of the Sun relative to its convective regions. The rotational profile obtained by Espinosa Lara & Rieutord (2007) for fully radiative stars curiously enters into the category of solutions shown in Fig. 2c₂. We note that for all cases shown in Fig. 2 we have used constant values of Υ_H all over the star. However, they are valid only for the convective regions, which lie roughly above the shaded sector and where, as already noted, Υ_H is fairly uniform.

According to a strong claim that shear instabilities widely prevail over all other mechanisms that can act to redistribute the angular velocity, most if not all models of stellar evolution of early-type stars with rotation were carried out by assuming that Ω has a shellular nature. This strongly simplifies the numerical aspect of the study, but does not necessarily preclude that other distributions of Ω might coexist with the shellular one, or even dominate in wide domains of the stellar interior. The S = S(Ω²) − case with Υ_Ω > 1.6 clearly depicts shellular-like rotation over a wide interior stellar region.

Although the present discussion should not be considered an argument in favor of non-conservative-like profiles of the angular velocity in the envelopes of early-type rapid rotators, the success that the displayed arguments do have in explaining the solar rotational profile in the convective layers, strongly motivate us to address the study of the surface differential rotation in rapidly rotating early-type stars. Whatever the observational or theoretical indication that ∂Ω/∂θ ≠ 0 can exist in the surface of rapidly rotating early-type stars, they would immediately imply that profiles of non-shellular rotation exist, perhaps of some type similar to that in Fig. 2. They would also immediately enable us to achieve clearer understandings of the transport of angular momentum in the stars and the related mixing processes of chemical elements in the stellar surface layers.

2.4. Induced gradient in Ω(θ) in the surface of early-type rapid rotators

In rapidly rotating stars, there is a strong latitudinal temperature gradient, which can induce a latitudinal gradient on the surface angular velocity Ω_s. This effect can be estimated from Eq. (4) by calculating the latitudinal variation in the temperature over an isobar. Giving to the surface temperature a functional form T = T(P,θ), where P is the pressure and θ the colatitude, we obtain $\mathrm{\nabla }\mathit{T}\mathrm{\left(}\mathit{P,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathit{\partial T}}{\mathit{\partial P}}\mathrm{\nabla }\mathit{P}\mathrm{+}\frac{\mathit{\partial T}}{\mathit{\partial \theta }}\frac{\mathrm{ê}\mathit{\theta }}{\mathit{r}}\mathit{,}$(20)where ê_θ is the colatitude unit vector and r is the radial spherical coordinate. Making use of the equation of state P ∝ ρT, Eq. (20) and the r-component of Eq. (5), Eq. (4) transforms into (Espinosa Lara & Rieutord 2007) $\mathit{\varpi }\frac{\mathit{\partial }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{{\mathit{g}}_{\mathit{r}}}{\mathit{r}}{\left(\frac{\mathit{\partial }\ln \mathit{T}}{\mathit{\partial \theta }}\right)}_{\mathrm{P}}\mathit{,}$(21)where the notation (...)_P indicates that the partial θ-derivative is along a barotropic surface. As in rapidly rotating stars, it is ∂lnT/∂θ < 0 and g_r < 0, from Eq. (21) it appears that ∂Ω/∂z < 0, which means that Ω should increase from the pole to the equator. Although in rapid rotators there are convective zones in the envelope, the dominant radiation flux can still be estimated using the diffusion approximation F_rad = −χ∇T, where χ is the coefficient of radiative conductivity. For simplicity, we adopt von Zeipel’s approximation (von Zeipel 1924) to represent the stellar surface temperature distribution as a function of the colatitude $\mathit{T}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathit{p}}{\left(\frac{\mathit{g}}{{\mathit{g}}_{\mathit{p}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(22)where T_p is the polar temperature of the star, g =  [(g_r)² + (g_θ)²] ^1/2 with g_r and g_θ given by Eq. (10), and g_p is the value of the gravity in the pole. Using the Roche approximation for the gravitational potential, Eq. (21) written in spherical coordinates becomes $\begin{array}{ccc}\frac{\mathit{\partial }\ln {\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\partial r}}\mathrm{-}\frac{\tan \mathit{\theta }}{\mathit{r}}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\mathrm{\left)}{\mathit{f}}_{\mathrm{1}}\mathrm{\right(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\left)}\mathrm{\right]}\frac{\mathit{\partial }\ln {\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\partial \theta }}\mathrm{=}& & \\ & & \mathrm{-}\frac{\mathrm{1}}{\mathit{r}}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}\end{array}$(23)with the following definitions $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}}\left[\frac{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}{\sin }^{\mathrm{2}}\mathit{\theta }}{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\mathrm{(}\mathrm{2}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\mathrm{)}{\sin }^{\mathrm{2}}\mathit{\theta }}\right]\\ {\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}& \mathrm{=}& \left[\frac{\mathrm{2}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\eta }}_{\mathit{r}}\mathrm{)}{\sin }^{\mathrm{2}}\mathit{\theta }}{\sin \mathit{\theta }\cos \mathit{\theta }}\right]\frac{\mathit{\partial }\ln \mathit{r}}{\mathit{\partial \theta }}\mathrm{+}\mathrm{2}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\\ {\mathit{\eta }}_{\mathit{r}}& \mathrm{=}& \frac{{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}}{\mathit{GM}}\mathrm{·}\end{array}$(24)The characteristic equation describing the dependence of the surface angular velocity with θ is then $\mathrm{-}\frac{\mathrm{d}\mathit{\theta }}{\frac{\tan \mathit{\theta }}{\mathit{r}}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathit{r}}\mathrm{\left)}{\mathit{f}}_{\mathrm{1}}\mathrm{\right(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\left)}\mathrm{\right]}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}\ln {\mathrm{\Omega }}^{\mathrm{2}}}{\frac{\mathrm{1}}{\mathit{r}}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathit{r}}\mathrm{\right)}}\mathit{,}$(25)which is solved by numerical integration. Figure 3 shows the θ-dependence of the angular velocity given as ΔΩ(θ) =  [Ω(θ,η_r) − Ω(40°,η_r)] /Ω(40°,η_r), which was obtained from Eq. (25) for different values of η. In this calculation, the gradient ∂lnr/∂θ in the stellar surface was estimated using the Roche potential of a rigid rotator characterized by average rotational parameters η. In Fig. 3, we see that in general it is ∂Ω/∂θ > 0, in particular for η_o > 0.7 near the equator. This acceleration seems to be similar in character to the one previously calculated by Espinosa Lara & Rieutord (2007). As η → 0, the curves ΔΩ(θ) → 0.

3. Equation of the stellar surface

3.1. The angular velocity distribution in the stellar surface

Since the shape of a rapidly rotating star depends explicitly on its surface angular velocity distribution, we adopt a law to study the possible effects it can produce. In Sect. 2.4, we have seen that the temperature gradient induced by the gravity darkening effect introduces a small acceleration of the angular velocity from the pole towards the equator, which is not simple to represent analytically. The physical properties of the layers beneath the surface could perhaps enforce a stronger surface angular velocity gradient than inferred in Sect. 2.4. If in rapidly rotating massive stars rotational profiles of the class obtained in Sect. 2.2 actually existed, a first approach to describing their surface rotation could rely on the use of a solar-like surface angular velocity. The solar surface differential rotation depends on the colatitude θ as (Snodgrass 1984) $\frac{{\mathrm{\Omega }}_{\mathrm{\odot }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{2}\mathit{\pi }}\mathrm{=}\mathrm{451.5}\mathrm{-}\mathrm{65.3}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}\mathrm{66.7}{\cos }^{\mathrm{4}}\mathit{\theta }\mathrm{nHz}\mathit{.}$(26)This expression carries three coefficients that need to be determined empirically. At the moment, it seems unrealistic to use a similar expression for other stars since their surface cannot be resolved. However, within the errors smaller than 1.6%, Eq. (26) can be reduced to the simplest one $\frac{{\mathrm{\Omega }}_{\mathrm{\odot }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{2}\mathit{\pi }}\mathrm{=}\mathrm{459.3}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0.29}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathrm{nHz}\mathit{,}$(27)where there are only two quantities to determine from observations. A first inquiry about the differential rotation on the surface of massive and intermediate mass stars can then be justified by using the simplified expression (Maunder formula) $\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{o}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\alpha }{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathit{.}$(28)In Eq. (28), α < 0 indicates that the pole rotates slower than the equator and vice-versa if α > 0. If massive and intermediate mass stars had surface rotations similar to the Sun, their differential rotation parameter would be expected to be at least of the order of α ~ −0.3.

In general, to study stars with unresolved surfaces we should probably use expressions of the type $\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{o}}\mathrm{[}\mathrm{1}\mathrm{+}\mathit{\alpha f}\mathrm{(}\mathit{\theta }\mathrm{\left)}\mathrm{\right]}\mathit{,}$(29)where the function f(θ) needs to be somehow specified in advance.

3.2. Formulation of the equation of the stellar surface

The shape of a rapidly rotating star is generally described using the total potential of a self-gravitating system with a conservative rotation law, i.e., an internal cylindrical angular velocity, Ω = Ω(ϖ), of which the rigid rotation is a particular case. When a star has a non-conservative internal rotational law Ω = Ω(ϖ,z), as used in models of stellar evolution and/or suggested by the developments given in Sect. 2.2, it is no longer possible to define a rotational potential. However, according to Maeder (2009) the surface of a star with a non-conservative rotation law can be identified as the region where an arbitrary displacement ds does not imply any work done by the effective gravity g_eff${g}\mathrm{eff}\mathit{.}\mathrm{d}{s}\mathrm{=}\mathrm{0.}$(30)The effective gravity g_eff written in cylindrical coordinates has the form ${g}\mathrm{eff}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{+}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}\mathit{\varpi }\mathit{.}\mathrm{ê}\mathit{\varpi }\mathit{,}$(31)where Φ_G is the gravitational potential. Using the function Ψ(ϖ,z) defined as $\mathrm{\Psi }\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}{\mathit{\varpi }}^{\mathrm{2}}\mathit{,}$(32)the effective gravity can be expressed in the form ${g}\mathrm{eff}\mathrm{=}\mathrm{-}\mathrm{\left[}\mathrm{\nabla \Psi }\mathrm{\right(}\mathit{\varpi ,z}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varpi }}^{\mathrm{2}}\mathrm{\nabla }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{]}\mathit{,}$(33)so that the condition given in Eq. (30) can be rewritten as $\mathrm{d\Psi }\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varpi }}^{\mathrm{2}}\mathrm{\left(}\mathrm{\nabla }{\mathrm{\Omega }}^{\mathrm{2}}\mathit{.}\mathrm{d}{s}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(34)The function Ψ(ϖ,z) can be a potential if and only if it is ∂Ω/∂z = 0. It has been shown by Meynet & Maeder (1997) that in stars with shellular internal rotational profiles, the surfaces Ψ(ϖ,z) =  constant are parallel to isobar surfaces. However, since dΨ(ϖ,z) is not an exact differential, the integration of Eq. (34) depends on the chosen path. To define the shape of a star, it seems natural to integrate Eq. (34) over a meridian curve, from the pole (spherical coordinate θ = 0) towards an arbitrary point R_s(θ), where R_s(θ) represents the function describing the stellar surface. Thus, by virtue of Eq. (32) we obtain ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}{\mathit{\varpi }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\theta }}{\mathit{\varpi }}^{\mathrm{2}}\mathrm{\left(}\mathrm{\nabla }{\mathrm{\Omega }}^{\mathrm{2}}\mathit{.}\mathrm{d}{s}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{,}$(35)which is in principle the sought equation to calculate the shape of a star having non-conservative rotational laws. Nevertheless, since the integration indicated in Eq. (35) depends on the angular velocity distribution defined only over the stellar surface, the following change in the coordinates describing the stellar surface ${\mathit{\varpi }}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathit{R}\mathrm{(}\mathit{\theta }{\mathrm{)}}_{\mathrm{s}}\sin \mathit{\theta }\mathrm{;}{\mathit{z}}_{\mathrm{s}}\mathrm{(}\mathit{\theta }\mathrm{)}\mathrm{=}{\mathit{R}}_{\mathrm{s}}\mathrm{(}\mathit{\theta }\mathrm{)}\cos \mathit{\theta ,}$(36)readily transforms ∇Ω²(ϖ,z).ds into a total differential that enables us to integrate Eq. (35) by parts and derive the sought relation ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\theta }}{\mathrm{\Omega }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\left(\frac{\mathrm{d}{\mathit{\varpi }}^{\mathrm{2}}}{\mathrm{d}\mathit{\theta }}\right)\mathrm{d}\mathit{\theta }\mathrm{=}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(37)In Eq. (37), we wrote ${\mathrm{\Omega }}_{\mathrm{s}}^{\mathrm{2}}$ to indicate explicitly that the angular velocity concerns only the stellar surface. From Eq. (37), we can see that the stellar shape depends on both the rotational profile in the surface Ω_s(θ) and the gravitational potential Φ_G, which in turn depends on the internal rotational law by means of the centrifugally distorted density distribution in the stellar interior. However, this last dependence is of second order or negligible, as will be shown in Sect. 5.4.

3.2.1. Shape of stars with a Maunder rotation law in the surface

According to the discussion in Sect. 2.4 and the simplified solar surface rotational law (28), the differential rotation parameter in fast rotating early-type stars is expected to be α < 0. However, because of the still poorly known physical characteristics of the convective external envelope layers of these objects, and possible effects linked to magnetic fields, we cannot exclude that α > 0 in some objects. In this respect, we note that hydrodynamic calculations in stars with 8.75   M_⊙ carried out by Deupree (1998) show that the rotation law in the convective core becomes of conservative type Ω ~ ϖ^-0.7 (ϖ being the distance from the rotation axis). If convection in the envelope forced it to rotate in a similar way, we should expect an external imprint of this law given by an α > 0. Many calculations of the stellar structure with fast differential rotation have been done in the past assuming that Høilnad’s criterion for dynamical stability is satisfied and that radiative viscosity has significant effects on the surface layers. The resulting rotation law implies that the specific angular momentum increases with the mass contained in cylinders with radius ϖ (Tassoul 1978), which also leads to rotation laws that imply α > 0. However, the final shape of the rotation law in the envelope must certainly be the consequence of the interplay among many hydrodynamic and magneto-hydrodynamic instabilities, whose very final combined effect is still highly unknown. We adopt then the surface rotational law with the functional form given by Eq. (27), and assume that the parameter α can imply either equator-ward or polar-ward accelerated angular velocities.

In a first step, we separate the effects produced on the external stellar shape by the surface differential rotation from those induced by the internal rotation on the mass distribution. This can simply be done, as we justify in Sect. 5.4, by using a Roche approximation for the gravitational potential ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}$(38)Equation (37) can then be given the form $\frac{{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{\eta }}_{\mathrm{o}}\mathrm{\left[}\mathit{I}\mathrm{\right(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}\mathrm{-}\mathit{I}\mathrm{(}\mathit{\theta }\mathrm{\left)}\mathrm{\right]}}\mathit{,}$(39)where we have written $\mathit{I}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\theta }}{\left[\frac{{\mathrm{\Omega }}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathrm{\Omega }}_{\mathrm{o}}}\right]}^{\mathrm{2}}\left(\frac{\mathrm{d}{\mathit{\varpi }}^{\mathrm{2}}}{\mathrm{d}\mathit{\theta }}\right)\mathrm{d}\mathit{\theta }\mathit{.}$(40)Since in relation (39) the function of the stellar surface R_s(θ) also appears in the integrand in terms of ϖ by (36), R_s(θ) is obtained by iteration. We start the iteration with a first estimate of I(π/2) based on the shape of a star with a rigid rotation characterized by the rotational parameter η_o${\mathit{\eta }}_{\mathrm{o}}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{o}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}}{\mathit{GM}}\mathit{,}$(41)so that $\frac{{\mathit{R}}_{\mathrm{s}}^{\mathrm{\left(}\mathrm{o}\mathrm{\right)}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}_{\mathrm{o}}\left\{\mathrm{1}\mathrm{-}{\left[\frac{{\mathit{R}}_{\mathrm{s}}^{\mathrm{\left(}\mathit{o}\mathrm{\right)}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\right]}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }\right\}}\mathrm{·}$(42)It is also obvious that in relation (39) the equatorial-to-polar radii ratio R_e/R_p is the solution of the iteration, which finally gives $\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{\eta }}_{\mathrm{o}}\mathit{I}\left(\frac{\mathit{\pi }}{\mathrm{2}}\mathit{,}\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\mathit{,\alpha }\right)\mathit{.}$(43)We recall that for rigid rotation in the Roche approximation we have ${\left(\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\right)}_{\mathrm{rigid}}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}_{\mathrm{o}}\mathit{.}$(44)In Table 1, we present the ratios of equator-to-polar radii R_e/R_p as a function of η_o and α used to obtain the stellar shapes shown in Fig. 4, which shows the shapes of stars computed with equatorial acceleration ratios η_o = 0.7 and 0.9, whose surface angular velocity is given by Eq. (28), where the parameter α has several positive and negative values. In this figure, we can see that the stellar shape is a function of α that sensitively differs from that of homologous rigid rotators with the same parameter η_o. We note that for Ω(θ) with α > 0 (acceleration from the equator towards the pole), we always find that $\frac{{\mathit{R}}_{\mathrm{e}}\mathrm{\left(}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{p}}\mathrm{\left(}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha }\mathrm{\right)}}\mathit{>}{\left[\frac{{\mathit{R}}_{\mathrm{e}}\mathrm{\left(}{\mathit{\eta }}_{\mathrm{o}}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{p}}\mathrm{\left(}{\mathit{\eta }}_{\mathrm{o}}\mathrm{\right)}}\right]}_{\mathrm{rigid}}\mathit{.}$(45)For modest α > 0 (recalling that in the Sun, it is α ≃ ₀.3), we can have R_e/R_p > 1.5 even for η_o < 1, although the stellar surface has no polar dimples as in models with high rotational energy content (See Sect. 3.2.2). This is due both to a rotational stretching of R_e and to the concomitant shrinkage of the polar radius R_p, where Ω_pole > Ω_equator. When Ω_s(θ) is accelerated from the pole towards the equator (α < 0), we have R_e(η_o, − |α|)/R_p(η_o, −|α|) < (R_e(η_o)/R_p(η_o))_rigid. Here, the equator/pole radii ratios respond mainly to the rotational stretching of R_e.

We note that in stars where α ≠ 0, the Vsini parameter does not only depend on the rotation in the equator. A straightforward interpretation of the quantity Vsini can produce an incorrect value of the acceleration ratio η_o. Since we are dealing here with fast rotators (η_o ≳ 0.5), the formation of spectral lines used to determine the Vsini needs to be treated properly by taking into account the surface velocity fields, the stellar deformation, and the gravity darkening effect, all dependent on η_o and α. We also note that if the interferometric measurements produce radii ratios R_e/R_p > 1.5, it does not necessarily mean that the star is in a state of “supercritical” rotation. From these comments and the results shown in Table 1, we clearly see that the study of the surface differential rotation prefigures the need for two sources of information: a) spectroscopy, which can help us to estimate α and the inclination angle factor sini; b) interferometry, which provides information related to the stellar geometrical deformation and incidentally with α. We briefly discuss these issues in Sects. 6.1 and 6.2.

3.2.2. Stars with polar dimples

Owing to the centripetal acceleration, the fraction of the gravitational force that must be supported by the pressure gradient increases towards the center of the stars. This means that the pressure gradient is more centrally concentrated in a rotating star than a non-rotating one. The rotation tends then to decrease the moment of inertia, which is accompanied by an increase in the density gradient towards the center. The gravitational attraction towards the pole is thus stronger, which produces a flattening of the pole or a dimple. In models of stars with high differential rotations in depth that enables the star to have a significant content of rotational energy, i.e. when the energy fraction τ = K/|W| (K =  rotational kinetic energy; W =  gravitational potential energy) becomes τ ≳ 0.03, the polar region is either strongly flattened or hollower (Bodenheimer 1971; Clement 1974, 1979; Zorec 1986; Jackson et al. 2005). Similar dimples can also be produced by a surface angular velocity accelerated enough towards the polar latitudes. To illustrate this effect, we adopt the following angular velocity profile, which is sometimes used in the literature to model internal conservative (cylindrical) laws (Jackson et al. 2005) for stars with high rotational energies (0.02 ≲ τ ≲ 0.22) ${\mathrm{\Omega }}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{o}}}{\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{\left[}{\mathit{R}}_{\mathrm{s}}\mathrm{\right(}\mathit{\theta }\mathrm{)}\mathit{/}{\mathit{R}}_{\mathrm{e}}{\mathrm{]}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}\mathit{,}$(46)where R_s(θ) is the equation of the stellar surface and β is a free parameter. A relation similar to Eq. (46) may suit the rotational profiles shown in Fig. 2a if the angular velocity is accelerated from the equator to the pole. Figure 5 shows the shape of model stars with a surface angular velocity profile given by expression (46) for β = 5.0 and several values of ratio η_o in the equator calculated with Eq. (39). It can be shown that using the Roche approximation and the rotational law in Eq. (46), the ratio of equatorial to polar radii is then given by $\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathit{\eta }}_{\mathrm{o}}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{.}$(47)Writing ${\mathit{R}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}$, the condition for a flattened/hollowed polar region to appear is dz/dx ≥ 0 at x = 0. The condition to obtain a hollowed polar region is then given by ${\mathit{\eta }}_{\mathrm{o}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }{\mathrm{)}}^{\mathrm{2}}\mathrm{\ge }{\left[\mathrm{1}\mathrm{+}{\mathit{\eta }}_{\mathrm{o}}\left(\frac{\mathrm{1}\mathrm{+}\mathit{\beta }}{\mathrm{2}}\right)\right]}^{\mathrm{3}}\mathit{.}$(48)Making u = R_e/R_p = 1 + η_o(1 + β)/2, relation (48) becomes u³/(u − 1) ≥ 2(1 + β). Since it is  [u³/(u − 1)] _min for u = 1.5, condition (48) requires that β ≥ 2.375 independently of the value of η_o.

Since the simplest equation to represent the rotationally deformed surfaces shown in Fig. 5 is $\frac{{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\mathrm{=}\mathrm{1}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{e}}}\right)\cos \mathit{\theta ,}$(49)an alternative expression to relation (28) for the surface angular velocities accelerated from the equator to the pole could be ${\mathrm{\Omega }}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{o}}}{\mathrm{1}\mathrm{+}\mathit{\beta }{\left[\mathrm{1}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{e}}}\right)\cos \mathit{\theta }\right]}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }}\mathit{,}$(50)would require us to infer the ratio R_p/R_e from interferometry, which according to Eq. (47) is related to η_o and β, and β derived from spectroscopy. However, before using expressions of the type given by either Eq. (46) or (50), it would perhaps be preferable to await observational confirmation that stellar shapes similar to those depicted in Fig. 5 actually exist.

4. The gravity darkening effect in stars with non-conservative rotation laws

Owing to the radiative equilibrium that prevails in the atmosphere of massive and intermediate-mass stars, the emerging bolometric radiative flux F is proportional to the temperature gradient ∇T in the outer layers of the stellar envelope. Since these layers are also in hydrostatic equilibrium, their pressure gradient ∇P ∝ ∇T is balanced by the surface effective gravity g_eff. We have then the following phenomenological relation that describes the principle of von Zeipel’s relation $\mathit{F}\mathrm{\propto }\mathrm{\nabla }\mathit{T}\mathrm{\equiv }\mathrm{\nabla }\mathit{P}\mathrm{\propto }{g}\mathrm{eff}\mathit{.}$(51)As $\mathit{F}\mathrm{\propto }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}$, where T_eff is the effective temperature at a given point $\overline{\mathit{r}}$ in the stellar surface, it follows the known von Zeipel’s relation (von Zeipel 1924)

${\mathit{T}}_{\mathrm{eff}}\mathrm{=}C\mathrm{\left(}\overline{\mathit{r}}\mathrm{\right)}\mathrm{\times }{g}\begin{array}{}\mathrm{1}\mathit{/}\mathrm{4}\\ \mathrm{eff}\end{array}\mathit{.}$(52)However, according to Poincaré-Wavre’s theorem (Tassoul 1978) only when the rotational law in the external layers is conservative do the constant equipotential, density, and temperature surfaces coincide to consider that $C\mathrm{\left(}\overline{\mathit{r}}\mathrm{\right)}\mathrm{=}\overline{C}\mathrm{=}$ constant. Otherwise, as for rotational laws discussed in Sect. 2.2, $C\mathrm{\left(}\overline{\mathit{r}}\mathrm{\right)}$ is not constant with the colatitude θ. Its detailed expression for the non-conservative case with shellular rotation was obtained by Maeder (1999).

An equivalent relation to Eq. (52) for spectroscopic and interferometric analysis can be given the form ${\mathit{T}}_{\mathrm{eff}}\mathrm{=}\overline{C}\mathrm{\times }{g}\begin{array}{}{\mathrm{fi}}_{\mathrm{GD}}\\ \mathrm{eff}\end{array}\mathit{,}$(53)where fi_GD = 0.25 + δ, so that δ masks the variation in $C\mathrm{\left(}\overline{\mathit{r}}\mathrm{\right)}$ over the stellar surface. The particular notation fi_GD used for the gravity darkening power in relation (53) is introduced to avoid confusions with the β parameter that appears in the internal rotation law (46). Two-dimensional (2D) models of rotating stars and some observations suggest that δ ≲ 0 (cf. Lucy 1967; Lovekin et al. 2006; van Belle et al. 2006; Monnier et al. 2007; Zhao et al. 2009). We insist on δ not being in principle the function of taking into account the change in the stellar surface from radiative to possible convective equilibrium-dominated layers because of the strong change in T_eff with θ, as seems to be understood in Aufdenberg et al. (2006) and Zhao et al. (2009), but mainly the non-conservative nature of the rotational law in the external layers.

5. Models of rotating stars

To study the effect of the centrifugally distorted internal mass distribution on the gravitational potential in the surface of a star, to test the validity of the Roche approximation, it is enough to calculate models of stellar structure where only the global dynamical aspects induced by the rotation are considered. To this end, we differentiate the primary effects produced by the rotation from those induced by the evolution. The primary thermodynamic effects due to the stellar evolution are taken into account by making use of the relations between the pressure and density calculated with one-dimensional models of stellar evolution without rotation. We assume thus that the changes produced on the barotropic relation P = P(ρ) at a given evolutionary stage of a star by the several instabilities and the diffusion of chemical elements unleashed through the stellar evolution by the rotation, have second order effects on the establishment of the dynamical equilibrium of the rotating star. In principle, nothing prevents us from using the barotropic relations derived with models of stellar evolution with shellular rotation, but the results will probably not be more reliable. They all are calculated for low energy ratios τ = K/|W|, which do not correspond to the model inputs regarding the rotational law and energy tested in the present attempt.

Since we are not interested in the precise description of all non-linear time-dependent phenomena associated with the viscosity and internal flows in rotating stars, and because the total energy carried by the meridional circulation is small, our models are axisymmetric, steady state, and circulation free. Owing to these assumptions and Poincaré-Wavre’s theorem (Tassoul 1978), our model stars behave as barotropes. We then adopt internal rotational laws of conservative form, Ω = Ω(ϖ), where ϖ is the distance to the rotation axis. According to the discussion in Sect. 5.2, the conservative laws are expected to produce stronger stellar geometrical deformations for a given amount of rotational kinetic energy than the homologous non-conservative ones.

For conservative rotation laws, the gravitational potential Φ(ϖ,z) and the density distribution ρ(ϖ,z) in the rotating star are simultaneous solutions to the hydrostatic equilibrium equation ${\mathit{\rho }}^{-1}\mathrm{\nabla }\mathit{P}\mathrm{=}\mathrm{\nabla \Phi }\mathrm{+}{\mathit{j}}^{\mathrm{2}}{\mathit{\varpi }}^{-3}\mathit{ê}\mathit{\varpi }$(54)and the Poisson equation $\mathrm{\Delta \Phi }\mathrm{=}\mathrm{4}\mathit{\pi G\rho ,}$(55)where (ϖ,φ,z) are the cylindrical coordinates with z containing the rotation axis, ê_ϖ is the unit vector perpendicular to the z-axis, P is the pressure, and j = Ωϖ² is the specific angular momentum.

Equations (54) and (55) are solved with the adopted complementary barotropic relation $\mathit{P}\mathrm{=}\mathit{P}\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathit{.}$(56)The P = P(ρ) relations used in this work have a two-component polytropic character $\mathit{P}\mathrm{=}\mathit{a}{\mathit{\rho }}^{{\mathit{\gamma }}_{\mathit{a}}}\mathrm{+}\mathit{b}{\mathit{\rho }}^{{\mathit{\gamma }}_{\mathit{b}}}\mathit{,}$(57)where the constants a, b, γ_a, and γ_b were adjusted so as to: a) reproduce the pressure P_c and the density ρ_c in the center of the non-rotating star of given mass and evolutionary stage; b) ensure a continuous distribution of the pressure-density relation at the radius of the stellar core; c) obtain the correct stellar mass at the stellar radius as tabulated by Schaller et al. (1992) for 1-D evolutionary models for the initial metallicity Z = 0.02. The function (57) is continued in the stellar atmosphere by another pressure-density relation calculated by Castelli & Kurucz (2003) for stellar atmospheres as a function of the parameters (T_eff,log g).

The first order effects due to the stellar evolution are thus accounted for by the pressure-density relations at the center of the star and the ∂P/∂ρ gradients. An additional term in relation (57) could also in principle take into account the presence of the convective regions in the envelope induced by the rapid rotation, but this was not done in the present approach. The only rotational effect on the P = P(ρ) relation considered here is by means of the mass-compensation effect (Sackmann 1970), which increases the density ρ_c at the center of the star. To this end, we iterated ρ_c until the nominal stellar mass M was attained. This iteration also implied that the central pressure P_c was changed in accordance.

The gravitational potential Φ(ϖ,z) is obtained by solving the Poisson equation given in Eq. (55) with the cell-method adapted by Clement (1974) to stellar structure calculations. The density distribution ρ(ϖ,z) is derived from the integrated form of Eq. (54) ${\mathrm{\int }}_{{\mathit{\rho }}_{\mathrm{c}}}^{\mathit{\rho }}\frac{\mathrm{d}\mathit{P}}{\mathit{\rho }}\mathrm{=}\mathrm{\Phi }\mathrm{\left(}\mathit{\varpi ,z}\mathrm{\right)}\mathrm{-}{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\varpi }}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\varpi }\mathrm{\right)}\mathit{\varpi }\mathrm{d}\mathit{\varpi ,}$(58)where Φ_c is the gravitational potential at the stellar center. Given a rotational law Ω(ϖ) and a barotropic relation in Eq. (56), the common solution to Eqs. (54) and (55) is performed over the entire space. The iteration of Φ and ρ is stopped when the highest density difference in the (ϖ,z)-space is max(δρ/ρ) ≲ 10^-6. In our iterations, the virial relation δ =  [2(K + U) − W]/|W| = 0, where K = kinetic energy, U = internal energy, and W = total gravitational potential energy, is verified to better than δ ≈ 2 × 10^-4 in the ZAMS models, and δ ≈ 6 × 10^-3 by the TAMS models. Since in the framework of conservative rotational laws, the surfaces of constant pressure, density, and total potential are parallel, the rotationally distorted shape of our models is defined by the total equipotential surface that contains the polar “photospheric” radius R_p. This radius is identified by the layer whose density satisfies the model-atmosphere relation τ_Ross(ρ) = 2/3 in the stellar atmosphere models of Castelli & Kurucz (2003). The local effective temperature in the pole needs to satisfy also the gravity darkening effect. We modified accordingly the effective temperature given by Schaller et al. (1992) for the given mass M using von Zeipel’s approximation (von Zeipel 1924). The transformation to the rotation-dependent effective temperature was performed following the procedure given in Frémat et al. (2005).

5.1. Rigid rotation

The simplest conservative rotational law is that of rigid rotation. The models of stars with rigid rotation were calculated assuming that at each evolutionary phase they experience instantaneous total redistribution of their internal angular momentum. The characteristics of these models are given in Table 2. In this table, t/t_MS is the fractional age of the star, where t_MS is the time that a non-rotating star of mass M spends on the main sequence, Ω_cr is the critical angular velocity, Ω/Ω_cr represents the angular velocity for which the model was calculated, ρ_c is the core density of the rotating object, R_e/R_⊙ is the equatorial radius of the model star in solar units and R_e/R_p the equatorial-to-polar radii ratio, V_eq is the equatorial linear velocity in km s^-1, J/M is the total specific angular momentum, $\mathit{\eta }\mathrm{=}{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}\mathit{/}\mathit{GM}$ is the ratio of centrifugal to the gravitational acceleration in the equator, and K/|W| is the ratio of the kinetic rotational energy (K) to the absolute value of the gravitational potential energy (W). Some of these models are shown in Fig. 6, where in all cases the iso-density surfaces are for the same ρ/ρ_c density ratios. Similar models were also calculated for other masses and age ratios t/t_MS. In spite of the simple approach used to calculate them, the radius ratios compare very closely with those calculated by Ekström et al. (2008) for the same K/|W| energy ratios. We note, however, that the models obtained by these authors cannot be compared directly with ours because the distribution of the internal rotational velocity is not the same. Starting from a quasi-rigid rotation in the ZAMS, Ekström et al. (2008) accounted for a consistent evolution of the angular momentum distribution inside the star throughout the calculated stellar evolution span.

5.2. Shellular differential rotation

Zahn et al. (2010) studied the shapes of stars with internal shellular differential rotation and concluded that in 7   M_⊙ stars the radius ratio R_e/R_p at critical equatorial rotation can be enlarged from 1% to 4%, depending on the internal rotational energy content and the evolutionary phase. However, it can easily be seen that when two rotation laws are described with the same function and both imply the same central to surface angular velocity ratios Ω_center/Ω_surface, but one is shellular and the other cylindrical, the law that is shellular may have a weaker effect on the global internal density distribution than the cylindrical one. For simplicity, we change the independent variable P (pressure) in the shellular law by the radius r. The comparison is depicted in Fig. 8: the specific centrifugal force acting at a point p in the stellar interior is $\begin{array}{ccc}{\mathit{F}}_{\mathrm{c}}^{\mathrm{shellular}}& \mathrm{=}& {\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\varpi }\\ {\mathit{F}}_{\mathrm{c}}^{\mathrm{cylindrical}}& \mathrm{=}& \end{array}$(59)for shellular and cylindrical angular velocity distribution laws, respectively (ϖ = rsinθ is the distance to the rotation axis). At the equator, if we assume that Ω^cylindrical(ϖ,π/2) = Ω^shellular(r,π/2), for a given point p(z,ϖ), it will be $\begin{array}{ccc}\mathrm{\Omega }\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{cylindrical}}& \mathit{>}& {\mathrm{\Omega }}^{\mathrm{shellular}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{,}\\ {\mathit{F}}_{\mathrm{c}}^{\mathrm{cylindrical}}\mathrm{\left(}\mathit{p}\mathrm{\right)}& \mathit{>}& \end{array}$(60)This holds in particular for the points in the stellar surface. However, since it has been shown in Sect. 3.2 that the external geometrical deformation of the star is mainly a function of the angular velocity law on the surface, cylindrical rotation laws will carry stronger geometrical deformations than shellular ones for a similar total kinetic energy. Moreover, in Sect. 2.2 we have anticipated that strong cylindrical components of the angular velocity distributions in the external layers probably exist. It is then important to determine the order of magnitude of the quadrupole factor γ in rotating stars using conservative internal rotational laws. To this end, in the next section we calculate two-dimensional models of rotating stars with conservative internal rotational laws.

5.3. Cylindrical differential rotation

The models with differential rotation were calculated using the same barotropic relations P = P(ρ) as for the solid-body rotation. We adopted the internal angular velocity distributed as $\mathrm{\Omega }\mathrm{\left(}\mathit{\varpi }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{o}}}{\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{(}\mathit{\varpi }\mathit{/}{\mathit{R}}_{\mathrm{e}}{\mathrm{)}}^{\mathit{n}}}\mathit{,}$(61)where β > 0 is a free parameter and Ω_o is the angular velocity at the axis of rotation. For simplicity, we considered n = 2. Nevertheless, it can be shown that all models with n ≤ 2 obey the Solberg-Høiland stability criterion, and that at whatever distance ϖ the energy ratio τ(ϖ) is K(ϖ)/|W(ϖ)| ≲ 0.138. This ensures that no region in the stellar interior is unstable to the secular instability that can carry it to a three-axial Jacoby ellipsoidal configuration.

In Table 3, we indicate the characteristics of the calculated models. The entries have the same meaning as in Table 2. The higher K/|W| values given in Table 3 are about the highest we could obtain for the given masses, ages, and β parameters. Figure 7 shows some of these calculated models.

Using the same algorithm, we can obtain models that approach more closely the limit of dynamical stability K/|W| ~ 1/4 if the objects are on the ZAMS and have masses higher than M = 15   M_⊙. However, Clement (1979) in his analysis of the secular stability of differentially rotating objects already pointed out the unstable character of model stars with K/|W| ≳ 0.1.

5.4. The Roche approximation

With the models thus calculated, which give us the shape of a rotationally distorted star, we can test to what degree the gravitational potential Φ_G in the stellar surface may deviate from the simple central-field expression used in the Roche approximation.

Hubbard et al. (1975) showed that the gravitational potential of rotating centrally condensed objects can be given by a multipol expansion ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\left[\mathrm{1}\mathrm{-}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\left(\frac{{\mathit{R}}_{\mathrm{o}}}{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\right)}^{\mathrm{2}\mathit{n}}{\mathit{J}}_{\mathrm{2}\mathit{n}}{\mathit{P}}_{\mathrm{2}\mathit{n}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\right]\mathit{,}$(62)where R_o is the radius of the rotationally undistorted star, J_2n are the zonal harmonic coefficients, and P_2n(cosθ) are the Legendre polynomials. To derive relation (62), Hubbard et al. (1975) assumed that the density can be expanded in Legendre polynomials. The J_2n coefficients are then obtained as integrals over the volume of the axisymmetric object which generally produce values so that J₂ > 0, J₄ < 0, J₆ > 0, J₈ < 0, and so on. Since we are not interested in developing a detailed theory of the gravitational potential of rotating stars, but only to test the validity of the Roche approximation, we do not calculate the volume integrals over the rotationally deformed star, but use the harmonic coefficients coefficients as mere parameters giving a quantitative indication of the deviation from the central-field form of the gravitational potential. We determine the quantities J_2n in expression (62) by simply fitting the gravitation potential predicted by the 2D models of rotating stars calculated in Sects. 5.1 and 5.3 using a least squares method. In this algebraic estimate of J_2n, their order of magnitude are preserved, but not always the expected sign.

In relation (62), we retain only the quadrupole moment to estimate the effects in rigid rotators ${\mathrm{\Phi }}_{\mathrm{G}}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathrm{1}\mathrm{-}{\left[\frac{{\mathit{R}}_{\mathrm{o}}}{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\right]}^{\mathrm{2}}{\mathit{J}}_{\mathrm{2}}{\mathit{P}}_{\mathrm{2}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{.}$(63)The shape of the surface thus becomes $\frac{{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}\mathit{\gamma }{\mathrm{(}{\frac{{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{2}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{)}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}_{\mathrm{o}}\left[\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}}^{\mathrm{)}}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }\right]}\mathit{,}$(64)where we have introduced the notation $\mathit{\gamma }\mathrm{=}{\mathit{J}}_{\mathrm{2}}{\left(\frac{{\mathit{R}}_{\mathrm{o}}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{2}}\mathrm{·}$(65)The largest difference between the estimate given by Eq. (64) to the ratio R_s(θ)/Re and that produced by the Roche approximation in Eq. (44), is expected for θ = 0. Thus according to (64) we have $\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{\eta }}_{\mathrm{o}}\mathrm{+}\mathit{\gamma }\mathrm{)}}{\mathrm{1}\mathrm{-}\mathit{\gamma }{\left(\frac{{\mathit{R}}_{\mathrm{e}}}{{\mathit{R}}_{\mathrm{p}}}\right)}^{\mathrm{2}}}\mathrm{·}$(66)Since it is always true that R_e/R_p ≥ 1, relation (66) implies that for whatever 0 ≤ η_o ≤ 1 and γ ≠ 0 the equator to polar radii ratio is slightly larger than obtained from the sheer Roche approximation.

In Table 4, we indicate the quadrupole factors γ for 5 M_⊙ and 15 M_⊙ stars in the ZAMS and TAMS evolutionary phases obtained by fitting relation (64) with surface gravitational potential obtained with the model calculation. In this table, we see that the smaller J₂ the more evolved is the star, simply because the star is more centrally condensed such that the central-field approximation for Φ_G holds better. Zahn et al. (2010) calculated the same factors for rigid and shellular differential rotators. Reducing their J₂ estimates by (R_p/R_e)² to be able to approach the ratio (R_o/R_e)² more closely than using the average stellar radius for R_o, we see that for the ZAMS and TAMS epochs of rigid critical rotators we obtain the same values for the harmonic coefficient J₂ than Zahn et al. (2010). Having then γ ≲ 0.004, from (44) we see that neglecting the quadrupole term in (64), in rigid rotators at critical rotation the R_e/R_p can be underestimated by less than 2%. For shellular rotators, Zahn et al. (2010) found that at critical equatorial rotation the flattening changes by 4% because Ω_center/Ω_surface = 4 in the ZAMS and is smaller near the TAMS.

Finally, we note that within approximation (64), the actual ratio of the centrifugal to gravity acceleration at the equator should be ${\mathit{\eta }}_{\mathit{J}}\mathrm{=}\frac{{\mathit{\eta }}_{\mathrm{o}}}{\mathrm{1}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\mathit{\gamma }}\mathrm{·}$(67)In stars with conservative differential rotation laws of the type given in Eq. (61) that imply rotational energies τ = K/|W| > τ_cr, the deformations of the internal mass distribution and the deviations to the central-field gravitational potential are expected to be stronger than in stars with rigid rotation. To test the use of the Roche approximation in these cases, we calculated the zonal harmonic coefficients γ_2n = J_2n(R_o/R_e)² from n = 1 to 4. Table 5 lists the zonal multipolar coefficients γ_2n obtained by fitting relation (62) with the calculated gravitational potential for several values of the coefficient β and angular velocity ratio Ω_o/Ω_cr, where Ω_o is the angular velocity at the rotation axis, and Ω_cr is the critical angular velocity of an homologous star (same mass and evolutionary stage) behaving as a rigid rotator. We insist on the algebraic nature of the γ_2n coefficients, since they were derived as fitting parameters to model-calculated gravitational potential. In this Table, we also give the fractional deviations of the radius vectors ΔR(θ)/R(θ) =  [R(θ)_model − R(θ)_γ2n] /R(θ)_model × 10⁴ for different colatitude angles θ, where R(θ)_model is obtained with the models calculated in Sect. 5.3, while the R(θ)_γ2n values refer to those estimated with the total potential that includes the gravitational potential given by Eq. (62). A quick inspection of Table 5 shows that in most cases a fairly precise description of the stellar surface can be obtained by considering only γ₂ and γ₄ in the gravitational potential given by Eq. (62). Nevertheless, in all cases the Roche approximation can be valid to better than 5% for ratios Ω_center/Ω_surface ≃ 4. The approximation works better the more evolved the star is and the higher the rotational energy parameter τ.

6. Attainable stellar parameters

According to Vogt-Russell’s theorem, “the complete structure of a star in hydrostatic and thermal equilibrium is uniquely determined by the total mass M and the run of chemical composition throughout the star, provided the structure equations are function of local parameters” (Cox & Giuli 1968). Then, for a given chemical composition, denoted here in a generic way by Q, it follows that parameters that conventionally specify the physical properties of a star are the luminosity L, the mass M, and the radius R (or the stellar age t).

An object in permanent rotation can still be considered to be in hydrostatic and thermal equilibrium. However, owing to the numerous effects the rotation induces, the number of parameters needed to characterize an object can be high: mass, total luminosity, age, total angular momentum, internal distribution of the angular momentum, dimensions or stellar geometry, inclination angle, degree of non-uniform distribution of the temperature in the surface, mixing of chemical elements, and degree of their non-uniform internal and surface distribution, etc. This number is certainly higher than we can tackle and/or determine from observations. We consider here only parameters we can actually determine with the diagnostic tools we presently have at our disposal. However, since many of these quantities are model-dependent, the resulting characterization of rotating stars would probably be limited or incomplete.

We may classify the accessible stellar parameters into three classes: 1) those related to the stellar internal structure; 2) parameters relevant to the structure of the atmosphere over the observed stellar hemisphere; 3) quantities related to the apparent stellar geometry.

Models of stellar structure: the existing models of stars with rotation start their MS evolution with low rotational energies in the ZAMS, i.e. $\mathit{\tau }\mathrm{=}\mathit{K}\mathit{/}\mathrm{\left|}\mathit{W}\mathrm{\right|}\mathrm{\le }{\mathit{\tau }}_{\mathrm{cr}}^{\mathrm{rig}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathit{,}$(68)where K is the rotational kinetic energy, W is the gravitational potential energy, and ${\mathit{\tau }}_{\mathrm{cr}}^{\mathrm{rig}}$ is the energy ratio for a rigid rotator at a critical equatorial rotation in the ZAMS. For stars with masses ranging from 3   M_⊙ to some 60   M_⊙, it is on average ${\mathit{\tau }}_{\mathrm{cr}}^{\mathrm{rig}}\lesssim \mathrm{0.015}$ (Zorec et al. 1988a). Fom Table 2, we see that τ_ZAMS is a function of M and V_eq. For rigid rotators in the ZAMS, V_eq and Ω/Ω_c are equivalent parameters.

Owing to mass and angular momentum losses, the energy ratio τ becomes a function of the time, so that $\mathit{\tau }\mathrm{\left(}\mathit{M,t}\mathrm{\right)}\mathit{<}{\mathit{\tau }}_{\mathrm{ZAMS}}^{\mathrm{rig}}\mathrm{\left(}\mathit{M}\mathrm{\right)}$. On account of possible redistribution phenomena of angular momentum that could take place in the pre-MS evolutionary phases, we should characterize the stellar initial rotational law with at least one single quantity. This single parameter should also correspond to the expected value of the kinetic energy ratio when ${\mathit{\tau }}_{\mathrm{ZAMS}}\mathit{>}{\mathit{\tau }}_{\mathrm{ZAMS;cr}}^{\mathrm{rig}}$. However, existing models of rotating stars that provide us today with evolutionary paths to infer the stellar masses and ages, are calculated with initial rigid rotations in the ZAMS. In this case, τ, Ω/Ω_c, and V_eq can be assumed to carry the same information. We must keep in mind that the masses and ages we infer from these models are depend on the assumption given in Eq. (68). The models provide us with quantities averaged over the stellar surface (cf. Meynet & Maeder 2000; Ekström et al. 2008), namely surface-averaged bolometric luminosity  ⟨ L ⟩ , effective temperature  ⟨ T_eff ⟩ , surface effective gravity  ⟨ log g_eff ⟩  (occasionally), and equatorial linear rotational velocity V_eq$\begin{array}{ccc}\mathrm{⟨}\mathit{L}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}\mathit{L}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\tau }}_{\mathrm{ZAMS}}\mathrm{\right)}\mathrm{⟩}\mathit{,}\\ \mathrm{⟨}{\mathit{T}}_{\mathrm{eff}}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}{\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\tau }}_{\mathrm{ZAMS}}\mathrm{\right)}\mathrm{⟩}\mathit{,}\\ \mathrm{⟨}\log {\mathit{g}}_{\mathrm{eff}}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}\log {\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\tau }}_{\mathrm{ZAMS}}\mathrm{\right)}\mathrm{⟩}\mathit{,}\\ {\mathit{V}}_{\mathrm{eq}}& \mathrm{=}& {\mathit{V}}_{\mathrm{eq}}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\tau }}_{\mathrm{ZAMS}}\mathrm{\right)}\mathit{.}\end{array}$(69)From Maeder & Meynet (2000), we can see that because of rotation, stars behave as though they have lower effective mass than in reality, which mainly affects the emitted bolometric luminosity. According to this mass-compensation effect (Sackmann 1970) and the core fuelling with hydrogen by the meridional circulation, rotating stars have longer MS life-spans. Owing to the change in the opacity of the envelope by the diffusion of chemical elements from the core, mainly He, stars become more luminous than their mass-homologous rotationless objects. Evolutionary timescales in the post-MS phases are also sensitively changed. If stars were actually neat differential rotators, i.e. with rotational energies ${\mathit{\tau }}_{\mathrm{ZAMS}}\mathit{>}{\mathit{\tau }}_{\mathrm{cr}}^{\mathrm{rig}}$ (Bodenheimer 1971; Clement 1979), the evolutionary tracks needed to interpolate stellar masses would also differ from those existing today.

Models of stellar atmospheres: high resolution spectra and spectral energy distributions of rapidly rotating stars are both i-dependent (i is the inclination angle). They can be modeled by taking into account the rotationally-induced stellar geometrical deformation, but mainly the concomitant latitude-dependent surface gravity and temperature distribution: gravity darkening (Frémat et al. 2005; Lovekin et al. 2006; Gillich et al. 2008).

Since a priori we do not have any indication on what may be the function describing the surface angular velocity of a star, in Sect. 3.1 we adopted a simple Maunder-like expression given by Eq. (28), i.e., Ω(θ) = Ω(Ω_o,α,θ), where Ω_o is the angular velocity of the equator, α is the surface differential rotation parameter, and θ is the colatitude angle. Instead of using Ω_o we prefer to use the force ratio parameter η_o defined in Eq. (41). Since α ≠ 0, the iso-radial velocity lines responsible for the Doppler shifts and causing the rotationally broadening of spectral lines are no longer straight lines but curves, so that the true velocity in Vsini does not necessarily correspond to the equatorial linear velocity (Zorec et al. 1988b, 2004). To this adds the uncertainty in the Vsini parameter related to the less effective contribution to the line broadening by the equatorial regions affected by the gravitational darkening (Stoeckley 1968b; Townsend et al. 2004; Frémat et al. 2005).

For the colatitude-dependent distribution of the effective temperature, we concluded in Sect. 4 that ${\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\overline{C}\mathrm{\times }{g}\begin{array}{}{\mathrm{fi}}_{\mathrm{GD}}\\ \mathrm{eff}\end{array}$, where $\overline{C}$ can be assumed to be caused by the average stellar geometrical deformation and where fi_GD = 0.25 + δ has a free parameter δ that should be δ ≤ 0 according to what is suggested in the literature (cf. Lucy 1967; Lovekin et al. 2006; van Belle et al. 2006; Monnier et al. 2007; Zhao et al. 2009). For models that begin evolving in the ZAMS as rigid rotators, τ and η_o are also synonymous. The apparent, hemisphere-dependent bolometric luminosity, effective temperature, effective gravity and rotation parameter are then functions of at least the seven unknowns of mass M, age t, initial chemical composition Q, equatorial rotation parameter η_o, surface differential parameter α, power fi_GD in the von Zeipel relation, and the inclination angle i$\begin{array}{ccc}{\mathit{L}}^{\mathrm{app}}& \mathrm{=}& \mathit{L}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathit{,i}\mathrm{\right)}\mathit{,}\\ {\mathit{T}}_{\mathrm{eff}}^{\mathrm{app}}& \mathrm{=}& {\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathit{,i}\mathrm{\right)}\mathit{,}\\ \log {\mathit{g}}^{\mathrm{app}}& \mathrm{=}& \log \mathit{g}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathit{,i}\mathrm{\right)}\mathit{,}\\ \mathrm{(}\mathit{V}\sin \mathit{i}{\mathrm{)}}^{\mathrm{app}}& \mathrm{=}& \mathit{V}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathrm{\right)}\mathrm{\times }\sin \mathit{i}\mathit{.}\end{array}$(70)We also include the dependence on fi_GD in log g^app and (Vsini)^app, since their prediction is based on spectral lines produced by gravity-darkened models. The above seven-fold dependence is valid also for the equivalent widths (W_λ) of spectral lines, the line FWHM (full width at half maximum), the zeroes q_n of the Fourier transforms of spectral lines, as well as the energy distributions and spectrophotometric parameters such as (λ₁,D) of the BCD system, which describe the Balmer discontinuity and are highly sensitive to ($\log {\mathit{g}}^{\mathrm{app}}\mathit{,}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{app}}$) (Zorec et al. 2009). Referring to all these observable quantities with a generic letter O_λ, they are formally $O_{\mathit{\lambda }}\mathrm{=}{O}_{\mathit{\lambda }}\mathrm{\left(}\mathit{M,t,Q,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathit{,i}\mathrm{\right)}\mathit{.}$(71)We note that for stars with a given initial chemical composition Q and a prescription for their internal rotation, the parametric couples (⟨ L ⟩ , ⟨ T_eff ⟩), and (M,t) are equivalent quantities. From observations, we derive L^app, ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{app}}$, and log g^app, which are i-dependent, while the models of stellar evolution provide us with  ⟨ L ⟩ ,  ⟨ T_eff ⟩ , and sometimes  ⟨ log g ⟩ . This means that M and t can be estimated by iteration, where we also iterate simultaneously the inclination angle i (Zorec et al. 2005).

Stellar geometry: since it is possible to use the Roche approximation to describe the rotationally induced stellar deformation, from Eq. (37) it comes that the inclination angle-dependent shape of rapidly rotating star depends only on η_o, α and i. If we are not interested in carrying out imaging of fast rotating objects, interferometry can provide us with data that to a first approximation do not depend on M, t, Q, and fi_GD, i.e. the apparent shape of the star (y,x), where x is the coordinate measured along the equatorial line projected onto the sky, and y is: $\mathit{y}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathit{R}\mathrm{\left(}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha ,i,x}\mathrm{\right)}}{{\mathit{R}}_{\mathrm{e}}}\mathit{,}$(72)so that the apparent polar-to-equatorial radius ratios are y(x = 0) = R_p(i)/R_e and y(x = 1) = 0.0. We examine this in greater detail in Sect. 6.2.

If models are calculated for a given metallicity Z, which stands for the initial generic quantity Q used above, we conclude that a massive and intermediate-mass fast rotating star with surface differential rotation can be described as a function of six parameters $\begin{array}{ccc}& \mathit{M,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,}\mathit{t,}\mathit{\alpha ,}{\mathrm{fi}}_{\mathrm{GD}}\mathit{,}\mathrm{and}\mathit{i,}& \end{array}$(73)which we should be able to derive by interpreting the observed data. We briefly review some of the data and techniques that can help us to determine some of these parameters, in particular η_o, α, and i.

6.1. Information on the differential rotation from spectral lines

As stated in Sect. 4, the distribution of the effective temperature in the stellar surface deviates from von Zeipel’s prescription, i.e. ${\mathit{T}}_{\mathrm{eff}}\mathrm{\propto }{\mathit{g}}_{\mathrm{eff}}^{{\mathrm{fi}}_{\mathrm{GD}}}$ with fi_GD ≠ 1/4 if α ≠ 0. However, in what follows of the present discussion we shall assume the validity of the original von Zeipel relation, but use a modified expression of the effective gravity g_eff according to the latitudinal dependence of the surface differential angular velocity Ω(η,α,θ) given by Eq. (28).

In this section, we report first estimates of the effects produced by surface rotational laws Ω(η,α,θ) on some spectral lines and envision with them the way we can obtain estimates of α and the inclination angle i from spectroscopy .

We modified our FASTROT code (Frémat et al. 2005), initially written for a uniform angular velocity on the stellar surface, to take into account an external Maunder-like rotational profile and the related changes in the gravity and the effective temperature as a function of the colatitude θ in the calculation of spectral line profiles. In Fig. 9a, we show the Si ii 6347 line profiles obtained with the new FASTROT program, which are produced by MS stars with T_eff = 19   500 K and log g = 4.0 at rest, whose rotation is characterized by several ratios of the angular velocity in the equator Ω_o/Ω_c without surface differential rotation (α = 0). In Fig. 9b, we show the Si ii 6347 line profiles for stars with the same (T_eff,log g) rest parameters as before, but now for Ω_o/Ω_c = 0.9 and several degrees of differential rotation (α ≠ 0). All line profiles in Figs. 9a and 9b were obtained with gravity-darkened models seen at the inclination angle is i = 45°. We note that the kinematic effect on the line broadening with surface differential rotation is determined by curves of constant Doppler shift, which are not straight lines as in the solid body rotation, but curves whose shapes depend on the value of α and the inclination angle i (Domiciano de Souza et al. 2004b). Other glimpses into the sensitivity of different spectral lines in early-type rapid rotators to the effects carried by a surface rotation law Ω(η,α,θ) can be seen in Zorec et al. (2007).

The Fourier transform can give a quantitative description of the effects induced on the spectral lines by the latitudinal differential rotation and perhaps on the associated gravity darkening effect. This information can be carried by the ratio q₂/q₁ of the second (q₂) to the first zero (q₁) of the Fourier lobes, as suggested by Reiners & Schmitt (2004, and references therein). In Fig. 10a, we see the ratio q₂/q₁ determined by the Fourier analysis of the lines shown in Fig. 9b, as a function of the ratio Ω_o/Ω_c for stars without surface differential rotation (α = 0). We see in this figure that the ratio q₂/q₁ becomes sensitive to the gravity darkening only when Ω_o/Ω_c approaches a critical value. In Fig. 10b, we show the ratios q₂/q₁ for the profiles in Fig. 9b, i.e. Ω_o/Ω_c = 0.9 and α ≠ 0, where the corresponding gravity darkening is also taken into account. As in Fig. 9, the inclination in Fig. 10 is i = 45°. From the results shown in Fig. 10b, we see that the differential rotation in fast rotating early-type star produces fairly large effects on the ratio q₂/q₁, which are then expected to be measurable. It is also expected that a series of relations of the type ${\mathit{q}}_{\mathrm{2}}\mathit{/}{\mathit{q}}_{\mathrm{1}}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{\alpha ,i}\mathrm{\right)}$(74)can be obtained, whose functional dependence on α and i differ from line to line. This is because the dominant excitation mechanism of spectral lines is not the same for all of them, and that they can react more or less in different ways, or be more or less sensitive to the non-uniformity of the surface effective temperature and gravity. This precludes our being able to define a prototype line-broadening function that could be used for all spectral lines analyzed in a spectrum, which will be studied in detail in a following paper.

We note that for this type of study, only spectral lines that do not undergo strong Stark broadening can be used. Accordingly, in early-type stars the He i 4471 and Mg ii 4481 lines are more frequently used to determine the Vsini parameter (cf. Slettebak 1982; Chauville et al. 2001; Frémat et al. 2005). To test the sensitivity of these lines to the effects related to the gravitational darkening, we calculated the line profiles of He i 4471, Mg ii 4481, and Si ii 6371, and obtained the diagrams shown in Fig. 11a, where for each studied line the inclination of the rotation axis is presented against the equivalent width and the FWHM (full width at half maximum). For the calculations shown in Fig. 11, we set the differential rotation parameter to α = 0. However, similar behaviors occur when α ≠ 0. The curves in Fig. 11 show the differentiated sensitivity of the He i, Mg ii, and Si ii lines to the effects carried by the rotation.

Using both the measured equivalent widths and FWHM of lines, and the diagrams of the type shown in Fig. 11a, we obtain intersections as shown in Fig. 11b that determine estimates of Ω/Ω_c and i. From Eq. (74), we can then infer a value of α and use it to obtain a new series of curves similar to those in Fig. 11a. This iteration can be followed until convergence for Ω/Ω_c, i, and α is attained. To infer a first estimate of uncertainties related to this method, we note that the curves of Fig. 11 can be reasonably fitted with relations such as $\mathit{i}\mathrm{=}\mathit{A}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{W}}^{\mathit{B}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}\mathrm{+}\mathit{C}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{,}$(75)where ω = Ω/Ω_c. It then follows that uncertainties δW/W in the equivalent widths carry errors in the inclination angle i which are on the order of: $\frac{\mathit{\delta i}}{\mathit{i}}\mathrm{\simeq }\mathit{B}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\left[\mathrm{1}\mathrm{-}\frac{\mathit{C}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathit{i}}\right]\frac{\mathit{\delta W}}{\mathit{W}}\mathrm{·}$(76)For an order of magnitude estimate, we consider Achernar, for which measurements of the apparent radii ratio R_e/R_p were obtained (Domiciano de Souza et al. 2003; Kervella & Domiciano de Souza 2006) and another analysis of the inclination angle was performed (Vinicius et al. 2006). Assuming for this star that η_o ~ 0.8−0.9 and i ≲ 75°, from the He i4471 and Mg ii4481 we can infer, respectively, that B₄₄₇₁ = 3.3 and C₄₄₇₁ = 11.6, B₄₄₈₁ = 2.5 and C₄₄₈₁ = 13.7, and δW/W ~ 0.03 that lead to δi ~ 5°, which is an upper limit, since as seen in Fig. 11b, the Si ii lines reduce this uncertainty.

A simple polynomial fit of α as a function of the ratio R_q = q₂/q₁ shown in Fig. 10 leads to uncertainties of δα ≃ 0.15 for R_q = 1.65 and δR_q ≃ 0.1, and smaller values of δα for other ratios R_q. However, the determination of the parameters Ω/Ω_c (or η_o) and i, and then α, can be greatly improved using curves similar to Fig. 11 calculated for many other spectral lines of different excitation character, i.e., collision-dominated or photoionization-dominated, which enable then to infer the inclination angle factor sini more precisely (Hutchings & Stoeckley 1977; Ruusalepp 1989). A future paper will be devoted to the discussion of this method and the related uncertainties.

6.2. Information on the differential rotation from interferometry

Interferometry can help provide different types of information on rotating stars: a) data related to the geometrical shape of stars; b) independent estimates of the differential rotation parameter α and the inclination i; c) imaging of the non-uniform hemisphere of a rotationally deformed star. In this introductory paper, we do not address the imaging techniques.

a) Stellar geometry. The main goal is to measure the geometrical deformations carried by the stellar rotation, and in particular those added by the surface differential rotation described in Sect. 3.2.1. Using the integrated fluxes in different spectral bands, the geometrical shape of the sources can be determined from the visibility measurements using the widest possible coverage of the plane (u,v) (van Belle et al. 2001; Domiciano de Souza et al. 2003; Aufdenberg et al. 2006; van Belle et al. 2006; McAlister et al. 2005; Kervella & Domiciano de Souza 2006; Monnier et al. 2007; Zhao et al. 2009). The resulting description of the geometry of the source is however i-inclination angle dependent.

As in Eq. (72), we write as y(x,i) the coordinates measured perpendicular to the x axis drawn in the direction of the stellar equator projected onto the sky. The y(x,i) ordinates determined from the interferometry are then related to the true shape of the star by $\begin{array}{ccc}\mathit{y}\mathrm{\left(}\mathit{i,x}\mathrm{\right)}& \mathrm{=}& \mathrm{\left[}{\mathit{R}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\right(}\mathit{\theta ,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha }\mathrm{)}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{-}{\mathit{x}}^{\mathrm{2}}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{2}}\cos \mathit{i}\\ & & \mathrm{+}{\mathit{R}}_{\mathrm{s}}\mathrm{\left(}\mathit{\theta ,}{\mathit{\eta }}_{\mathrm{o}}\mathit{,\alpha }\mathrm{\right)}\cos \mathit{\theta }\sin \mathit{i,}\end{array}$(77)where R_s(θ,η_o,α) is the equation of the stellar surface normalized to R_e given by Eq. (39) and θ is the colatitude angle measured in the stellar reference system. Equations (77) and (39) are two independent equations that for a given pair of prameters (i,α) enable us to obtain an x = x(η,θ) relation and determine the sought-after parameter η. This leads us to the final objective, i.e., to derive the true radii ratio R_p/R_e = R_s(θ = 0)/R_e of the star from the apparent stellar shape y(i,x)/R_e. To achieve this purpose, the inclination angle i can be inferred using several iteration and modeling methods: 1) approximate ones such as those tested by Vinicius et al. (2006) and Carciofi et al. (2008); 2) the spectroscopic inferences described in Sect. 6.1; 3) spectrophotometric constraints imposed by the (λ₁,D) parameters of the Balmer discontinuity (Zorec et al. 2005); 4) far and near-UV spectrophotometry (Frémat et al. 2002, 2005); 5) differential interferometry as described below (see item b).

A first insight into the uncertainties incurred when using the VINCI/VLTI interferometric data can be inferred by assuming that the apparent star corresponds to the projection of an axisymmetric ellipsoid. Equation (77) then becomes (Kanaan et al. 2008) ${\mathit{\rho }}_{\mathrm{ap}}\mathrm{=}{\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}_{\mathrm{tr}}^{\mathrm{2}}\mathrm{)}{\sin }^{\mathrm{2}}{\mathit{i}}^{\mathrm{]}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(78)where we have ρ_ap = (R_p/R_e)_apparent and ρ_tr = (R_p/R_e)_true. For the VINCI/VLTI data on Achernar analyzed by Kervella & Domiciano de Souza (2006), we infer that δρ_ap ≃ 0.02. Using a Monte Carlo simulation for the propagation of errors, we obtain δρ_tr ≃ 0.03, which from Table 1 implies that the differential rotation parameter α can be estimated with an uncertainty of δα ≲ 0.15.

Improved estimates of ρ_ap are expected from the VEGA/CHARA instrument (Mourard et al. 2009). To determine the lowest stellar fattenings that we can reliably measure and the related uncertainties, the important quantity we need to estimate is the visibility ratio $\frac{\mathrm{d}{\mathit{V}}^{\mathrm{2}}}{\mathrm{⟨}{\mathit{V}}^{\mathrm{2}}\mathrm{⟩}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{-}{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{V}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}$(79)where ${\mathit{V}}_{\mathrm{p}}^{\mathrm{2}}$ is the visibility in the polar direction of the apparent shape of the star and ${\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}$ is the visibility in its equatorial direction. Mourard et al. (2009) demonstrated that highly reliable measurements with the VEGA/CHARA instrument can be performed of the angular diameter when the relative uncertainty in the squared visibility curve is σ_V²/V² ≲ 0.02. We then expect that visibility curves obtained in directions towards the apparent stellar polar and equatorial directions, respectively, be able to be differentiated and thus a difference between R_e and R_p be found as soon as dV²/ ⟨ V² ⟩  ≳ 0.02. We then tested the ratio in Eq. (79) for B0, B2, and B8-type objects of luminosity classes V and III with apparent magnitudes V = 5 mag and V = 4 mag. The magnitude V = 5 mag corresponds to an upper conservative limiting magnitude for the VEGA/CHARA interferometer in the medium-resolution spectral mode. The tested B sub-spectral types were chosen to be B0 to include massive objects for which rapid rotators are frequent, B2 spectral type which represent the highest frequency of Be stars, and B8 to represent Bn stars, which are rapid rotators without emission lines and most currently found among late B-type stars.

Using the fundamental parameters given in Zorec & Briot (1991), Divan & Zorec (1982), and Zorec et al. (2009) for the visible absolute magnitude M_V, the absolute bolometric magnitude M_bol, and the effective temperature T_eff, respectively, we inferred the circular angular diameters θ_o of the tested B sub-spectral types associated with their apparent average radius R_o = (R_e + R_p)/2. We calculated the visibilities for two baselines, B = 156 m and B = 300 m, attainable with the VEGA/CHARA instrument at the effective wavelength λ = 5500   Å, assuming that stellar discs have uniform brightness. The visibilities were calculated separately for the polar and equatorial directions and test radius ratios ranging from R_e/R_p = 1.05 to R_e/R_p = 1.25, to account for the lower deformations of rotators characterized by parameters (η_o,α) reported in Table 1. The results obtained are displayed in Table 6. We note that the fundamental parameters given in Table 6 are not for specific stars, but correspond to hypothetical objects having the quoted MK spectral types.

From the results displayed in Table 6 we conclude that provided the objects are observed with large enough baselines, even a modest flattening R_e/R_p ~ 1.05 could be detected at the limiting magnitude V = 5 mag. Brighter objects with slightly larger flattenings are found to have visibility ratios dV²/⟨ V² ⟩  that are easily beyond the uncertainty range of 3σ_V²/V² ~ 0.06. Hence, rotational parameters (η_o,α) can be inferred with the VEGA/CHARA instrument to much better than 2% in rotators brighter than V = 5 mag with rates η_o ≳ 0.4.

For the sake of completeness, we derive in Appendix A an expression for the visibility ratio dV²/⟨ V² ⟩  as a function of the apparent flattening R_e/R_p and the argument x = πθ_o(B/λ), where θ_o = (2/D) [(R_e + R_p)/2] is the average angular diameter, given by $\begin{array}{c}\\ \frac{\mathrm{d}{\mathit{V}}^{\mathrm{2}}}{\mathrm{⟨}{\mathit{V}}^{\mathrm{2}}\mathrm{⟩}}& \mathrm{=}& \mathrm{-}\mathrm{4}\left\{\frac{{\mathit{\delta }}_{\mathrm{ep}}\mathrm{\left[}\mathrm{\Gamma }\mathrm{\right(}\mathit{x}\mathrm{)}\mathit{/}\mathit{V}\mathrm{(}\mathit{x}\mathrm{\left)}\mathrm{\right]}}{\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathrm{ep}}^{\mathrm{2}}\mathrm{\left[}\mathrm{\Gamma }\mathrm{\right(}\mathit{x}\mathrm{)}\mathit{/}\mathit{V}\mathrm{(}\mathit{v}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}\right\}\\ \mathrm{\Gamma }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathrm{d}\mathit{V}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{/}\mathrm{d}\ln \mathit{x}\\ {\mathit{\delta }}_{\mathrm{ep}}& \mathrm{=}& \frac{\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathit{/}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathit{/}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{+}\mathrm{1}}\end{array}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}$(80)which gives a clear insight into the dependence of dV²/⟨ V² ⟩  on R_e/R_p, provided that V_e and V_p are both in the respective first lobes of J₁(x), i.e. x_p ≲ 3.8327 and x_e ≲ 3.8327. In Fig. 12, we compare the ratios dV²/⟨ V² ⟩  given as a function of the apparent flattening R_e/R_p and calculated for several average arguments x (circles). In the same figure, we also superimpose the estimates derived from the approximation given in Eq. (80) (dashed lines). Figure 12 may help us to sharpen the observation strategy. From the figure, it is clear that when R_e ~ R_p, it is dV²/⟨ V² ⟩  → 0 and when for the same baseline one of the visibilities is significantly higher than the other the visibility ratio saturates, dV²/⟨ V² ⟩  → 2, so that the resolution for the determination of high flattenings becomes less obvious.

Finally, we note that for η_o ≃ 0.4 and αt0 the flattening is R_e/R_p ≃ 1.2, which also implies that Ω_oΩ_c ≃ 0.63 and V_eq/V_crit ≃ 0.7. This translates into V_eq ~ 260 [(M/M_⊙)/R_o/R_⊙)]  km s^-1 (Frémat et al. 2005), which means that the equatorial velocities range from 300 to 220 km s^-1 in dwarf stars with spectral types going from B0 to A0, and from 230 to 180 km s^-1 in giant stars of the same spectral type. Searching in the catalog by Glebocki & Stawikowski (2000) for emission-less objects hotter than A2, with Vsini ≳ 200 km s^-1, apparent magnitudes V ≤ 5 mag, and declinations δ > −20°, we found 63 OB stars that can be observed with VEGA/CHARA in the framework of the present study. Since their average inclination angle is probably i ≳ 52° ~ arcsin(π/4), we may expect that many of them have true equatorial velocities in the range 250   km   s^-1 ≲ V_eq ≲ V_crit.

A detailed account of the iteration procedure used to determine R_e/R_p, η_o, α, and i from interferometric data will be given elsewhere for actual stars observed with VEGA/CHARA.

b) Parameters (α,i). Chelli & Petrov (1995a,b) and Petrov et al. (1996) highlighted the advantages of combining Fourier transforms with the differential interferometry. Domiciano de Souza et al. (2004b) showed that using the wavelength-dependent photometric barycenter, measured in angular units of interferometric intensity maps, it is possible to produce independent observational estimates of the differential rotation parameter α and the inclination angle i. The method is fairly independent of the characteristics of the observed spectral line and gives well defined indications of the sign of α. It also provides two estimates of i. Combining then data from the above class a) with either interferometric or spectroscopic information, it is possible to infer the right value of i. The method developed by Domiciano de Souza et al. (2004b) is used for slow rotators. Its implementation to gravity-darkened and α-dependent differential rapidly rotating early-type stars is currently being developed.

c) Imaging. Recent imaging of fast rotators using interferometric data has detected deviations from the standard von Zeipel gravity-darkening prescription (Monnier et al. 2007; Zhao et al. 2009). The imaging technique can therefore provide complementary information to solve a possible parameter degeneracy that may include the differential rotation in the stellar surface.

We note that while in the imaging techniques parameters such as fi_GD and the inclination angle i are inferred mostly from confidence levels in the χ² fitting surfaces, in the present attempt we try to infer α and i using direct methods. This does not preclude, however, the use of confidence levels to decide the choice of other parameters, in particular the gravity darkening exponent fi_GD.

6.3. Method to derive the parameters

We developed GIRFIT, a powerful tool to fit observed spectra of low and high dispersion (Frémat et al. 2006) with theoretical fluxes based on model atmospheres calculated with TLUSTY (Hubeny & Lanz 1995) and/or with ATLAS9 (Kurucz 1993; Castelli et al. 1997). Adopting parameterized values of fi_GD, we can thus infer the apparent quantities given in Eq. (70). Corrections for first order rotation effects can then be attempted as in Frémat et al. (2005), Zorec et al. (2005), Martayan et al. (2006), and Martayan et al. (2007) to derive the average parameters in Eq. (69) and obtain estimates of M and t. Since in most remaining calculations the mass enters as M^1/2, we can proceed to derive the other quantities with approximate values of M.

In principle, from Eq. (74) and relations of the type shown in Fig. 11b we can derive α) and i. By minimizing the differences between Eq. (77) and interferometric observations, another set of independent estimates of (η_o,α,i) can be derived.

The choice of the iterated parametric set (M,t,η_o,α,i) as a function of the gravity darkening power fi_GD can be finally decided by applying a GIRFIT-like residual minimization procedure between the observed and the FASTROT model spectra.

7. The rotational profile in the envelope

From the discussion above, it follows that if α ≠ 0, it will probably be impossible to know anything more about the characteristics of the rotation law in the convective layers beneath the stellar surface. Apart from the rigid or shellular-type rotation profiles being able to be excluded, spectroscopic and interferometric data do not provide any further information. Higher radii ratios R_e/R_p than allowed by the inferred values of α may be indicative of large values of the energy ratio τ = K/|W|, but no detailed description can be derived about the internal rotation law itself. To establish the properties of this law near the surface, at least whether they are to be likely of the speculated types generated by the condition S = S(j²), S = S(Ω²), or S = S(ϵ_Ω), no information can be obtained otherwise than from the analysis of non-radial pulsation modes excited in these regions, provided they can be unambiguously detected from time series of spectral line variations. Ando (1980) showed that the “effective depths” of p_n-modes, where n = 1−5, for stars with masses from 5   M_⊙ to 20   M_⊙ range from 0.7 to 0.9R_e, which means that the detected pulsations, mainly the splitting of modes, may carry some information on the characteristics of the rotation law in the upper envelope layers. In all cases, the interpretation of these modes rely on the knowledge of the stellar fundamental parameters, in particular the inclination angle i of the rotation axis (Floquet et al. 2000, 2002; Neiner et al. 2002; Levenhagen et al. 2003; Frémat et al. 2006) whose determination can be attempted with the methods outlined in the present paper.

8. Conclusions

Recent theoretical works suggest that rapidly rotating early-type stars should have rather deep convective layers in the envelope. This situation may favor some coupling between the convective region and the differential rotation, as happens in the Sun. We have considered the possible properties the iso-rotation curves can have in the envelope immediately under the stellar surface by assuming several conditions for the specific entropy S as a function of the: 1) specific angular momentum j, as implied by marginal Solberg-Høiland stability condition; 2) squared angular velocity Ω², which accounts for the Solar differential rotation in the convection zone; 3) specific kinetic energy ϖ²Ω², which reproduces the calculated rotational profiles in radiative envelopes. In all cases, the function that describes the gradient ∂Ω/∂θ in the stellar surface must be specified in advance. We thus assumed a simple Maunder formula with a unique differential rotation-parameter α to represent the imprint on the stellar surface produced by the differential rotation beneath the surface. In spite of the simplicity of this relation, for many theoretical reasons it would be highly interesting to infer from observations reliable orders of magnitude of α and its sign. While a value α ≠ 0 can imply that the rotation beneath the surface is neither rigid nor shellular, neither spectroscopy nor interferometry can say anything about their actual functional nature beneath the surface. Other piece of information is thus needed, perhaps non-radial pulsation data of pulsation modes excited in the intermediate and upper stellar layers.

In early-type fast rotators, we have shown that the surface temperature gradient might induce a positive gradient in the surface angular velocity, i.e., ∂Ω/∂θ > 0.

In this paper, we have summarized the kind of information that it is possible to extract from the spectroscopic and interferometric data about the differential rotation of rapidly rotating early-type stars. The equation of the surface of a star with a non-conservative rotation law was discussed. It has been shown that a differential rotation in the surface induces measurable stellar deformations. These deformations carry a surface gravitational darkening effect that needs to be studied consistently with the induced geometrical deformations. We have shown that the effects induced by the surface differential rotation can be studied consistently and in a reliable way by a combined analysis of spectroscopic and interferometric data. From spectroscopy, estimates of the differential rotation parameter α with uncertainties δα ≲ 0.15 and of the inclination angle i with δi ≲ 5° can be obtained. The interferometry can provide information about the true ratio R_e/R_p of the equatorial to the polar radii to within less than 3%, provided that a reliable estimate of the inclination angle is given from spectroscopy. The VEGA/CHARA interferometric instrument is capable of providing useful information about rotationally induced stellar deformations for nearly 60 OB rapidly rotating stars (Vsini > 200 km s^-1) whose spectra are not marred by circumstellar emission/absorption. Differential interferometry can help distinguish the α and i parameters and provide reliable independent estimates of each.

By calculating simplified models of stars at different evolutionary stages in the MS with internal rigid and differential conservative rotation laws, we have tested the use of the Roche approximation to represent the gravitational potential of rotating stars. We have concluded that we can safely use the Roche approximation for most, if not all stellar objects. Only in highly centrifugally deformed objects may this approximation not be applicable, but doubts exist as to whether these cases actually exist in Nature.

Uncertainties and constraints carried by the spectroscopic and interferometric modeling as well as the application of these methods to real stars will be presented in separate papers.

Online material

Acknowledgments

We thank the referee whose careful reading, remarks and criticisms allowed to improve the presentation of our results. We are particularly indebted to Dr. Claire Halliday for her careful language editing of this paper. Y.F. aknowledges the FNRS (Belgium) for travel assistance in the framework of contract 1.5.211.09 “Crédit aux chercheurs”. L.C. acknowledges financial support from the Agencia de Promoción Científica y Tecnológica (BID 1728 OC/AR PICT 111), from CONICET (PIP 0300) and the Programa de Incentivos G11/089 of the National University of La Plata, Argentina.

References

Appendix A: Approximate analytic representation of dV²/⟨ V² ⟩

We assume that the observed visibility of an oblate source having an equator to polar radii ratio R_e/R_p to be given for an equivalent uniform disc with angular diameter θ_o as $\mathit{V}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\left|\frac{\mathrm{2}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{x}}\right|\mathit{,}$(A.1)where J₁(x) is the Bessel function of the first kind, x = πθ_o(B/λ) where B is the baseline of the interferometer projected onto the sky and λ is the effective wavelength at which observations are carried out. To achieve an accuracy of better than 0.2%, we can use the development $\mathit{V}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\simeq }\mathrm{1}\mathrm{+}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{3}}{\mathit{a}}_{\mathit{n}}{\mathit{x}}^{\mathit{n}}$(A.2)with a₁ = 0.0173261, a₂ = −0.160464, and a₃ = 0.229252, which is valid for 0.0 ≤ x ≤ 3.8327 (first lobe of J₁(x)). Defining the argument x as a function of the average angular diameter of the source θ_o = 2 ×  [(R_e + R_p)/2] /D (where D is the distance of the source), the corresponding arguments of the visibility function in the equatorial and polar direction of the same source can be written as $\begin{array}{ccc}{\mathit{x}}_{\mathrm{p}}& \mathrm{=}& \mathit{x}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathrm{ep}}\mathrm{)}\\ {\mathit{x}}_{\mathrm{e}}& \mathrm{=}& \mathit{x}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\delta }}_{\mathrm{ep}}\mathrm{)}\\ {\mathit{\delta }}_{\mathrm{ep}}& \mathrm{=}& \frac{\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathit{/}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{-}\mathrm{1}}{\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathit{/}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{+}\mathrm{1}}\mathrm{·}\end{array}$(A.3)Since 0.0 ≤ δ_ep ≤ 1 we can use a convergent development of Eq. (A.2) in a series of powers of δ_ep. Conserving only the first order terms in δ_ep, Eq. (A.2) becomes ${\mathit{V}}_{\mathrm{\pm }\mathit{\delta }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\simeq }\mathit{V}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\pm }{\mathit{\delta }}_{\mathrm{ep}}\mathrm{\Gamma }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(A.4)where $\mathrm{\Gamma }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathrm{3}}\mathit{n}{\mathit{a}}_{\mathit{n}}{\mathit{x}}^{\mathit{n}}\mathrm{=}\mathrm{d}\mathit{V}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{/}\mathrm{d}\ln \mathit{x}$. Writing V_p = V [x × (1 + δ_ep]  and V_e = V [x × (1 − δ_ep] , we obtain $\frac{\mathrm{d}{\mathit{V}}^{\mathrm{2}}}{\mathrm{⟨}{\mathit{V}}^{\mathrm{2}}\mathrm{⟩}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{-}{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{V}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}\mathrm{)}}\mathrm{=}\mathrm{-}\mathrm{4}\begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{{\mathit{\delta }}_{\mathrm{ep}}\mathrm{\left[}\mathrm{\Gamma }\mathrm{\right(}\mathit{x}\mathrm{)}\mathit{/}\mathit{V}\mathrm{(}\mathit{x}\mathrm{\left)}\mathrm{\right]}}{\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathrm{ep}}^{\mathrm{2}}\mathrm{\left[}\mathrm{\Gamma }\mathrm{\right(}\mathit{x}\mathrm{)}\mathit{/}\mathit{V}\mathrm{(}\mathit{x}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{,}$(A.5)which is assumed to be used only when both V_e and V_p are in the respective first lobes of J₁(x), i.e. x_p ≲ 3.8327 and x_e ≲ 3.8327.

All Tables

All Figures