© ESO, 2017

1. Introduction

Rotation induces a geometrical deformation to the stars characterized by a polar flattening and an equatorial stretching. Since the radiation tends to emerge isotropically from the object, the radiation flux becomes a function of the stellar latitude known as the gravitational-darkening effect.

According to Poincaré-Wavre’s theorem, in barotropic systems any of the following statements implies the three others (Tassoul 1978): 1) the angular velocity depends only on the distance ϖ to the rotational axis, Ω = Ω(ϖ); 2) the effective gravity is derived from a total gravitation-rotational potential, Φ; 3) the effective gravity is normal to the isopycnic surfaces; 4) the isobaric and isopycnic surfaces coincide. If the stellar atmospheres of barotropic stars are also in hydrostatic and radiative equilibrium, a relation holds between the bolometric radiation flux F and the effective gravity, g_eff, known as von Zeipel’s theorem (von Zeipel 1924) $\begin{array}{ccc}\mathit{F}\mathrm{=}\mathit{c}\mathrm{\left(}\mathrm{\Phi }\mathrm{\right)}{\mathit{g}}_{\mathrm{eff}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathit{,}& & \end{array}$(1)where Φ is the total gravitation-rotational potential describing the stellar surface and c(Φ) is thus a constant. This relation holds when the radiation flux can be written in the diffusion approximation, which is an asymptotic solution for the radiation field valid at great depths in a semi-infinite atmosphere. In this case, the gravity-darkening exponent takes the value β₁ = 1.0.

It was recognized very early on that in rotating objects simultaneous hydrostatic and radiative equilibrium contradict each other (Osaki 1966). However, if radiative equilibrium is maintained forcibly in baroclinic stars, meaning that their rotation laws are non-conservative, Ω = Ω(ϖ,z) (z coordinate parallel to the rotation axis), authors have shown that β₁ = β₁(θ) ≤ 1, where θ is the co-latitude (Smith & Worley 1974; Kippenhahn 1977; Maeder 1999; Lovekin et al. 2006; Claret 2012).

For the particular case of a strict solid-body rotation, Espinosa Lara & Rieutord (2011) have shown that Eq. (1) can take another formal aspect. Using full two-dimensional stellar models (e.g. Rieutord et al. 2016), they noted that, to a high degree of approximation, the emerging radiation flux-vector F emitted by an axially symmetric stellar atmosphere in radiative equilibrium is anti-parallel to the vector of the local effective gravity g_eff. Thus, assuming this parallelism, they write $\begin{array}{ccc}{F}\mathrm{=}\mathrm{-}F\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}{{g}}_{\mathrm{eff}}\mathit{,}& & \end{array}$(2)where (r,θ) are the spherical coordinates. The function F(r,θ) is hereinafter called the gravity-darkening function (GDF). We note that before Espinosa Lara & Rieutord (2011), F(r,θ) was simply the constant c(Φ) in Eq. (1). It entered the formulation of the known paradox of von Zeipel, which stipulates the impossibility of having stable barotropic stellar models in radiative equilibrium (e.g. Roxburgh 1966, p. 211).

To determine F(r,θ) in the case of a rigidly rotating star, Espinosa Lara & Rieutord (2011) maintain the condition of radiative equilibrium written in its general form (Mihalas 1979) and solved the following differential equation $\begin{array}{ccc}\mathrm{\nabla }\mathrm{·}{F}\mathrm{=}{{g}}_{\mathrm{eff}}\mathrm{·}\mathrm{\nabla }F\mathrm{+}F\mathrm{\nabla }\mathrm{·}{{g}}_{\mathrm{eff}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(3)and obtained that β₁ depends on the co-latitude angle θ and that it is everywhere a decreasing function of the the stellar flattening, ε = 1 − R_p/R_e (R_p and R_e are polar and equatorial radii, respectively). We note that current interferometric imaging and modelling of rapidly rotating atmospheres produce inferences of the gravity-darkening exponent in intermediate-mass stars ${\mathit{\beta }}_{\mathrm{1}}^{\mathrm{obs}}\lesssim \mathrm{1}$ (Monnier et al. 2007, 2012; Zhao et al. 2009; Che et al. 2011; Domiciano de Souza et al. 2014), although within the current measurement uncertainties, these values are ${\mathit{\beta }}_{\mathrm{1}}^{\mathrm{obs}}\gtrsim \overline{{\mathit{\beta }}_{\mathrm{1}}}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}$ systematically, where $\overline{{\mathit{\beta }}_{\mathrm{1}}}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}$ is an average relation between β₁ and ε (see Eq. (40) in Sect. 7.1).

Apart from the above mentioned radiation transfer reasons, the exponent β₁ determined observationally can also be dependent on: 1) the line-of-sight angle i of stars, because β₁ is a function of θ; 2) the differential rotation of the stellar surface, which introduces a stronger dependence of effective temperature with co-latitude than in stars with solid-body rotation. In this respect, Delaa et al. (2013) noted that differential rotation introduces additional contrast on the brightness distribution over the stellar disc, which increases the dependence of β₁ with the line-of-sight angle.

We can think of the surface differential rotation as the external imprint of the internal rotation law in baroclinic stars (Espinosa Lara & Rieutord 2007, 2013). Statistical inferences concerning intermediate-mass stars suggest that their atmospheres may have differential rotation (Zorec & Royer 2012). The rotational increase of the radiative gradient in the envelope of massive rapid rotators significantly enlarges the external convective zones (Clement 1979; Maeder et al. 2008). Like in the Sun, the interaction of rotation with convection may keep driving differential rotation in these regions that in turn can be responsible for the surface differential rotation (Zorec et al. 2011, 2017).

The aim of the present work is to obtain solutions for the function F(r,θ) by imposing the condition of radiative equilibrium in the stellar atmosphere (cf. Eq. (3)), while the stellar surface rotates differentially. We explore the incidence of this rotation on the observed stellar fundamental parameters, in particular the measured Vsini parameter and the so-called gravity-darkening exponent β₁.

Fig. 1 Surface angular velocity Ω as a function of colatitude θ as given by Eqs. (4) and (5). Panel a: Υ(θ) = coskθ for α = + 0.5 and − 0.5, and several values of k from 0.5 to 10. Panel b: Υ(θ) = sinkθ for α′ = − 1 / 3 and α′ = + 1.0[ α′ = − α/ (1 + α ], and the same values of k as in (a). Colours identify the values of k.

Fig. 2 Geometrical deformation of stellar surfaces produced by differential rotation laws given by Eqs. (4) and (5). Panel a: Stellar geometries calculated using Υ(θ) = coskθ, rotational ratio η = 0.8, differential rotation parameters α = − 0.5, 0.0 (solid body rotation) and 0.5, and powers k = 0.5, 2.0 and 6.0. Panel b: Calculations done with Υ(θ) = sinkθ and the same values of η and powers k, but for α′ = − α/ (1 + α). In all cases the colours identify the values of the power k, while the line-styles characterize the values of the differential-rotation parameter α, and the corresponding α′.

2. The surface differential rotation law

Having no a priori indications on a general form of the function describing the surface differential rotation law Ω(θ), but inspired by the simplified Maunder’s relation in the solar example (e.g. Maunder & Maunder 1905, p. 819), we propose to use the following relation $\begin{array}{ccc}\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{o}}\mathrm{[}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{(}\mathit{\theta }\mathrm{\left)}\mathrm{\right]}\mathit{,}& & \end{array}$(4)where α is a free parameter called the differential rotation parameter. The case α = 0 corresponds to rigid rotation. Accordingly, we adopt two exploratory forms for the function Υ(θ), namely $\begin{array}{ccc}\mathrm{{\rm Y}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(5)where k is also a free quantity that takes any real positive value. Figure 1 shows the behaviour of the assumed angular velocity laws given by Eqs. (4) and (5) for α = + 0.5 (α′ = − 1 / 3), α = 0.0 and − 0.5 (α′ = 1.0), and several values of k from 0.5 to 10.0. The higher the value of k the larger the domain sketched as nearly a solid body rotation. This domain occurs near the equator for Υ(θ) = cos^kθ and near the pole for Υ(θ) = sin^kθ. When Υ(θ) = cos^kθ we have Ω_o = Ω_e, while for Υ(θ) = sin^kθ, Ω_o = Ω_e(1 + α)^-1. In both cases Ω_e is the angular velocity at the equator.

In order to grasp the differences inherent to the use of one of the two Υ-functions given in Eq. (5), it is necessary to compare the effects produced by rotation laws implying the same ratio Ω_p/ Ω_e. Thus, for a given value of α used with Υ(θ) = cos^kθ, we have to take α′ = − α/ (1 + α) associated with Υ(θ) = sin^kθ so that Ω(θ) = Ω_o [ 1 + α′sin^kθ ]. Obviously, when k = 2 both Υ(θ) represent the same rotation law. Nevertheless, keeping the reciprocity between α and α′, the functions Υ(θ) = cos^kθ and Υ(θ) = sin^k′θ can produce resembling behaviours of predicted parameters (although not exactly the same), because $\begin{array}{ccc}\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(6)independently of the value of α. Taking the average of the two thirds of powers k′(k,θ) calculated in the middle points of the θ-interval – neglecting the two halves of the remaining third of points situated in the extremes of the interval where in any case the predicted functions fit the right values because of the boundary conditions, it comes, for example that Υ(θ) = cos^kθ with k = 0.5 will lead to resembling results as Υ(θ) = sin^k′θ when k′ ≃ 7.0, or cos^kθ with k = 6.0 similar results as sin^k′θ where k′ ≃ 0.5 (this effect is illustrated in many figures below, as: Figs. 4a and f, or Figs. 4c and d; Figs. 5a and f, or Figs. 5c and d; Figs. 6a, d and 6b, d, or Figs. 7c, f and 7a, c; Figs. 12a and f, or Figs. 12c and d).

In principle we could let − ∞ <α< + ∞, but we shall limit the values of the differential rotation parameter to − 1 ≤ α< ∞ so as to prevent extreme cases where the pole and the equator rotate in opposite senses. Furthermore, according to statistical suggestions based on a study of Be stars (Zorec et al. 2017), in the present work we limit the values of α to the interval − 0.5 <α ≤ + 0.5. We note that the solar surface differential rotation can be sketched with α ≃ − 0.3 and Υ(θ) = cos²θ. On the other hand, Espinosa Lara & Rieutord (2013) predicted α ≳ − 0.2 for stars with masses 2 ≲ M/M_⊙ ≲ 4.

Instead of using Ω_e/ Ω_c to characterize the rotational velocity at the equator, in this paper we use the ratio η between the centrifugal and gravitational acceleration at the equator defined as follows $\begin{array}{ccc}\mathit{\eta }\mathrm{=}{\mathrm{\Omega }}_{\mathrm{e}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{e}}^{\mathrm{3}}\mathit{/}\mathit{GM}\mathrm{=}{\left({\mathrm{\Omega }}_{\mathrm{e}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{c}}\right)}^{\mathrm{2}}{\mathrm{(}{\mathit{R}}_{\mathrm{e}}\mathit{/}{\mathit{R}}_{\mathrm{c}}\mathrm{)}}^{\mathrm{3}}\mathit{,}& & \end{array}$(7)where M is the stellar mass, G the graviational constant; Ω_c and R_c are the critical angular velocity and the critical radius at the equator, respectively. To make easier any comparison of the expressions derived in the present paper with those in Espinosa Lara & Rieutord (2011) where the ratio V_e/V_K is used as the non-dimensional parameter characterizing the stellar rotation (V_e and V_K are the actual and Keplerian linear equatorial velocities, respectively), we note that η = (V_e/V_K)², so that it holds 0 ≤ η ≤ V_e/V_K≤ Ω_e/ Ω_c ≤ 1.

3. The stellar geometry and the surface gravity

When the rotation law is conservative, the geometry of the stellar surface can be described with a total gravity-rotation potential. As noted in Sect. 1, the rotation law given in Eq. (4)represents the boundary condition in the surface of an internal non-conservative rotational law Ω = Ω(r,θ). Unlike for barotropic models, in these baroclinic stars it is no longer possible to define a rotational potential. According to Maeder (2009), the surface in such objects should be the region where an arbitrary displacement ds does not imply any work done by the effective gravity g_eff, i.e. g_eff ds = 0. In the present paper we adopted this approach as it was previously done also in Zorec et al. (2011). Figure 2 shows the geometries obtained with this method when the rotation law is given by Eqs. (4) and (5), and taking into account several values of α and k. The rotation at the equator is characterized by the rotational ratio η = 0.8. In all cases of Fig. 2 the colours identify the values of the power k, while the line-styles characterize the values of the differential-rotation parameter α (and the corresponding α′). The green curves are for α = 0.

In Table 1 we give the radii ratios R_e/R_p for η = 0.2 (Ω / Ω_c = 0.69; V/V_c = 0.52) and η = 0.8 (Ω / Ω_c = 0.99; V/V_c = 0.92), for several values of k and α. In each series of parameters (η,k,α) the radii ratios were calculated also for Υ = sin^kθ and α′ = − α/ (1 + α).

Zahn et al. (2010) have shown that the level surfaces of models with shellular rotation laws have the same shape as those for solid-body rotation, including the stellar surface. Moreover, the effect of the internal mass distribution on the surface layer at equatorial critical rotation deviates from a genuine central gravitational field by no more than some 2%. The same conclusion was also put forward by Zorec et al. (2011) in models calculated using Clement’s like conservative rotation laws for kinetic-energy parameters higher than κ = K/ | W | = 0.10 (K is the kinetic rotational energy stored by the star; W is the gravitational potential energy). Thus, we can safely use the Roche approximation for the surface effective gravity g_eff and write $\begin{array}{ccc}& & {{g}}_{\mathrm{eff}}\mathrm{=}{\mathit{g}}_{\mathit{r}}\hspace{0.17em}{ê}\mathit{r}\mathrm{+}{\mathit{g}}_{\mathit{\theta }}\hspace{0.17em}{ê}\mathit{\theta }\\ & & {\mathit{g}}_{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{r}{\sin }^{\mathrm{2}}\mathit{\theta }\\ & & \end{array}$(8)where G is the gravitational constant; M is the stellar mass; Ω = Ω(θ) is the angular velocity given by Eqs. (4) and (5); ê_r and ê_θ are the unit vectors associated with the spherical coordinates (r,θ). In what follows, we use the dimensionless radial coordinate r = r/R_e, where R_e is the rotationally modified stellar equatorial radius. To simplify notation we write g_eff = | − g_eff | and introduce the dimensionless expression for the effective gravity γ(r,θ) as follows $\begin{array}{ccc}& & {\mathit{g}}_{\mathrm{eff}}\mathrm{=}{\mathrm{(}{\mathit{g}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{g}}_{\mathit{\theta }}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{⟨}\mathit{g}\mathrm{⟩}\mathit{\gamma }\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\\ & & \mathit{\gamma }\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}{\left\{\mathrm{1}\mathrm{+}\left[{\mathrm{(}\mathrm{1}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{\eta }\end{array}{\mathit{r}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right]{\sin }^{\mathrm{2}}\mathit{\theta }\right\}}^{\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \begin{array}{c}\mathrm{˜}\\ \mathit{\eta }\end{array}\mathrm{=}\mathit{\eta }{\left[\frac{\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathrm{\Omega }}_{\mathrm{e}}}\right]}^{\mathrm{2}}\\ & & \end{array}$(9)where ${\mathit{g}}_{\mathrm{o}}\mathrm{=}\mathit{GM}\mathit{/}{\mathit{R}}_{\mathrm{o}}^{\mathrm{2}}$ is the surface gravity of the stellar parent non-rotating counterpart having a spherical radius R_o and the same mass M as the rotating object.

From Eq. (9) it follows that the extreme values of γ(r,θ) are: γ_e = γ(r,π/ 2) =1 − η and γ_p = γ(r,0) =(R_e/R_p)², which are functions of (η,α) and are dependent on the function Υ(θ) adopted (see Table 1).

4. The gravity-darkening function

The gravity darkening function F(r,θ), as defined by Eq. (2), is constrained by the boundary condition obtained when r → 0, namely $\begin{array}{ccc}\underset{\mathit{r}\mathrm{\to }\mathrm{0}}{\lim }F\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathit{L}\mathit{/}\mathrm{4}\mathit{\pi GM,}& & \end{array}$(10)where L is the stellar bolometric luminosity emitted by the stellar core, whose physical properties depend on the total amount and distribution of the angular momentum inside the star (e.g. Sackmann 1970; Bodenheimer 1971; Clement 1979; Eriguchi & Mueller 1991; Uryu & Eriguchi 1994, 1995; Maeder & Meynet 2000; Deupree 2001; Jackson et al. 2005; Espinosa Lara & Rieutord 2007, 2013; Ekström et al. 2008; Maeder 2009; Fujisawa 2015). In what follows, we use the dimensionless form of F$\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}F\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\frac{\mathrm{4}\mathit{\pi GM}}{\mathit{L}}\mathrm{·}& & \end{array}$(11)From Eq. (10) we have $\begin{array}{ccc}\underset{\mathit{r}\mathrm{\to }\mathrm{0}}{\lim }\mathit{f}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathrm{1.}& & \end{array}$(12)Using spherical coordinates, the relation in Eq. (3) becomes the following linear partial differential equation: $\begin{array}{ccc}\mathit{P}\frac{\mathit{\partial f}}{\mathit{\partial r}}\mathrm{+}\mathit{Q}\frac{\mathit{\partial f}}{\mathit{\partial \theta }}\mathrm{=}\mathrm{2}\mathit{f,}& & \end{array}$(13)where the notations stand for $\begin{array}{ccc}& & \mathit{P}\mathrm{=}\mathrm{-}\frac{{\mathit{g}}_{\mathit{r}}}{\mathrm{\Omega }\left[\mathrm{\Omega }\mathrm{+}\mathrm{(}\mathit{\partial }\mathrm{\Omega }\mathit{/}\mathit{\partial \theta }\mathrm{)}\sin \mathit{\theta }\cos \mathit{\theta }\right]}\\ & & \end{array}$(14)The integration of Eq. (13) can be reduced to the integration of Lagrange’s system of ordinary linear differential equations, called also method of characteristics. Let us follow Rieutord (2016). We first eliminate the RHS of Eq. (13) by writing the equation for lnf and setting: $\begin{array}{ccc}\ln \mathit{f}\mathrm{=}\ln \mathit{G}\mathrm{+}\mathit{A}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}& & \end{array}$(15)After some calculations we find that $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}\ln \mathrm{[}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\tan }^{\mathrm{2}}{\mathit{\theta }}^{\mathrm{]}}\mathit{.}& & \end{array}$(16)Now, Eq. (13) may be rewritten $\begin{array}{ccc}\mathit{P}\frac{\mathit{\partial G}}{\mathit{\partial r}}\mathrm{+}\mathit{Q}\frac{\mathit{\partial G}}{\mathit{\partial \theta }}\mathrm{=}\mathrm{0.}& & \end{array}$(17)The characteristic lines of G in the (r,θ)-plane, are the lines where G is constant. On these lines ∂_rGdr + ∂_θGdθ = 0, so that on these lines we have $\begin{array}{ccc}\mathit{Q}\mathrm{d}\mathit{r}\mathrm{-}\mathit{P}\mathrm{d}\mathit{\theta }\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(18)which is the equation of characteristic lines. After some simplifications it reads $\begin{array}{ccc}\mathrm{-}\left({\mathit{g}}_{\mathit{\theta }}\mathit{/}\mathit{r}\right)\mathrm{d}\mathit{r}\mathrm{+}{\mathit{g}}_{\mathit{r}}\mathrm{d}\mathit{\theta }\mathrm{=}\mathrm{0.}& & \end{array}$(19)To solve (19), we need to transform it into an exact differential. We can indeed multiply the equation by any function of (r,θ). It turns out that if we multiply (19) by $\begin{array}{ccc}\mathit{M}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathit{N}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\frac{{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}}{\mathit{GM}}\mathit{,}& & \end{array}$(20)it can be integrated. Indeed, demanding that $\begin{array}{ccc}\frac{\mathit{\partial }\left(\mathrm{-}\mathit{M}{\mathit{g}}_{\mathit{\theta }}\mathit{/}\mathit{r}\right)}{\mathit{\partial \theta }}\mathrm{=}\frac{\mathit{\partial }\left(\mathit{M}{\mathit{g}}_{\mathit{r}}\right)}{\mathit{\partial r}}\mathit{,}& & \end{array}$(21)leads to the following expression for N(θ): $\begin{array}{ccc}\mathit{N}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathit{A}\frac{{\cos }^{\mathrm{2}}\mathit{\theta }}{{\mathrm{\Omega }}^{\mathrm{4}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\sin \mathit{\theta }}\mathit{,}& & \end{array}$(22)where A is an arbitrary constant.

Thus G is constant on the lines τ(r,θ) = Cst, if τ verifies: $\begin{array}{ccc}& & \frac{\mathit{\partial \tau }}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\frac{\mathit{M}{\mathit{g}}_{\mathit{\theta }}}{\mathit{r}}\mathrm{=}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{o}}^{\mathrm{2}}}{\mathit{GM}}{\cos }^{\mathrm{3}}\mathit{\theta }\\ & & \end{array}$(23)where we chose $\mathit{A}\mathrm{=}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{o}}^{\mathrm{2}}$. Since the RHS of the previous equations have been chosen such that dτ is an exact differential, the integration is straightforward and gives $\begin{array}{ccc}\mathit{\tau }\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathit{\eta }{\left(\frac{{\mathrm{\Omega }}_{\mathrm{o}}}{{\mathrm{\Omega }}_{\mathrm{e}}}\right)}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}{\cos }^{\mathrm{3}}\mathit{\theta }\mathrm{+}T\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}& & \end{array}$(24)The function T(θ), here referred to as the integral of characteristic curves, is finally given by $\begin{array}{ccc}T\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{\pi }\mathit{/}\mathrm{2}}^{\mathit{\theta }}\frac{\cos \mathit{\phi }\cot \mathit{\phi }}{{\left[\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\right]}^{\mathrm{2}}}\mathrm{d}\hspace{0.17em}\mathit{\phi }& & \end{array}$(25)that can be calculated numerically for any real positive value of k. For some low natural number k, the integral in Eq. (25) admits analytical expressions that are given in Appendix A.

Fig. 3 Function lnf(θ) given by Eq. (11) against the colatitude θ, calculated for η = 0.2 and 0.8, with both functions Υ(θ) given in Eq. (5), differential-rotation parameters α = − 0.5, 0.0, and 0.5, and for powers k = 0.5, 2.0 and 6.0. Blocks (a) and (b) are for Υ(θ) = coskθ. Blocks (c) and (d) are for Υ(θ) = sinkθ. Colours indicate the power k and they are the same for all blocks. The line-styles identify the values of α and they are the same for α′ = − α/ (1 + α) where α are the same as in (a) and (b).

Fig. 4 Panels a–c: Effective temperature Teff(θ,η,α) / ⟨ T ⟩ given by Eq. (34) as a function of the co-latitude θ when Υ(θ) = coskθ, for η = 0.8, k = 0.5, 2.0 and 6.0, and α = − 0.5, 0.0 and +0.5. Panels d–f: Same for Υ(θ) = sinkθ and the same values of k as before, but for α′ = − α/ (1 + α) where α is the same as in (a) to (c). The black-dashed lines correspond to the classic von Zeipel approximation normalized in θ = 0 at the curve corresponding to α = 0.0.

Fig. 5 Panels a–c: Effective temperature ratios ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{eq}}\mathit{/}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{p}}$ as a function of η and different values k and α when Υ(θ) = coskθ. Panels d–f: Same for Υ(θ) = sinkθ and for the same values of k as before, but for α′ = − α/ (1 + α) where α is the same as in (a) to (c). The black-dashed lines correspond to the classic von Zeipel approximation.

In this work we used only numerical estimates of T(θ). The numerical estimates of the integral in Eq. (25) were obtained proceeding to the change of variable u = lntan(φ/ 2) and using the 16-point Gauss-Legendre quadrature rule. The function T(θ) thus obtained was also calculated for all values of k for which we could obtain the analytical expressions given in Appendix A. We thus noted that the errors of the numerical estimates of Eq. (25) are δT ≲ 10^-6 at θ ~ 10^-5 rad (), and they become δT ≲ 10^-12 as soon as θ ≳ 0.07 rad (θ ≳ 8°). The precision of the numerical estimates increases even more as θ → π/ 2.

The limit expression of T(θ) for α = 0 comes directly from Eq. (25), which is also obeyed by all analytical expressions given in Appendix A: $\begin{array}{ccc}\underset{\mathit{\alpha }\mathrm{\to }\mathrm{0}}{\lim }T\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\ln \tan \frac{\mathit{\theta }}{\mathrm{2}}\mathrm{+}\cos \mathit{\theta ,}& & \end{array}$(26)which is the form of T(θ) obtained by Espinosa Lara & Rieutord (2011, see their Eq. (20)).

Because τ ≡ τ(r,θ), where r and θ are independent variables, we can use instead the pairs (r,τ) or (θ,τ) as independent variables to obtain two simplified forms of Eq. (13) valid over the τ-characteristic curves $\begin{array}{ccc}\mathit{Q}{\left(\frac{\mathit{\partial f}}{\mathit{\partial \theta }}\right)}_{\mathit{\tau }}\mathrm{=}\mathrm{2}\mathit{f}\mathrm{and}\mathit{P}{\left(\frac{\mathit{\partial f}}{\mathit{\partial r}}\right)}_{\mathit{\tau }}\mathrm{=}\mathrm{2}\mathit{f}\mathit{.}& & \end{array}$(27)The first differential equation in Eq. (27) can be easily integrated to obtain $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Psi }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}{{\tan }^{\mathrm{2}}\mathit{\theta }\mathrm{[}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{(}\mathit{\theta }\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(28)where the factor Ψ(τ) is an integration constant that carries the dependence of f with r through a transcendental function ϑ(r,θ) defined, using Eq. (24), as $\begin{array}{ccc}T\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathit{\eta }{\left(\frac{{\mathrm{\Omega }}_{\mathrm{o}}}{{\mathrm{\Omega }}_{\mathrm{e}}}\right)}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}{\cos }^{\mathrm{3}}\mathit{\theta }\mathrm{+}T\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}& & \end{array}$(29)which assumes that there is a unique relation between ϑ and τ. From Eq. (29) it appears that ϑ → θ as r → 0. Therefore, noting that τ(r,θ) is constant over a characteristic curve, the value of Ψ(τ) can be specified at any point (r,θ) on the characteristic, in particular at r = 0. Thus, making use of the boundary condition in Eq. (12) we obtain $\begin{array}{ccc}\mathrm{\Psi }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\tan }^{\mathrm{2}}\mathit{\vartheta }\hspace{0.17em}{\left[\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}\right]}^{\mathrm{2}}\mathit{.}& & \end{array}$(30)Introduced into Eq. (28), it leads to $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}{\left(\frac{\tan \mathit{\vartheta }}{\tan \mathit{\theta }}\right)}^{\mathrm{2}}{\left[\frac{\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{{\rm Y}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\right]}^{\mathrm{2}}\mathrm{·}& & \end{array}$(31)Equation (31) is the generalized version of that obtained by Espinosa Lara & Rieutord (2011), which is valid for α = 0 (their Eq. (26)).

Equations (24) and (29), together with Eq. (25), or the forms for T(θ) given in Appendix A, define ϑ(r,θ) for every chosen point (r,θ), in particular for the stellar surface represented by the function r = r(θ). To this end it is worth noting that because ϑ ≠ θ except at the extremes ϑ = θ = 0 and ϑ = θ = π/ 2, it is impossible to derive reliable values of ϑ = ϑ(r,θ) from Eq. (29) if the stellar geometry, i.e. r = r(θ), is not determined consistently with the rotation law on which depends T(θ).

We have calculated GDF numerically using Eq. (31), where r(θ) is the radius-vector that describes the stellar surface. The function ϑ(θ) is determined with Eqs. (24) and (29) by interpolation, where T(ϑ) at a given θ is taken as a linear function of lnθ. For values of ϑ when θ → π/ 2 we consider T(ϑ) as a linear function of its expression for α = 0, so that ϑ is derived by iterating Eq. (26).

The behaviour of the function f(θ) calculated for Υ = cos^kθη = 0.2, 0.8, and for several values of α and k is shown in Fig. 3. In this figure the function f(θ) calculated with Υ = sin^kθ and α′ =− α/ (1 + α) is also shown. It shows that f [ cos^kθ,α ] ≠f [ sin^kθ,α′(α) ], except for k = 2 where f(θ) is the same for both Υ(θ) functions. The function f(θ) obtained for α = 0.0 (green curves in Fig. 3) corresponds to the solution previously obtained by Espinosa Lara & Rieutord (2011).

Fig. 6 Panels a–c: Parameters Vsini determined with the FT method using the He i 4471 line produced in atmospheres having a rotation law given by Eq. (4)with Υ = coskθ and characterized by different values of k and α, against the actual Veqsini in stars having M = 9 M⊙ and t/tMS = 0.6, and rotating with several equatorial ratios η. Panels d–f: Same for the Mg ii 4481 line.

Fig. 7 Similar as Fig. 6 but for Υ = sinkθ. The cases with k = 2 are omitted as they are identical to those in Fig. 6.

From Eq. (15) we can derive the limits of f(r,θ) when θ → 0 and θ → π/ 2$\begin{array}{ccc}& & \underset{\mathit{\theta }\mathrm{\to }\mathrm{0}}{\lim }\mathit{f}\mathrm{=}{\mathit{f}}^{\mathrm{pole}}\mathrm{=}\exp \left[\frac{\mathrm{2}}{\mathrm{3}}\mathit{\eta }{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{3}}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{p}}}{{\mathrm{\Omega }}_{\mathrm{e}}}\right)}^{\mathrm{2}}\right]\\ & & \end{array}$(32)where R_p and R_e are the polar and equatorial radii, respectively. Ω_p is the polar angular velocity, which imply that Ω_p/ Ω_e = 1 + α when Υ(θ) = cos^kθ and Ω_p/ Ω_e = (1 + α′)^-1 = (1 + α) for Υ(θ) = sin^kθ. The dependence of f^pole with α and k is given through the radii ratio R_p/R_e, which is a function of them both and on the function Υ(θ) chosen. However, f^eq depends only on the ratio η. It may be worth noting that for η = 0, it is ϑ ≡ θ so that in Eq. (31) the GDF becomes f(r,θ) = 1. This also makes that from Eq. (31) we have f^pole = f^eq = 1.

For a given ratio η we obtain the same value f^eq whatever the function Υ(θ). We may wonder then whether the same value for f^pole can be derived whatever the chosen function Υ(θ) provided we use the appropriate reciprocal gravity-darkening parameter, i.e. α or α′ = − α/ (1 + α), as it is apparent form the expression of f^pole given in Eq. (32). Actually, this reciprocity is valid only for k = 2. In the next section we show that for k ≠ 2 this apparent reciprocity does not produce the same ratio R_p/R_e, and the geometry of the stellar surface is not the same either. This is due to the fact that R_e/R_p = 1 + ηΓ, with $\mathrm{\Gamma }\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\mathrm{\left[}\mathrm{\Omega }\mathrm{\right(}\mathit{\theta }\mathrm{)}\mathit{/}{\mathrm{\Omega }}_{\mathrm{e}}{\mathrm{]}}^{\mathrm{2}}\mathit{...}\mathrm{d}\mathit{\theta }$ (see Zorec et al. 2011) where, as seen from Fig. 1, Ω(θ) / Ω_e depends on Υ(θ) for k ≠ 2, even if we use α and its reciprocal α′ according to each case. Consequently, the behaviour of the resulting GDF, f(θ), over the interval 0 <θ<π/ 2 will depend in a strict sense on the chosen Υ(θ) function if k ≠ 2.

5. The effective temperature as a function of the colatitude

The magnitude of the flux-vector in Eq. (2) can be rewritten as $\begin{array}{ccc}\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi GM}}\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}& & \end{array}$(33)where σ_SB is the Štefan-Boltzmann constant. Using Eqs. (9) and (11) this relation can be translated into $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathit{T}\mathrm{⟩}{\left[\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{\gamma }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}& & \end{array}$(34)where the radial variable r = r(θ) was used as the radius-vector representing the stellar surface, so that γ becomes a function of θ only, apart from its dependence on (η,α). It is obvious that for η = 0 (no rotation), f(θ) = γ(θ) = 1 and T_eff(θ) = ⟨ T ⟩. The factor ⟨ T ⟩ is the uniform effective temperature over a sphere of radius R_e(η) radiating the luminosity L$\begin{array}{ccc}\mathrm{⟨}\mathit{T}{\mathrm{⟩}}^{\mathrm{4}}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{SB}}}\mathrm{·}& & \end{array}$(35)It must be understood that when interpreting observations the quantity ⟨ T ⟩ is a free parameter, because the bolometric luminosity L depends on physical conditions that characterize the core of rotating stars. This luminosity depends, among other factors, on the total angular momentum stored in the star and on its internal distribution, which are unknown. In Figs. 4a to c T_eff(θ,η,α) / ⟨ T ⟩ is shown as a function of the co-latitude θ when Υ(θ) = cos^kθ, η = 0.8, k = 0.5, 2.0 and 6.0, and for several values of α ranging from − 0.5 to + 0.5. Figures 4d to f show T_eff(θ,η,α) / ⟨ T ⟩ against θ for Υ(θ) = sin^kθ and (k,α′), where α′ = − α/ (1 + α). The green curves correspond to T_eff(θ) / ⟨ T ⟩ for α = α′ = 0, previously obtained for rigid rotators by Espinosa Lara & Rieutord (2011). In Fig. 5 we plot the effective temperature ratios ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{eq}}\mathit{/}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{p}}$ as a function of η and different values k and α, when Υ(θ) = cos^kθ or Υ(θ) = sin^kθ.

Fig. 8 Effective temperature given as ln [ Teff(θ) / ⟨ T ⟩ ] 4 against ln [ geff(θ) / ⟨ g ⟩ ] for η = 0.8. Panels a–c: are for Υ(θ) = coskθ with k = 0.5, 2.0 and 6.0, respectively and α ranging from − 0.5 to + 0.5. Panels d–f: are for Υ(θ) = sinkθ with k = 0.5, 2.0 and 6.0, respectively and α′(α) ranging from − 1.0 to + 1.0. In all cases the left-bottom corner corresponds to the equator and the right-upper corner to the pole.

6. The V sini parameter

Vsini is a parameter that mainly reflects the global broadening of a spectral line. It was shown in Zorec et al. (2017) that this broadening is produced by monochromatic contributions from curved loci of points over the observed stellar hemisphere having the same Doppler displacement. The shape of these curves depend on the differential rotation law and on the inclination angle, while for rigid rotation these curves are straight lines which are independent of the inclination angle. Due to the differentiated sensitivity of spectral lines to the physical formation conditions (Thomas 1983; Zorec et al. 2017), according to the spectral line, the different regions on the observed stellar hemisphere do not contribute with the same efficiency to the rotational broadening because of the non uniform distribution of temperature and gravity. To illustrate this effect, we have calculated He i 4471 and Mg ii 4481 lines broadened by differential rotators with rotation laws given by Eq. (4) and both forms of Υ(θ) given in Eq. (5). For this exercise we used a model star with mass M = 9 M_⊙ and fractional age t/t_MS = 0.6 (t_MS is the time a rotating star spends in the main sequence evolutionary phase). Rotation laws are characterized by several parameters k, α, and different equatorial velocities parametrized according to η. Applying the classic Fourier transform method we determined the Vsini parameters shown in Figs. 6 and 7 where we compare the resulting V_eqsini with the actual input V_eqsini corresponding to the given values of η and inclination angles. We note that deviations from the identity relation (y = x line) also exist for rigid rotators. They are due to the gravitational darkening effect, which is not taken into account in the classic rigid-rotation broadening function G_R(λ) (Gray 1992).

Because it has become of common use today, the classic method based on the FT was employed to determine the Vsini. This method assumes, however, that the observed rotationally broadened line profile F(λ) can be represented as the convolution with an analytical rotation broadening-function G_R(λ) of a flux line-profile F(λ) emitted by a non-rotating star having uniform effective temperature and gravity. The specific intensity I(λ,μ) contributing to the observed line flux F(λ) is thus the same over the entire stellar disc, which obviously does not occur in rotating stars. The derivation of an analytical expression for G_R(λ) that is independent of the inclination angle requires that the star be spherical and behaves as a rigid rotator (Gray 1975, 1992), i.e. the monochromatic Doppler displacements are produced over straight strips of constant radial velocities. The Vsini is then obtained by comparing the zeros of the flux line-profile FT F(λ) with the corresponding zeros of the FT of the function G_R(λ).

In Zorec et al. (2017) it has been shown that depending on the spectral line and on the value of α, the FT of lines may have unusual shapes, making the zeros difficult to identify and to interpret. This difficultly is a direct consequence of the inconsistency raised by the use an analytic expression G_R(λ) adapted for uniform rigid rotators, which does not take into account: a) the actual shape of curves of constant radial velocity contributing to the individual rotational Doppler-displacements, and b) the non uniformity of the effective temperature and surface effective gravity determining the local monochromatic specific intensities I(λ,μ).

Figures 6 and 7 reveal that the Vsini obtained using the FT-method can be overestimated when α> 0, and underestimated for α< 0, as compared to its value for rigid rotators. We notice that for Υ(θ) = cos^kθ and a given | α | the deviations are larger when α< 0 than for α> 0. For Υ(θ) = sin^kθ the deviations reverse according to the α-corresponding reciprocal parameter α′. Finally, according to the spectral line used, the obtained Vsini is not exactly the same. To apprehend what can be the possible deviations that the current Vsini determinations may incur, let us recall that in the Sun α ≃ − 0.3 and Υ(θ) = cos²θ, and that predictions made by Espinosa Lara & Rieutord (2013) for rapid rotators of masses 2 ≲ M ≲ 4 M_⊙ foresee α ≳ − 0.2. Facing the likely possibility that stars are differential rotators, these results warn us that the current interpretations based on the rigid rotation hypothesis could be misleading. However, detailed modelling of spectral lines can be made today to account for the noted differences in the Vsini obtained from different spectral lines. They could help us to derive indications on the value of the α parameter.

Needless to say that the classic method of determining Vsini based on the FWHM (full-width at half-maximum) produces the same type of deviations, because the comparisons are made also with the rotational broadening of spectral lines produced by classic spherical rigid rotators.

Fig. 9 Panels a–c: Gravity-darkening slope ${\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$ against θ defined in Eq. (38) for Υ(θ) = coskθ with k = 0.5, 2.0, 6.0, and α = − 0.5, − 0.25, 0.0 (green lines), +0.25 and +0.5; d)–f) same, but for Υ(θ) = sinkθ and α′ = + 1.0, + 1 / 3, 0.0 (green lines), − 1 / 5 and − 1 / 3; g)–i) as from (a) to (c), but for the term ${\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$; j)–l) as from (d) to (f), but for the term ${\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$.

7. The gravity-darkening exponent

The gravity-darkening exponent was introduced with the purpose of representing in a simple way the distribution of the local effective temperature T_eff(θ) over a rotating stellar surface. As noted in Sect. 1, the known von Zeipel’s relation that holds for barotropic objects $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{F}}{\mathit{g}}_{\mathrm{eff}}^{{\mathit{\beta }}_{\mathrm{1}}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}& & \end{array}$(36)is characterized by the constant gravity-darkening exponent β₁, and the constant C_F which depends on the gravity-rotation equipotential of the stellar surface. The notation β₁ is used here to distinguish it from the exponent β = β₁/ 4 that is also frequently employed to characterize the gravity-darkening effect. In the following section, we will comment on some drawbacks of using Eq. (36).

Fig. 10 Average gravitational-darkening slope ⟨ B1 ⟩ given by Eq. (40) as function of η, calculated respectively for: a) Υ(θ) = coskθ; b) sinkθ. The colours identify the differential rotation parameter α, while k is identified using different line-styles.

7.1. Dependence of the gravity-darkening exponent with colatitude

Up to now, the relation in Eq. (36) has been used to interpret the observations of rapidly rotating stars regardless of whether they are barotropic or baroclinic. However, β₁ and C_F cannot be assumed independent of the colatitude-angle θ even in barotropic model-stars having atmospheres in radiative equilibrium (Espinosa Lara & Rieutord 2011; Rieutord 2016).

To inquire on the effects that a surface differential rotation can induce on β₁ and C_F, Fig. 8 shows the relation between lnT_eff(θ) calculated with Eq. (34) against lng_eff(θ) using both forms of Υ(θ) in Eq. (5). The model-stars are assumed to rotate with η = 0.8 with several values of α and k. Obviously, similar behaviours are noted for other values of η.

Because the curves in Fig. 8 are not straight lines, neither C_F nor β₁ can be considered as genuine constants, even though a quasi-linear relation between lnT_eff(θ) and lng_eff(θ) seems to exist over some intervals of θ. To emphasize more clearly the lack of constancy of β₁ and C_F, we define a pseudo-exponent or gravity-darkening slope b₁(θ) as $\begin{array}{ccc}\ln {\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\ln {\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}\ln \left[{\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{/}{\mathit{\sigma }}_{\mathrm{SB}}\right]\mathit{,}& & \end{array}$(37)and consider that between two successive, close discrete points θ_i and θ_{i + 1}, ${\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$ and ${\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$ can be assumed constant; i.e.: ${\mathit{b}}_{\mathrm{1}}\mathrm{(}{\overline{{\mathit{\theta }}_{\mathit{i}}}}^{\mathrm{)}}\mathrm{=}\mathrm{4}\frac{\ln \left[{\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathit{/}{\mathit{T}}_{\mathrm{eff}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\right]}{\ln \left[{\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathit{/}{\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathrm{\right)}\right]}\mathrm{;}\overline{{\mathit{\theta }}_{\mathit{i}}}\mathrm{=}\mathrm{(}{\mathit{\theta }}_{\mathit{i}}\mathrm{+}{\mathit{\theta }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{2.}$(38)Inserting the values ${\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\overline{{\mathit{\theta }}_{\mathit{i}}}\mathrm{\right)}$ in Eq. (37), the corresponding $\ln {\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{{\mathit{\theta }}_{\mathit{i}}}\mathrm{\right)}$ can be estimated. The inferred quantities ${\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\overline{{\mathit{\theta }}_{\mathit{i}}}\mathrm{\right)}$ and $\ln {\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{{\mathit{\theta }}_{\mathit{i}}}\mathrm{\right)}$ are shown in Fig. 9. They clearly demonstrate that both b₁ and C_F are dependent on θ and that consequently they cannot be used as constants in Eq. (36).

7.2. The observed gravity-darkening slope

7.2.1. Neglecting the viewing angle

Replacing in Eq. (37) the expressions for g_eff and T_eff given respectively by Eqs. (9) and (34), and imposing that b₁ and C_F do not depend on θ, it follows that $\begin{array}{ccc}{\mathit{B}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\ln \mathrm{[}\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{/}{\mathit{f}}^{\mathrm{pole}}\mathrm{]}}{\ln \left[\mathit{\gamma }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{/}{\mathit{\gamma }}^{\mathrm{pole}}\right]}\mathit{,}& & \end{array}$(39)where f^pole and γ^pole are the values of the functions f(θ) and γ(θ) in θ = 0. The function B₁(θ) behaves in a similar way as b₁(θ), so that it is needless to reproduce it graphically. Instead, we can replace f(θ) and γ(θ) in Eq. (39) by their values at θ = π/ 2 (see Eq. (32)) to obtain a kind of average ⟨ B₁(θ) ⟩ that takes the form $\begin{array}{ccc}\mathrm{⟨}{\mathit{B}}_{\mathrm{1}}\mathrm{⟩}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\left[\frac{\ln \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}\mathrm{+}\mathit{\eta }{\mathrm{(}{\mathit{R}}_{\mathrm{p}}\mathit{/}{\mathit{R}}_{\mathrm{e}}\mathrm{)}}^{\mathrm{3}}{\mathrm{(}{\mathrm{\Omega }}_{\mathrm{p}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{e}}\mathrm{)}}^{\mathrm{2}}}{\ln [\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}\mathrm{+}\mathrm{2}\ln \mathrm{(}{\mathit{R}}_{\mathrm{p}}\mathit{/}{\mathit{R}}_{\mathrm{e}}\mathrm{)}}\right]\mathit{,}& & \end{array}$(40)which is a generalized version of ⟨ B₁ ⟩ obtained by Espinosa Lara & Rieutord (2011) for rigid rotators. In Fig. 10, ⟨ B₁ ⟩ is shown as a function of η for − 0.5 ≤ α ≤ 0.5, k = 0.5, 2, 6, and for both functions Υ(θ). The expression of ⟨ B₁ ⟩ for α = 0 was sometimes used in the literature to compare the exponents b₁ derived from interferometric observations.

Fig. 11 Relation between the apparent 4ln [ Teff(i) / ⟨ T ⟩ ] and ln [ geff(i) / ⟩ ⟨ g ⟩ ]; a)–c) for Υ(θ) = coskθ with α from − 0.5 to +0.5; d)–f) for sinkθ and α′ =− α/ (1 = α). All curves are for η = 0.8, with k = 0.5, 2.0, 6.0. The colour-codes according to α and α′ indicated in the respective blocks for k = 2 hold for all curves having the same Υ(θ). The equator-on seeing directions (i → π/ 2) are in the lower left corner of each block, while the polar-on directions (i → 0) are in the upper right corner.

7.2.2. Considering the viewing angle

Any parameter describing a physical property of a rotating stellar atmosphere that depends on the co-latitude θ, automatically becomes a function of the viewing angle i when it is determined from observational data of angularly unresolved stars. This must also happen to the gravity-darkening slope given by Eq. (39). Such a quantity should thus depend on the apparent average effective temperature and gravity that characterize the observed stellar hemisphere, which are named hereafter as ⟨ T_eff(i) ⟩ and ⟨ g_eff(i) ⟩, respectively.

Let us define ⟨ T_eff(i) ⟩ and ⟨ g_eff(i) ⟩ so as to obtain a slope that might be assimilated with an empirically derived gravity-darkening slope. The apparent effective temperature ⟨ T_eff(i) ⟩ can be written as $\begin{array}{ccc}\mathrm{⟨}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{(}\mathit{i}{\mathrm{)}}^{\mathrm{⟩}}\mathrm{=}\frac{\mathit{L}\mathrm{\left(}\mathit{i}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathrm{SB}}\mathrm{𝒮}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{,}& & \end{array}$(41)where is the effective emitting area of the i-dependent rotationally deformed stellar hemisphere. In this expression, the apparent bolometric luminosity L(i) emitted by the stellar disc projected towards the observer is given by (Collins 1965; Maeder & Peytremann 1970; Collins 1973) $\begin{array}{ccc}\mathit{L}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{e}}^{\mathrm{2}}{\mathrm{􏽉}}_{\mathit{S}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\left(\frac{{n̂}\mathrm{·}{\iota ̂}}{{n̂}\mathrm{·}{ê}\mathit{r}}\right){\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\sin \mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi ,}& & \end{array}$(42)where is the unit vector normal to the stellar surface; is the unit vector directed towards the observer; ê_r is the unit vector associated with the r-spherical coordinate; S(i) indicates the observed region of the stellar hemisphere over which the integration is carried out. μ is the directional cosine ; I(μ) is the bolometric specific intensity of the radiation for which we use a quadratic form $\begin{array}{ccc}\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{1}}\mathrm{[}\mathrm{1}\mathrm{-}{\mathit{ϵ}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }\mathrm{)}\mathrm{-}{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{)}}^{\mathrm{]}}\mathit{,}& & \end{array}$(43)whose limb-darkening coefficients ϵ_1,2 have been interpolated in the tables given by Claret (2000) for each pair of local fundamental parameters [T_eff(θ),g_eff(θ)]. Knowing that for each point over the stellar surface we have $\mathrm{\left|}{F}\mathrm{\right|}\mathrm{=}\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{1}}\mathit{I}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathit{\mu }\hspace{0.17em}\mathrm{d}\mathit{\mu }$, the intensity I₁ = I(μ = 1) can be written as $\begin{array}{ccc}{\mathit{I}}_{\mathrm{1}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}}{\mathit{\pi }}\left[\frac{\mathit{F}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{6}\mathrm{-}\mathrm{2}{\mathit{ϵ}}_{\mathrm{1}}\mathrm{-}\mathrm{3}{\mathit{ϵ}}_{\mathrm{2}}}\right]\mathit{,}& & \end{array}$(44)which depends on f(θ) and γ(θ) through F(θ) given in Eq. (33).

Consistently with the definition of ⟨ T_eff(i) ⟩, we can use Eq. (33) to estimate ⟨ g_eff(i) ⟩ using the expression $\begin{array}{ccc}\mathrm{⟨}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{(}\mathit{i}{\mathrm{)}}^{\mathrm{⟩}}\mathrm{=}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi GM}{\mathit{\sigma }}_{\mathrm{SB}}}\mathrm{⟨}\mathit{f}{\mathit{g}}_{\mathrm{eff}}{\mathrm{⟩}}_{\mathit{i}}\mathit{,}& & \end{array}$(45)that with Eqs. (41) and (42) leads to $\begin{array}{ccc}\mathrm{⟨}{\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\frac{\mathrm{⟨}\mathit{f}\hspace{0.17em}{\mathit{g}}_{\mathrm{eff}}{\mathrm{⟩}}_{\mathit{i}}}{\mathrm{⟨}\mathit{f}{\mathrm{⟩}}_{\mathit{i}}}\mathit{,}& & \end{array}$(46)where ⟨ f ⟩ _i is given by $\begin{array}{ccc}\mathrm{⟨}\mathit{f}{\mathrm{⟩}}_{\mathit{i}}& \mathrm{=}& {}^{\mathrm{􏽉}}\mathit{S}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{\phi }\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\left(\frac{{n̂}\mathrm{·}{\iota ̂}}{{n̂}\mathrm{·}{ê}\mathit{r}}\right){\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\sin \mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\mathrm{d}{\mathit{\phi }}^{\mathrm{/}}\mathrm{𝒮}\mathrm{(}\mathit{i}{\mathrm{)}}^{}\\ \mathit{\phi }\mathrm{\left(}\mathit{\mu }\mathrm{\right)}& \mathrm{=}& \end{array}$(47)For the numerical estimates of the apparent ⟨ T_eff(i) ⟩ and ⟨ g_eff(i) ⟩ the knowledge of absolute values of T_eff(θ) and g_eff(θ) are required. To get a rough insight on their amplitude variation, we have assumed model-stars with mass M = 9 M_⊙ and fractional age t/t_MS = 0.6 (t_MS is the time spent by a rotating star in the main sequence phase). The bolometric luminosity L was calculated with the rotationally induced mass-compensation effect parametrized as in Frémat et al. (2005), where it is assumed that the stellar core rotates rigidly. The equatorial radii ratio R_e/R_o (R_o is the equatorial radius of a parent non-rotating object having the same mass) is given by the models of Zorec et al. (2011).

The relations between ln [ ⟨ T_eff(i) ⟩ / ⟨ T ⟩ ] and ln [ ⟨ g_eff(i) ⟩ / ⟨ g ⟩ ] shown in Fig. 11 are calculated using Υ(θ) = cos^kθ and sin^kθ, for η = 0.8, k = 0.5, 2.0, 6.0, and − 0.5 ≤ α ≤ − 0.5. While Fig. 8 represents the variation of T_eff against g_eff over the stellar surface, Fig. 11 represents these parameters integrated over a stellar surface seen according to different viewing angles. The upper right corner corresponds to pole-on seeing directions and the lower left corners to the equator-on directions. For a given geometrical deformation of stars parameterized here with the ratio η, it turns out that the resolution among these curves according to the differential parameter α is a strong function of the power k used in Υ(θ).

From Fig. 11 it seems legitimate to assume relations of the type $\begin{array}{ccc}\ln \mathrm{⟨}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathit{♭}\mathrm{\left(}\mathit{i}\mathrm{\right)}\ln \mathrm{⟨}{\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\mathit{A}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{⟩}& & \end{array}$(48)for each [ Υ,k,α ] so as to represent the dependence of the slope ♭(i) with the stellar inclination, as we did for b₁(θ) with θ, i.e. by taking here ♭(i) and ⟨ A(i) ⟩ as constants between two successive inclinations i_j and i_{j + 1} in a given discretization of this variable in 0° ≤ i ≤ 90°. The values of ♭(i) obtained for three inclination angles: i = 0°, 45° and 90° against the ratio η are shown in Fig. 12. Figures 12a–c are for Υ = cos^kθ with k = 0.5, 2 and 6, while Figs. 12d–f are for Υ = sin^kθ and the same exponents k. In these figures we note that the ♭(i) curves have a global behaviour that is similar to ⟨ B₁ ⟩ in Fig. 10, except for their spread as a function of i.

In Fig. 13 we show the values of ${\mathit{b}}_{\mathrm{1}}^{\mathrm{obs}}$ determined with interferometric data (see Domiciano de Souza et al. 2014), where the interpretation of observations was done using Eq. (36). These points have been originally given as a function of the stellar flattening ε = 1 − R_p/R_e. We have converted ε into η considering, as the authors themselves did, that stars are rigid rotators, i.e. R_e/R_p = 1 + η/ 2, which imply that η = 2ε/ (1 − ε). In this diagram the curves ♭₁(i) are reproduced for Υ(θ) = cos^kθ with k = 2 and different values of α, in particular those suited for rigid rotators, i.e. α = 0 (green curves). In Sect. 8 we discuss the reason of using k = 2. Taking into account the observational uncertainties, we can conclude that half of the observed points are located outside the zone limited by ♭₁(i = 0) and ♭₁(i = π/ 2) of rigid rotators. Two of them are likely situated in the region of α> 0, while there is one in the zone of α< 0, which indicates that other solutions than rigid rotation might also suit these objects.

Fig. 12 Gravitational-darkening slope ♭1(i) against the rotational ratio η. The colours identify the values of the differential parameter α, and the line-styles are for the inclination angles. In a) to c) the slope ♭1(i) was calculated with Υ(θ) = coskθ and from d) to f) b1(i) was obtained with Υ(θ) = sinkθ. As throughout in this work, rigid rotation (α = 0) is identified with green lines.

Fig. 13 Comparison of the empirical gravitational-darkening slope ${\mathit{♭}}_{\mathrm{1}}^{\mathrm{obs}}$ (see reference in Domiciano de Souza et al. 2014) with the curves representing the ♭1(i) of rotators with k = 2, Υ(θ) = cosθ, α = − 0.5, 0.0 (green curves) and α = + 0.5, and inclination angles i = 0°, 45° and 90° (dentified with different line styles). The dashed black dashed line is for ♭1(i) = 1.

However, according to the diagrams in Fig. 12 it clearly appears that there is not a unique way to characterize the studied stars with β₁ as rotators, since many solutions can be envisioned in terms of [ Υ,k,α ]. On the other hand, there is a strong concern with the interpretation of the observed values ${\mathit{♭}}_{\mathrm{1}}^{\mathrm{obs}}$, because β₁ and C_F were considered as constants in the reduction of data. Due to these ambiguities, it would perhaps be wiser to abandon the parameter β₁ to model fluxes and rather concentrate on determining them using the differential rotation parameter α and the inclination angle i (Domiciano de Souza et al. 2004).

8. Discussion and conclusion

Continuing the work initiated by Espinosa Lara & Rieutord (2011) to represent the distribution of the effective temperature in rotating stars, where the only requirement is that the outermost layers of the envelope be in radiative equilibrium and the radiation flux vector be dominantly anti-parallel to the effective surface gravity, we have discussed the case where these layers are in differential rotation. We have studied the particular case where the angular velocity is given by Eqs. (4) and (5), and obtained relationships that generalize those previously derived by Espinosa Lara & Rieutord (2011) for rigid rotators. As a further step, we have studied the effects of the surface differential rotation on the determination of the Vsini parameter and concluded that strong deviations can be expected. They seem to be large enough to be exploited empirically. We also find that a constant gravity-darkening exponent is inappropriate for describing stars as rotators. In fact, as this exponent is a function of the stellar co-latitude (Rieutord 2016), it becomes a function of the viewing angle that induces ambiguities in the interpretation of observations. The modelling of observed parameters of rotating stars must then be done directly through the gravity-darkening-function (GDF), f(θ), which describes the distribution of the effective temperature as indicated by Eqs. (34), (31), (29) and (25), for which the stellar geometry must be calculated consistently. According to discussions carried out in this work, the interpretation of observations would then require one to determine the set of quantities (η,α,Υ) instead of the simple pair (η,β₁).

In the approach where differential rotation is considered, a significant uncertainty remains in the choice of the function Υ(θ). Indeed, up to now, and except for the Sun, no empirical indication exists on the actual nature of Υ(θ) in stars. The quantity that needs to be inferred first is a parameter related to the average surface gradient ⟨ d Ω / d θ ⟩. To this end, the simplest expression for Υ(θ) could be 2θ/π. This makes α the unique free parameter characterizing the surface differential rotation, but in this case the calculation of the integral T(θ) (Eq. (25)) has some annoying numerical subtleties. However, we have shown in Sect. 2 that for the same ratio Ω_p/ Ω_e, no significant differences in the prediction or analysis of observed parameters may arise if we use only Υ(θ) = cos^kθ or sin^k′θ, because reliable correspondences can be established between (α,k) and (α′,k′). According to this finding, we can adopt Υ(θ) = cos^kθ and show easily that the weighted average gradient over the stellar surface $\begin{array}{ccc}⟨\mathrm{d}\left(\mathrm{\Omega }\mathit{/}{\mathrm{\Omega }}_{\mathrm{o}}\right)\mathit{/}\mathrm{d}\mathit{\theta }⟩\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\left[\mathrm{d}\left(\mathrm{\Omega }\mathit{/}{\mathrm{\Omega }}_{\mathrm{o}}\right)\mathit{/}\mathrm{d}\mathit{\theta }\right]\sin \mathit{\theta }\mathrm{d}\hspace{0.17em}\mathit{\theta }& & \end{array}$(49)is equal to ⟨ d(Ω / Ω_o) / dθ ⟩ = 2α/π if k = 2.55. Although nothing prevents us from making numerical models with k = 2.55, the results for T_eff(θ), Vsini and other parameters, show that no significant differences occur on the predicted or analyzed parameters if we simply adopt k = 2. Moreover, the representation Υ(θ) = cos²θ of Maunder & Maunder (1905) accounts for the solar surface differential rotation. Two first order indications could then be reliably drawn from observations: a) whether the surface differential rotation in stars is significant or not; b) what is its character: accelerated towards the pole or accelerated towards the equator. We note that for some forms of Υ(θ,k), in particular for k = 2, analytical expressions of T(θ) exist and can help to proceed in a similar way as for rigid rotators with the formulas developed by Espinosa Lara & Rieutord (2011).

Acknowledgments

We are thankful to the referee for his/her attentive reading of the manuscript and for his/her suggestions to improve the presentation of our results.

References

Appendix A: The functions T(θ)

The analytical expressions obtained of the function T(θ) with Υ(θ) = cos^k(θ) by integrating Eq. (23) for the natural numbers k = 1, 2, 3 and 4 are:

k = 1

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(A.1)

k = 2

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(A.2)

k = 3

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(A.3)

k = 4

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$The corresponding analytical expressions of T(θ) with Υ(θ) = sin^k(θ) by integrating Eq. (23) for the natural numbers k = 1 and 2, are:

k = 1

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(A.5)

k = 2

$\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(A.6)In spite of the formal difference between Eqs. (A.2) and (A.6), they are equivalent. They can be converted to each other knowing that sin²θ + cos^kθ = 1, making the transformation α into − α/ (1 + α) and considering the sign of α.

Appendix B: Limit forms of f (r,θ)

Integrating Eq. (25) by parts, we obtain the following expression equivalent to Eq. (24) $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$(B.1)As θ → 0 and θ → π/ 2, also ϑ → 0 and ϑ → π/ 2, respectively. Moreover, θ → 0 both θ and ϑ are infinitesimals of the same order, so that we can write ϑ = pθ where p is constant. Expanding the function tanx for small values of x and knowing that αΥ < 1, we derive $\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(B.2)Taking the limit of this expression for θ → 0, and identifying r(θ = 0 with stellar surface it comes that $\begin{array}{ccc}\begin{array}{c}\\ \end{array}& & \end{array}$(B.3)which introduced into Eq. (31) leads to $\begin{array}{ccc}\begin{array}{c}\\ \end{array}& & \end{array}$(B.4)where Ω_p is the polar angulat velocity, which carries that Ω_p/ Ω_e = 1 + α when Υ(θ) = cos^k(θ) and Ω_p/ Ω_e = (1 + α)^-1 for Υ(θ) = sin^k(θ). We note that the dependence of f^pole with α and k is given through the radii ratio R_p/R_e.

To calculate the limit of f(r,θ) as θ → π/ 2 we make the change of variables ϑ → π/ 2 − ϵ_ϑ and θ → π/ 2 − ϵ_θ, where ϵ_ϑ → 0 and ϵ_θ → 0 are infinitesimals of the same order. This operation implies the following changes of functions $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \\ & & \\ & & \\ & & \end{array}& & \end{array}$(B.5)and similar expressions for the functions dependent on ϑ. Introducing this terms into Eq. (B.1), it is obtained $\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(B.6)which for ϵ_θ → 0 reduces to $\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(B.7)where $\begin{array}{ccc}\begin{array}{c}\\ & & \end{array}& & \end{array}$(B.8)and Ω_o/ Ω_e = 1 when Υ(θ) = cos^k(θ), or Ω_o/ Ω_e = 1 / (1 + α) if Υ(θ) = sin^k(θ).

Finally, with Eqs. (B.7) and (31) we obtain $\begin{array}{ccc}\begin{array}{c}\\ \end{array}& & \end{array}$(B.9)which is independent of α and k, and where for the stellar surface it is lim_{θ → π/ 2}r = 1.

The forms Eqs. (B.4) and (B.9) can also be derived using the analytical expressions from Eqs. (A.1) to (A.6) by using series expansions for lntanx, arctanx and ln(1 ± x) valid for x< 1 and relations in Eq. (B.5).

All Tables

All Figures

	Fig. 1 Surface angular velocity Ω as a function of colatitude θ as given by Eqs. (4) and (5). Panel a: Υ(θ) = coskθ for α = + 0.5 and − 0.5, and several values of k from 0.5 to 10. Panel b: Υ(θ) = sinkθ for α′ = − 1 / 3 and α′ = + 1.0[ α′ = − α/ (1 + α ], and the same values of k as in (a). Colours identify the values of k.
In the text

Fig. 2 Geometrical deformation of stellar surfaces produced by differential rotation laws given by Eqs. (4) and (5). Panel a: Stellar geometries calculated using Υ(θ) = coskθ, rotational ratio η = 0.8, differential rotation parameters α = − 0.5, 0.0 (solid body rotation) and 0.5, and powers k = 0.5, 2.0 and 6.0. Panel b: Calculations done with Υ(θ) = sinkθ and the same values of η and powers k, but for α′ = − α/ (1 + α). In all cases the colours identify the values of the power k, while the line-styles characterize the values of the differential-rotation parameter α, and the corresponding α′.

In the text

Fig. 3 Function lnf(θ) given by Eq. (11) against the colatitude θ, calculated for η = 0.2 and 0.8, with both functions Υ(θ) given in Eq. (5), differential-rotation parameters α = − 0.5, 0.0, and 0.5, and for powers k = 0.5, 2.0 and 6.0. Blocks (a) and (b) are for Υ(θ) = coskθ. Blocks (c) and (d) are for Υ(θ) = sinkθ. Colours indicate the power k and they are the same for all blocks. The line-styles identify the values of α and they are the same for α′ = − α/ (1 + α) where α are the same as in (a) and (b).

In the text

Fig. 4 Panels a–c: Effective temperature Teff(θ,η,α) / ⟨ T ⟩ given by Eq. (34) as a function of the co-latitude θ when Υ(θ) = coskθ, for η = 0.8, k = 0.5, 2.0 and 6.0, and α = − 0.5, 0.0 and +0.5. Panels d–f: Same for Υ(θ) = sinkθ and the same values of k as before, but for α′ = − α/ (1 + α) where α is the same as in (a) to (c). The black-dashed lines correspond to the classic von Zeipel approximation normalized in θ = 0 at the curve corresponding to α = 0.0.

In the text

Fig. 5 Panels a–c: Effective temperature ratios ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{eq}}\mathit{/}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{p}}$ as a function of η and different values k and α when Υ(θ) = coskθ. Panels d–f: Same for Υ(θ) = sinkθ and for the same values of k as before, but for α′ = − α/ (1 + α) where α is the same as in (a) to (c). The black-dashed lines correspond to the classic von Zeipel approximation.

In the text

	Fig. 6 Panels a–c: Parameters Vsini determined with the FT method using the He i 4471 line produced in atmospheres having a rotation law given by Eq. (4)with Υ = coskθ and characterized by different values of k and α, against the actual Veqsini in stars having M = 9 M⊙ and t/tMS = 0.6, and rotating with several equatorial ratios η. Panels d–f: Same for the Mg ii 4481 line.
In the text

	Fig. 7 Similar as Fig. 6 but for Υ = sinkθ. The cases with k = 2 are omitted as they are identical to those in Fig. 6.
In the text

Fig. 8 Effective temperature given as ln [ Teff(θ) / ⟨ T ⟩ ] 4 against ln [ geff(θ) / ⟨ g ⟩ ] for η = 0.8. Panels a–c: are for Υ(θ) = coskθ with k = 0.5, 2.0 and 6.0, respectively and α ranging from − 0.5 to + 0.5. Panels d–f: are for Υ(θ) = sinkθ with k = 0.5, 2.0 and 6.0, respectively and α′(α) ranging from − 1.0 to + 1.0. In all cases the left-bottom corner corresponds to the equator and the right-upper corner to the pole.

In the text

Fig. 9 Panels a–c: Gravity-darkening slope ${\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$ against θ defined in Eq. (38) for Υ(θ) = coskθ with k = 0.5, 2.0, 6.0, and α = − 0.5, − 0.25, 0.0 (green lines), +0.25 and +0.5; d)–f) same, but for Υ(θ) = sinkθ and α′ = + 1.0, + 1 / 3, 0.0 (green lines), − 1 / 5 and − 1 / 3; g)–i) as from (a) to (c), but for the term ${\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$; j)–l) as from (d) to (f), but for the term ${\mathit{C}}_{\mathrm{F}}\mathrm{\left(}\overline{\mathit{\theta }}\mathrm{\right)}$.

In the text

	Fig. 10 Average gravitational-darkening slope ⟨ B1 ⟩ given by Eq. (40) as function of η, calculated respectively for: a) Υ(θ) = coskθ; b) sinkθ. The colours identify the differential rotation parameter α, while k is identified using different line-styles.
In the text

Fig. 11 Relation between the apparent 4ln [ Teff(i) / ⟨ T ⟩ ] and ln [ geff(i) / ⟩ ⟨ g ⟩ ]; a)–c) for Υ(θ) = coskθ with α from − 0.5 to +0.5; d)–f) for sinkθ and α′ =− α/ (1 = α). All curves are for η = 0.8, with k = 0.5, 2.0, 6.0. The colour-codes according to α and α′ indicated in the respective blocks for k = 2 hold for all curves having the same Υ(θ). The equator-on seeing directions (i → π/ 2) are in the lower left corner of each block, while the polar-on directions (i → 0) are in the upper right corner.

In the text

	Fig. 12 Gravitational-darkening slope ♭1(i) against the rotational ratio η. The colours identify the values of the differential parameter α, and the line-styles are for the inclination angles. In a) to c) the slope ♭1(i) was calculated with Υ(θ) = coskθ and from d) to f) b1(i) was obtained with Υ(θ) = sinkθ. As throughout in this work, rigid rotation (α = 0) is identified with green lines.
In the text

Fig. 13 Comparison of the empirical gravitational-darkening slope ${\mathit{♭}}_{\mathrm{1}}^{\mathrm{obs}}$ (see reference in Domiciano de Souza et al. 2014) with the curves representing the ♭1(i) of rotators with k = 2, Υ(θ) = cosθ, α = − 0.5, 0.0 (green curves) and α = + 0.5, and inclination angles i = 0°, 45° and 90° (dentified with different line styles). The dashed black dashed line is for ♭1(i) = 1.

In the text