Gravity darkening in stars with surface differential rotation

TL;DR

This paper develops a model to analyze how surface differential rotation affects gravity darkening and effective temperature distribution in stars, emphasizing the importance of considering these effects for accurate stellar parameter interpretation.

Contribution

It introduces a comprehensive model linking surface differential rotation laws with gravity darkening, providing analytical solutions and highlighting the dependence of key parameters on rotation and viewing angle.

Findings

Gravity darkening exponent depends on rotation and viewing angle.

Analytical forms derived for specific differential rotation parameters.

Recommendations for analyzing stellar data considering differential rotation effects.

Abstract

The interpretation of stellar apparent fundamental parameters must be treated consistently with the characteristics of their surface rotation law. We develop a model to determine the distribution of the effective temperature and gravity, which depend on the surface differential rotation law and on the stellar external geometry. The basic assumptions are: a) the atmosphric layers are in radiative equilibrium; b) the bolometric flux is anti-parallel with the effective gravity; c) the angular velocity in the surface is Omega=Omega(theta,alpha,k)] where (alpha,k) are free parameters. The effective temperature varies with co-latitude `theta' and on the differential-rotation law through the effective gravity and the gravity-darkening function (GDF). Although the derived expressions can be treated numerically, for some low integer values of k, analytical forms of the `integral ofcharacteristic curves', on which the determination of the GDF relies, are obtained. The effects of the quantities (eta,alpha,k) (eta=ratio between centrifugal and gravitational accelerations at the equator) on the determination of the Vsini parameter and on the `gravity-darkening exponent' are studied. Depending on the values of (eta,alpha,k) the velocity V in the derived Vsini may strongly deviate from the equatorial rotational velocity. It is shown that the von Zeipel's-like gravity-darkening exponent $\beta_1$ depends on all parameters $(\eta,\alpha,k)$ and that its value also depends on the viewing-angle i. Hence, there no unique interpretation of this exponent determined empirically in terms of (i,alpha). The data on rotating stars should be analyzed by taking into account the rotational effects through the GDF, by assuming k=2 as a first approximation. Instead of the classical pair (eta,beta1), it would be more useful to determine the quantities (eta,alpha,i) to characterize stellar rotation.

Tables

Table 1. Table 1: Equatorial to polar radii ratio R e / R p subscript 𝑅 e subscript 𝑅 p R_{\rm e}/R_{\rm p} as a function of η 𝜂 \eta , k 𝑘 k and α 𝛼 \alpha , for Υ = cos k ⁡ θ Υ superscript 𝑘 𝜃 \Upsilon=\cos^{k}\theta and Υ = sin k ⁡ θ Υ superscript 𝑘 𝜃 \Upsilon=\sin^{k}\theta .

	${\rm Y}={\cos }^k\theta$				${\rm Y}={\sin }^k\theta$					${\rm Y}={\sin }^k\theta$
$𝜼=0.2$
$𝒌$	0.5	2.0	3.0	6.0	0.5	2.0	3.0	6.0		0.5	2.0	3.0	6.0
$𝜶=0.00$	1.100	1.100	1.000	1.000	1.100	1.100	1.000	1.000	${𝜶}^{\prime }=0.0$	1.100	1.100	1.000	1.000
-0.50	1.037	1.061	1.069	1.081	1.142	1.218	1.246	1.293	+1.0	1.083	1.061	1.053	1.043
-0.25	1.065	1.078	1.083	1.090	1.113	1.134	1.142	1.154	+1/3	1.091	1.078	1.074	1.068
+0.25	1.143	1.125	1.120	1.112	1.093	1.082	1.079	1.074	-1/5	1.109	1.125	1.130	1.139
+0.50	1.192	1.153	1.141	1.125	1.088	1.072	1.067	1.058	-1/3	1.120	1.153	1.165	1.185
$𝜼=0.8$
$𝒌$	0.5	2.0	3.0	6.0	0.5	2.0	3.0	6.0		0.5	2.0	3.0	6.0
$𝜶=0.00$	1.400	1.400	1.400	1.400	1.400	1.400	1.400	1.400	${𝜶}^{\prime }=0.00$	1.400	1.400	1.400	1.400
-0.50	1.159	1.274	1.306	1.347	1.515	1.719	1.794	1.918	+1.0	1.350	1.274	1.247	1.202
-0.25	1.268	1.333	1.350	1.372	1.437	1.499	1.523	1.564	+1/3	1.374	1.333	1.317	1.292
+0.25	1.545	1.474	1.455	1.431	1.379	1.346	1.333	1.312	-1/5	1.427	1.474	1.491	1.521
+0.50	1.693	1.552	1.515	1.466	1.366	1.312	1.292	1.259	-1/3	1.456	1.552	1.589	1.651
Note: ${\alpha }^{\prime }=-\alpha /\left(1+\alpha \right)$; $R_e\left(\alpha =0\right)/R_p\left(\alpha =0\right)=1+\eta /2$ given by the Roche gravity-rotation potential that identifies the surface of
rigid rotators. The ${\rm Y}={\sin }^k\theta$ function is used first with $\alpha$ and then with ${\alpha }^{\prime }$.

Equations

$$\displaystyle\begin{array}[]{l}\displaystyle F=c(\Phi)g^{\beta_{1}}_{\rm eff},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\@vec{F}=-\mathscr{F}({\rm r},\theta)\@vec{g}_{\rm eff},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\nabla\cdot\@vec{F}=\@vec{g}_{\rm eff}\cdot\nabla\mathscr{F}+\mathscr{F}\nabla\cdot\@vec{g}_{\rm eff}=0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\Omega(\theta)=\Omega_{\rm o}[1+\alpha\Upsilon(\theta)],\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\displaystyle\Upsilon(\theta)&=&\left\{\begin{array}[]{ll}\displaystyle\cos^{k}\theta\\ \displaystyle\sin^{k}\theta,\end{array}\right.\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\displaystyle\sin^{k^{\prime}}\theta&=&1-\cos^{k}\theta\\ \displaystyle k^{\prime}&=&k^{\prime}(k,\theta),\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\eta=\displaystyle\Omega^{2}_{\rm e}R^{3}_{\rm e}/GM=(\Omega_{\rm e}/\Omega_{\rm c})^{2}\!(R_{\rm e}/R_{\rm c})^{3},\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\displaystyle\@vec{g}_{\rm eff}&=&\displaystyle g_{\rm r}\,\hat{\@vec{e}}_{\rm r}+g_{\theta}\,\hat{\@vec{e}}_{\theta}\\ \displaystyle g_{\rm r}&=&\displaystyle-\frac{GM}{{\rm r}^{2}}+\Omega^{2}{\rm r}\sin^{2}\theta\\ \displaystyle g_{\theta}&=&\displaystyle\Omega^{2}{\rm r}\sin\theta\cos\theta\,,\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\displaystyle g_{\rm eff}&=&\displaystyle\left(g^{2}_{\rm r}+g^{2}_{\theta}\right)^{1/2}=\langle{g}\rangle\gamma(r,\theta)\\ \displaystyle\gamma(r,\theta)&=&\displaystyle\frac{1}{r^{2}}\left\{1+[(1-\tilde{\eta}r^{3})^{2}-1]\sin^{2}\theta\right\}^{1/2}\\ \displaystyle\tilde{\eta}&=&\displaystyle\eta\left[\frac{\Omega(\theta)}{\Omega_{\rm e}}\right]^{2}\\ \displaystyle\langle{g}\rangle&=&\displaystyle g_{o}\left(\frac{R_{o}}{R_{\rm e}}\right)^{2},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\lim_{r\to 0}\mathscr{F}({\rm r},\theta)=L/4\pi GM,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle f({\rm r},\theta)=\mathscr{F}({\rm r},\theta)\frac{4\pi GM}{L}.\\ \end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\lim_{r\to 0}f(r,\theta)=1.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle P\frac{\partial f}{\partial r}+Q\frac{\partial f}{\partial\theta}=2f,\end{array}$$

$$\displaystyle\begin{array}[]{rcl}\displaystyle P&=&\displaystyle-\frac{g_{r}}{\Omega\left[\Omega+(\partial\Omega/\partial\theta)\sin\theta\cos\theta\right]}\\ \displaystyle Q&=&\displaystyle-\frac{g_{\theta}}{r\Omega\left[\Omega+(\partial\Omega/\partial\theta)\sin\theta\cos\theta\right]}.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\ln f=\ln G+A(\theta)\;.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle A(\theta)=-\ln\left[\Omega^{2}(\theta)\tan^{2}\theta\right].\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle P\frac{\partial G}{\partial r}+Q\frac{\partial G}{\partial\theta}=0.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle Qdr-Pd\theta=0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle-(g_{\theta}/r)dr+g_{r}d\theta=0.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle M(r,\theta)=N(\theta)\frac{\Omega^{2}(\theta)r^{2}}{GM}\;,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\frac{\partial(-Mg_{\theta}/r)}{\partial\theta}=\frac{\partial(Mg_{r})}{\partial r},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle N(\theta)=A\frac{\cos^{2}\theta}{\Omega^{4}(\theta)\sin\theta},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\frac{\partial\tau}{\partial{r}}=-\frac{Mg_{\theta}}{r}=\frac{r^{2}\Omega^{2}_{\rm o}}{GM}\cos^{3}\theta\\ \displaystyle\frac{\partial\tau}{\partial\theta}=Mg_{r}=\frac{\cos\theta\cot\theta}{\left[1+\alpha\Upsilon(\theta)\right]^{2}}-\frac{r^{3}\Omega^{2}_{\rm o}}{GM}\sin\theta\cos^{2}\theta,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\tau(r,\theta)=\frac{1}{3}\eta\left(\frac{\Omega_{\rm o}}{\Omega_{\rm e}}\right)^{2}r^{3}\cos^{3}\theta+\mathscr{T}(\theta)\;.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\mathscr{T}(\theta)=\int_{\pi/2}^{\theta}\frac{\cos\phi\cot\phi}{\left[1+\alpha\Upsilon(\phi)\right]^{2}}d\,\phi\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\lim_{\alpha\to 0}\mathscr{T}(\theta)=\ln\tan\frac{\theta}{2}+\cos\theta,\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle Q\left(\frac{\partial f}{\partial\theta}\right)_{\tau}=2f\ \ \ {\rm and}\ \ \ P\left(\frac{\partial f}{\partial r}\right)_{\tau}=2f.\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle f(r,\theta)=\frac{\Psi(\tau)}{\tan^{2}\theta[1+\alpha\Upsilon(\theta)]^{2}},\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\mathscr{T}(\vartheta)=\frac{1}{3}\eta\left(\frac{\Omega_{\rm o}}{\Omega_{\rm e}}\right)^{2}r^{3}\cos^{3}\theta+\mathscr{T}(\theta),\end{array}$$

$$\displaystyle\begin{array}[]{l}\displaystyle\Psi(\tau)=\tan^{2}\vartheta\,\left[1+\alpha\Upsilon(\vartheta)\right]^{2}.\end{array}$$

Full text

*Aims. *We aim to develop a model to determine the distribution of the effective temperature and gravity, which explicitly depend on the surface differential rotation law and on the concomitant stellar external geometry.

*Methods. *The basic assumptions in this model are: a) the external stellar layers are in radiative equilibrium; b) the emergent bolometric flux is anti-parallel with the effective gravity; c) the angular velocity in the surface obeys relations like $\Omega(\theta)=\Omega_{\rm o}[1+\alpha\Upsilon(\theta,k)]$ where $\Upsilon(\theta,k)=\cos^{k}\theta$ or $\sin^{k}\theta$, and where $(\alpha,k)$ are free parameters.

*Results. *The effective temperature varies with co-latitude $\theta$, with amplitudes that depend on the differential-rotation law through the surface effective gravity and the gravity-darkening function (GDF). Although the derived expressions can be treated numerically, for some low integer values of $k$, analytical forms of the “integral of characteristic curves”, on which the determination of the GDF relies, are obtained. The effects of the quantities $(\eta,\alpha,k)$ ($\eta=$ ratio between centrifugal and gravitational accelerations at the equator) on the determination of the $V\!\sin i$ parameter and on the “gravity-darkening exponent” are studied. Depending on the values of $(\eta,\alpha,k)$ the velocity $V$ in the derived $V\!\sin i$ may strongly deviate from the equatorial rotational velocity. It is shown that the von Zeipel’s-like gravity-darkening exponent $\beta_{1}$ depends on all parameters $(\eta,\alpha,k)$ and that its value also depends on the viewing-angle $i$. Hence, there no unique interpretation of this exponent determined empirically in terms of $(i,\alpha$).

*Conclusions. *We stress that the data on rotating stars should be analyzed by taking into account the rotational effects through the GDF, by assuming $k=2$ as a first approximation. Instead of the classic pair $(\eta,\beta_{1})$, it would be more useful to determine the quantities ($\eta,\alpha,i$) to characterize stellar rotation.

††offprints: J. Zorec:

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 $\varpi$ to the rotational axis, $\Omega=\Omega(\varpi)$; 2) the effective gravity is derived from a total gravitation-rotational potential, $\Phi$; 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_{\rm eff}$, known as von Zeipel’s theorem (von Zeipel, 1924)

[TABLE]

where $\Phi$ is the total gravitation-rotational potential describing the stellar surface and $c(\Phi)$ 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 $\beta_{1}=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, $\Omega=\Omega(\varpi,z)$ ($z$ coordinate parallel to the rotation axis), authors have shown that $\beta_{1}=\beta_{1}(\theta)\leq 1$, where $\theta$ 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. (2) 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 $\@vec{F}$ emitted by an axially symmetric stellar atmosphere in radiative equilibrium is anti-parallel to the vector of the local effective gravity $\@vec{g}_{\rm eff}$. Thus, assuming this parallelism, they write

[TABLE]

where (${\rm r},\theta$) are the spherical coordinates. The function $\mathscr{F}({\rm r},\theta)$ is hereinafter called the gravity-darkening function (GDF). We note that before Espinosa Lara & Rieutord (2011), $\mathscr{F}({\rm r},\theta)$ was simply the constant $c(\Phi)$ in Eq. (2). 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 $\mathscr{F}({\rm r},\theta)$ 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

[TABLE]

and obtained that $\beta_{1}$ depends on the co-latitude angle $\theta$ and that it is everywhere a decreasing function of the the stellar flattening, $\varepsilon=1-R_{\rm p}/R_{\rm e}$ ($R_{\rm p}$ and $R_{\rm 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 $\beta^{\rm obs}_{1}\lesssim 1$ (Monnier et al., 2007; Zhao et al., 2009; Che et al., 2011; Monnier et al., 2012; Domiciano de Souza et al., 2014), although within the current measurement uncertainties, these values are $\beta^{\rm obs}_{1}\gtrsim\overline{\beta_{1}}(\varepsilon)$ systematically, where $\overline{\beta_{1}}(\varepsilon)$ is an average relation between $\beta_{1}$ and $\varepsilon$ [see Eq. (90) in Sect.7.1].

Apart from the above mentioned radiation transfer reasons, the exponent $\beta_{1}$ determined observationally can also be dependent on: 1) the line-of-sight angle $i$ of stars, because $\beta_{1}$ is a function of $\theta$; 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 $\beta_{1}$ 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 $\mathscr{F}({\rm r},\theta)$ by imposing the condition of radiative equilibrium in the stellar atmosphere (cf. Eq. 6), while the stellar surface rotates differentially. We explore the incidence of this rotation on the observed stellar fundamental parameters, in particular the measured $V\!\sin i$ parameter and the so-called “gravity-darkening exponent” $\beta_{1}$.

2 The surface differential rotation law

Having no a priori indications on a general form of the function describing the surface differential rotation law $\Omega(\theta)$, 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

[TABLE]

where $\alpha$ is a free parameter called the “differential rotation parameter”. The case $\alpha=0$ corresponds to rigid rotation. Accordingly, we adopt two “exploratory” forms for the function $\Upsilon(\theta)$, namely

[TABLE]

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. (8) and (12) for $\alpha=+0.5$ ($\alpha^{\prime}\!=\!-1/3$), $\alpha=0.0$ and $-0.5$ ($\alpha^{\prime}\!=\!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 $\Upsilon(\theta)=\cos^{k}\theta$ and near the pole for $\Upsilon(\theta)=\sin^{k}\theta$. When $\Upsilon(\theta)=\cos^{k}\theta$ we have $\Omega_{\rm o}=\Omega_{\rm e}$, while for $\Upsilon(\theta)=\sin^{k}\theta$, $\Omega_{\rm o}=\Omega_{\rm e}(1+\alpha)^{-1}$. In both cases $\Omega_{\rm e}$ is the angular velocity at the equator.

In order to grasp the differences inherent to the use of one of the two $\Upsilon$-functions given in Eq. 12, it is necessary to compare the effects produced by rotation laws implying the same ratio $\Omega_{\rm p}/\Omega_{\rm e}$. Thus, for a given value of $\alpha$ used with $\Upsilon(\theta)=\cos^{k}\theta$, we have to take $\alpha^{\prime}=-\alpha/(1+\alpha)$ associated with $\Upsilon(\theta)=\sin^{k}\theta$ so that $\Omega(\theta)=\Omega_{\rm o}[1+\alpha^{\prime}\sin^{k}\theta]$. Obviously, when $k=2$ both $\Upsilon(\theta)$ represent the same rotation law. Nevertheless, keeping the reciprocity between $\alpha$ and $\alpha^{\prime}$, the functions $\Upsilon(\theta)=\cos^{k}\theta$ and $\Upsilon(\theta)=\sin^{k^{\prime}}\theta$ can produce resembling behaviours of predicted parameters (although not exactly the same), because

[TABLE]

independently of the value of $\alpha$. Taking the average of the two thirds of powers $k^{\prime}(k,\theta)$ calculated in the middle points of the $\theta-$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 $\Upsilon(\theta)=\cos^{k}\theta$ with $k=0.5$ will lead to resembling results as $\Upsilon(\theta)=\sin^{k^{\prime}}\theta$ when $k^{\prime}\simeq 7.0$, or $\cos^{k}\theta$ with $k=6.0$ similar results as $\sin^{k^{\prime}}\theta$ where $k^{\prime}\simeq 0.5$ [this effect is illustrated in many figures below, as: Figs. 4a and 4f, or Figs. 4c and 4d; Figs. 5a and 5f, or Figs. 5c and 5d; Figs. 6a,d and 6b,d, or Figs. 7c,f and 7a,c.; Figs. 12a and 12f, or Figs. 12c and 12d].

In principle we could let $-\infty<\alpha<+\infty$, but we shall limit the values of the differential rotation parameter to $-1\leq\alpha<\infty$ 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 $\alpha$ to the interval $-0.5<\alpha\leq+0.5$. We note that the solar surface differential rotation can be sketched with $\alpha\simeq-0.3$ and $\Upsilon(\theta)=\cos^{2}\theta$. On the other hand, Espinosa Lara & Rieutord (2013) predicted $\alpha\gtrsim-0.2$ for stars with masses $2\lesssim M/M_{\odot}\lesssim 4$.

Instead of using $\Omega_{\rm e}/\Omega_{\rm c}$ to characterize the rotational velocity at the equator, in this paper we use the ratio $\eta$ between the centrifugal and gravitational acceleration at the equator defined as follows

[TABLE]

where $M$ is the stellar mass, $G$ the graviational constant; $\Omega_{\rm c}$ and $R_{\rm 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_{\rm e}/V_{\rm K}$ is used as the non-dimensional parameter characterizing the stellar rotation ($V_{\rm e}$ and $V_{\rm K}$ are the actual and Keplerian linear equatorial velocities, respectively), we note that $\eta=(V_{\rm e}/V_{\rm K})^{2}$, so that it holds $0\leq\eta\leq V_{\rm e}/V_{\rm K}$ $\leq\Omega_{\rm e}/\Omega_{\rm c}\leq 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. 8 represents the boundary condition in the surface of an internal non-conservative rotational law $\Omega\!=\!\Omega(r,\theta)$. 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 $d$s does not imply any work done by the effective gravity geff, i.e. geff.$d$s=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. (8) and (12), and taking into account several values of $\alpha$ and $k$. The rotation at the equator is characterized by the rotational ratio $\eta=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 $\alpha$ (and the corresponding $\alpha^{\prime}$). The green curves are for $\alpha=0$.

In Table 1 we give the radii ratios $R_{\rm e}/R_{\rm p}$ for $\eta=0.2$ ($\Omega/\Omega_{\rm c}=$ 0.69; $V/V_{\rm c}=$ 0.52) and $\eta=0.8$ ($\Omega/\Omega_{\rm c}=$ 0.99; $V/V_{\rm c}=$ 0.92), for several values of $k$ and $\alpha$. In each series of parameters ($\eta,k,\alpha$) the radii ratios were calculated also for $\Upsilon=\sin^{k}\theta$ and $\alpha^{\prime}=-\alpha/(1+\alpha)$.

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 $\kappa=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 $\@vec{g}_{\rm eff}$ and write

[TABLE]

where $G$ is the gravitational constant; $M$ is the stellar mass; $\Omega=\Omega(\theta)$ is the angular velocity given by Eqs. (8) and (12); $\hat{\@vec{e}}_{\rm r}$ and $\hat{\@vec{e}}_{\theta}$ are the unit vectors associated with the spherical coordinates $({\rm r},\theta)$. In what follows, we use the dimensionless radial coordinate $r={\rm r}/R_{\rm e}$, where $R_{\rm e}$ is the rotationally modified stellar equatorial radius. To simplify notation we write $g_{\rm eff}=|-\@vec{g}_{\rm eff}|$ and introduce the dimensionless expression for the effective gravity $\gamma(r,\theta)$ as follows

[TABLE]

where $g_{o}=GM/R^{2}_{o}$ 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. (26) it follows that the extreme values of $\gamma(r,\theta)$ are: $\gamma_{\rm e}=\gamma(r,\pi/2)=$ $1-\eta$ and $\gamma_{\rm p}=\gamma(r,0)=$ $(R_{\rm e}/R_{\rm p})^{2}$, which are functions of $(\eta,\alpha)$ and are dependent on the function $\Upsilon(\theta)$ adopted (see Table 1).

4 The gravity-darkening function

The gravity darkening function $\mathscr{F}({\rm r},\theta)$, as defined by Eq. (4), is constrained by the boundary condition obtained when ${\rm r}\to 0$, namely

[TABLE]

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; Ekström et al., 2008; Maeder, 2009; Espinosa Lara & Rieutord, 2013; Fujisawa, 2015). In what follows, we use the dimensionless form of $\mathscr{F}$

[TABLE]

From Eq. (28) we have

[TABLE]

Using spherical coordinates, the relation in Eq. (6) becomes the following linear partial differential equation:

[TABLE]

where the notations stand for

[TABLE]

The integration of Eq. (34) 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. (34) by writing the equation for $\ln f$ and setting :

[TABLE]

After some calculations we find that

[TABLE]

Now, Eq. (34) may be rewritten

[TABLE]

The characteristic lines of $G$ in the $(r,\theta)$-plane, are the lines where $G$ is constant. On these lines $\partial_{r}Gdr+\partial_{\theta}Gd\theta=0$, so that on these lines we have

[TABLE]

which is the equation of characteristic lines. After some simplifications it reads

[TABLE]

To solve (47), we need to transform it into an exact differential. We can indeed multiply the equation by any function of $(r,\theta)$. It turns out that if we multiply (47) by

[TABLE]

it can be integrated. Indeed, demanding that

[TABLE]

leads to the following expression for $N(\theta)$:

[TABLE]

where $A$ is an arbitrary constant.

Thus $G$ is constant on the lines $\tau(r,\theta)=$ Cst, if $\tau$ verifies:

[TABLE]

where we chose $A=-\Omega^{2}_{\rm o}$. Since the RHS of the previous equations have been chosen such that $d\tau$ is an exact differential, the integration is straightforward and gives

[TABLE]

The function $\mathscr{T}(\theta)$, here referred to as the integral of characteristic curves, is finally given by

[TABLE]

that can be calculated numerically for any real positive value of $k$. For some low natural number $k$, the integral in Eq. (60) admits analytical expressions that are given in appendix A.

In this work we used only numerical estimates of $\mathscr{T}(\theta)$. The numerical estimates of the integral in Eq. (60) were obtained proceeding to the change of variable $u=\ln\tan(\phi/2)$ and using the 16-point Gauss-Legendre quadrature rule. The function $\mathscr{T}(\theta)$ 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. (60) are $\delta\mathscr{T}\lesssim 10^{-6}$ at $\theta\sim 10^{-5}$ rad ($\theta\sim 0^{\circ}.001$), and they become $\delta\mathscr{T}\lesssim 10^{-12}$ as soon as $\theta\gtrsim 0.07$ rad ($\theta\gtrsim 8^{\circ}$). The precision of the numerical estimates increases even more as $\theta\to\pi/2$.

The limit expression of $\mathscr{T}(\theta)$ for $\alpha=0$ comes directly from Eq. (60), which is also obeyed by all analytical expressions given in Appendix A:

[TABLE]

which is the form of $\mathscr{T}(\theta)$ obtained by Espinosa Lara & Rieutord (2011) [see their Eq. (20)].

Because $\tau\equiv\tau(r,\theta)$, where $r$ and $\theta$ are independent variables, we can use instead the pairs $(r,\tau)$ or $(\theta,\tau)$ as independent variables to obtain two simplified forms of Eq. (34) valid over the $\tau-$characteristic curves

[TABLE]

The first differential equation in Eq. (64) can be easily integrated to obtain

[TABLE]

where the factor $\Psi(\tau)$ is an integration constant that carries the dependence of $f$ with $r$ through a transcendental function $\vartheta(r,\theta)$ defined, using Eq. (58), as

[TABLE]

which assumes that there is a unique relation between $\vartheta$ and $\tau$. From Eq. (68) it appears that $\vartheta\to\theta$ as $r\to 0$. Therefore, noting that $\tau(r,\theta)$ is constant over a characteristic curve, the value of $\Psi(\tau)$ can be specified at any point $(r,\theta)$ on the characteristic, in particular at $r=0$. Thus, making use of the boundary condition in Eq. (32) we obtain

[TABLE]

Introduced into Eq. (66), it leads to

[TABLE]

Eq. (72) is the generalized version of that obtained by Espinosa Lara & Rieutord (2011), which is valid for $\alpha=0$ [their Eq. (26)].

Equations (58) and (68), together with Eq. (60), or the forms for $\mathscr{T}(\theta)$ given in Appendix A, define $\vartheta(r,\theta)$ for every chosen point $(r,\theta)$, in particular for the stellar surface represented by the function $r=r(\theta)$. To this end it is worth noting that because $\vartheta\neq\theta$ except at the extremes $\vartheta=\theta=0$ and $\vartheta=\theta=\pi/2$, it is impossible to derive reliable values of $\vartheta=\vartheta(r,\theta)$ from Eq. (68) if the stellar geometry, i.e. $r=r(\theta)$, is not determined consistently with the rotation law on which depends $\mathscr{T}(\theta)$.

We have calculated GDF numerically using Eq. (72), where $r(\theta)$ is the radius-vector that describes the stellar surface. The function $\vartheta(\theta)$ is determined with Eqs. (58) and (68) by interpolation, where $\mathscr{T}(\vartheta)$ at a given $\theta$ is taken as a linear function of $\ln\theta$. For values of $\vartheta$ when $\theta\to\pi/2$ we consider $\mathscr{T}(\vartheta)$ as a linear function of its expression for $\alpha=0$, so that $\vartheta$ is derived by iterating Eq. (62).

The behaviour of the function $f(\theta)$ calculated for $\Upsilon=\cos^{k}\theta$ $\eta=0.2$, 0.8, and for several values of $\alpha$ and $k$ is shown in Fig. 3. In this figure the function $f(\theta)$ calculated with $\Upsilon=\sin^{k}\theta$ and $\alpha^{\prime}=$ $-\alpha/(1+\alpha)$ is also shown. It shows that $f[\cos^{k}\theta,\alpha]\neq$ $f[\sin^{k}\theta,\alpha^{\prime}(\alpha)]$, except for $k=2$ where $f(\theta)$ is the same for both $\Upsilon(\theta)$ functions. The function $f(\theta)$ obtained for $\alpha=0.0$ (green curves in Fig. 3) corresponds to the solution previously obtained by Espinosa Lara & Rieutord (2011).

From Eq. (LABEL:char) we can derive the limits of $f(r,\theta)$ when $\theta\to 0$ and $\theta\to\pi/2$

[TABLE]

where $R_{\rm p}$ and $R_{\rm e}$ are the polar and equatorial radii, respectively. $\Omega_{\rm p}$ is the polar angular velocity, which imply that $\Omega_{\rm p}/\Omega_{\rm e}=1+\alpha$ when $\Upsilon(\theta)=\cos^{k}\theta$ and $\Omega_{\rm p}/\Omega_{\rm e}=(1+\alpha^{\prime})^{-1}=(1+\alpha)$ for $\Upsilon(\theta)=\sin^{k}\theta$. The dependence of $f^{\rm pole}$ with $\alpha$ and $k$ is given through the radii ratio $R_{\rm p}/R_{\rm e}$, which is a function of them both and on the function $\Upsilon(\theta)$ chosen. However, $f^{\rm eq}$ depends only on the ratio $\eta$. It may be worth noting that for $\eta=0$, it is $\vartheta\equiv\theta$ so that in Eq. (72) the GDF becomes $f(r,\theta)=1$. This also makes that from Eq. (72) we have $f^{\rm pole}=f^{\rm eq}=1$.

For a given ratio $\eta$ we obtain the same value $f^{\rm eq}$ whatever the function $\Upsilon(\theta)$. We may wonder then whether the same value for $f^{\rm pole}$ can be derived whatever the chosen function $\Upsilon(\theta)$ provided we use the appropriate reciprocal gravity-darkening parameter, i.e. $\alpha$ or $\alpha^{\prime}=-\alpha/(1+\alpha)$, as it is apparent form the expression of $f^{\rm pole}$ given in Eq. (75). Actually, this reciprocity is valid only for $k=2$. In the next section we show that for $k\neq 2$ this apparent reciprocity does not produce the same ratio $R_{\rm p}/R_{\rm e}$, and the geometry of the stellar surface is not the same either. This is due to the fact that $R_{\rm e}/R_{\rm p}=1+\eta\Gamma$, with $\Gamma=\int_{0}^{\pi/2}[\Omega(\theta)/\Omega_{\rm e}]^{2}...d\theta$ [see Zorec et al. (2011)] where, as seen from Fig. 1, $\Omega(\theta)/\Omega_{\rm e}$ depends on $\Upsilon(\theta)$ for $k\neq 2$, even if we use $\alpha$ and its reciprocal $\alpha^{\prime}$ according to each case. Consequently, the behaviour of the resulting GDF, $f(\theta)$, over the interval $0<\theta<\pi/2$ will depend in a strict sense on the chosen $\Upsilon(\theta)$ function if $k\neq 2$.

5 The effective temperature as a function of the colatitude

The magnitude of the flux-vector in Eq. (4) can be rewritten as

[TABLE]

where $\sigma_{\rm SB}$ is the Štefan-Boltzmann constant. Using Eqs. (26) and (30) this relation can be translated into

[TABLE]

where the radial variable $r=r(\theta)$ was used as the radius-vector representing the stellar surface, so that $\gamma$ becomes a function of $\theta$ only, apart from its dependence on ($\eta,\alpha$). It is obvious that for $\eta=0$ (no rotation), $f(\theta)=\gamma(\theta)=1$ and $T_{\rm eff}(\theta)=\langle{T}\rangle$. The factor $\langle{T}\rangle$ is the uniform effective temperature over a sphere of radius $R_{\rm e}(\eta)$ radiating the luminosity $L$

[TABLE]

It must be understood that when interpreting observations the quantity $\langle{T}\rangle$ 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. 4(a) to (c) $T_{\rm eff}(\theta,\eta,\alpha)/\langle T\rangle$ is shown as a function of the co-latitude $\theta$ when $\Upsilon(\theta)=\cos^{k}\theta$, $\eta=0.8$, $k=0.5$, 2.0 and 6.0, and for several values of $\alpha$ ranging from $-0.5$ to $+0.5$. Figures 4(d) to (f) show $T_{\rm eff}(\theta,\eta,\alpha)/\langle T\rangle$ against $\theta$ for $\Upsilon(\theta)=\sin^{k}\theta$ and ($k,\alpha^{\prime}$), where $\alpha^{\prime}=-\alpha/(1+\alpha)$. The green curves correspond to $T_{\rm eff}(\theta)/\langle{T}\rangle$ for $\alpha=\alpha^{\prime}=0$, previously obtained for rigid rotators by Espinosa Lara & Rieutord (2011). In Fig. 5 we plot the effective temperature ratios $T^{\rm eq}_{\rm eff}/T^{\rm p}_{\rm eff}$ as a function of $\eta$ and different values $k$ and $\alpha$, when $\Upsilon(\theta)=\cos^{k}\theta$ or $\Upsilon(\theta)=\sin^{k}\theta$.

6 The $V\!\sin i$ parameter

$V\!\sin i$ 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. (8) and both forms of $\Upsilon(\theta)$ given in Eq. (12). For this exercise we used a model star with mass $M=9M_{\odot}$ and fractional age $t/t_{\rm MS}=0.6$ ($t_{\rm MS}$ is the time a rotating star spends in the main sequence evolutionary phase). Rotation laws are characterized by several parameters $k$, $\alpha$, and different equatorial velocities parametrized according to $\eta$. Applying the classic Fourier transform method we determined the $V\!\sin i$ parameters shown in Figs. 6 and 7 where we compare the resulting $V_{\rm eq}\!\sin i$ with the actual input $V_{\rm eq}\!\sin i$ corresponding to the given values of $\eta$ 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_{\rm R}(\lambda)$ (Gray, 1992).

Because it has become of common use today, the classic method based on the FT was employed to determine the $V\!\sin i$. This method assumes, however, that the observed rotationally broadened line profile $\mathscr{F}(\lambda)$ can be represented as the convolution with an analytical “rotation broadening-function” $G_{\rm R}(\lambda)$ of a flux line-profile $F(\lambda)$ emitted by a non-rotating star having uniform effective temperature and gravity. The specific intensity $I(\lambda,\mu)$ contributing to the observed line flux $F(\lambda)$ 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_{\rm R}(\lambda)$ 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 $V\!\sin i$ is then obtained by comparing the zeros of the flux line-profile FT $\mathscr{F}(\lambda)$ with the corresponding zeros of the FT of the function $G_{\rm R}(\lambda)$.

In Zorec et al. (2017) it has been shown that depending on the spectral line and on the value of $\alpha$, 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_{\rm R}(\lambda)$ 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(\lambda,\mu)$.

Figures 6 and 7 reveal that the $V\!\sin i$ obtained using the FT-method can be overestimated when $\alpha>0$, and underestimated for $\alpha<0$, as compared to its value for rigid rotators. We notice that for $\Upsilon(\theta)=\cos^{k}\theta$ and a given $|\alpha|$ the deviations are larger when $\alpha<0$ than for $\alpha>0$. For $\Upsilon(\theta)=\sin^{k}\theta$ the deviations reverse according to the $\alpha$-corresponding reciprocal parameter $\alpha^{\prime}$. Finally, according to the spectral line used, the obtained $V\!\sin i$ is not exactly the same. To apprehend what can be the possible deviations that the current $V\!\sin i$ determinations may incur, let us recall that in the Sun $\alpha\simeq-0.3$ and $\Upsilon(\theta)=\cos^{2}\theta$, and that predictions made by Espinosa Lara & Rieutord (2013) for rapid rotators of masses $2\lesssim M\lesssim 4M_{\odot}$ foresee $\alpha\gtrsim-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 $V\!\sin i$ obtained from different spectral lines. They could help us to derive indications on the value of the $\alpha$ parameter.

Needless to say that the classic method of determining $V\!\sin i$ 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.

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_{\rm eff}(\theta)$ over a rotating stellar surface. As noted in Sect.1, the known von Zeipel’s relation that holds for barotropic objects

[TABLE]

is characterized by the constant gravity-darkening exponent $\beta_{1}$, and the constant $C_{\rm F}$ which depends on the gravity-rotation equipotential of the stellar surface. The notation $\beta_{1}$ is used here to distinguish it from the exponent $\beta=\beta_{1}/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. (83).

7.1 Dependence of the gravity-darkening exponent with colatitude

Up to now, the relation in Eq. (83) has been used to interpret the observations of rapidly rotating stars regardless of whether they are barotropic or baroclinic. However, $\beta_{1}$ and $C_{\rm F}$ cannot be assumed independent of the colatitude-angle $\theta$ 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 $\beta_{1}$ and $C_{\rm F}$, Fig. 8 shows the relation between $\ln T_{\rm eff}(\theta)$ calculated with Eq. (79) against $\ln g_{\rm eff}(\theta)$ using both forms of $\Upsilon(\theta)$ in Eq. (12). The model-stars are assumed to rotate with $\eta=0.8$ with several values of $\alpha$ and $k$. Obviously, similar behaviours are noted for other values of $\eta$.

Because the curves in Fig. 8 are not straight lines, neither $C_{\rm F}$ nor $\beta_{1}$ can be considered as genuine constants, even though a quasi-linear relation between $\ln T_{\rm eff}(\theta)$ and $\ln g_{\rm eff}(\theta)$ seems to exist over some intervals of $\theta$. To emphasize more clearly the lack of constancy of $\beta_{1}$ and $C_{\rm F}$, we define a pseudo-exponent or “gravity-darkening slope” $b_{1}(\theta)$ as

[TABLE]

and consider that between two successive, close discrete points $\theta_{i}$ and $\theta_{i+1}$, $b_{1}(\overline{\theta})$ and $C_{\rm F}(\overline{\theta})$ can be assumed constant; i.e.:

[TABLE]

Inserting the values $b_{1}(\overline{\theta_{i}})$ in Eq. (85), the corresponding $\ln C_{\rm F}(\overline{\theta_{i}})$ can be estimated. The inferred quantities $b_{1}(\overline{\theta_{i}})$ and $\ln C_{\rm F}(\overline{\theta_{i}})$ are shown in Fig. 9. They clearly demonstrate that both $b_{1}$ and $C_{\rm F}$ are dependent on $\theta$ and that consequently they cannot be used as constants in Eq. (83).

7.2 The observed gravity-darkening slope

7.2.1 Neglecting the viewing angle

Replacing in Eq. (85) the expressions for $g_{\rm eff}$ and $T_{\rm eff}$ given respectively by Eqs. (26) and (79), and imposing that $b_{1}$ and $C_{\rm F}$ do not depend on $\theta$, it follows that

[TABLE]

where $f^{\rm pole}$ and $\gamma^{\rm pole}$ are the values of the functions $f(\theta)$ and $\gamma(\theta)$ in $\theta=0$. The function $B_{1}(\theta)$ behaves in a similar way as $b_{1}(\theta)$, so that it is needless to reproduce it graphically. Instead, we can replace $f(\theta)$ and $\gamma(\theta)$ in Eq. (88) by their values at $\theta=\pi/2$ [see Eq. (75)] to obtain a kind of average $\langle B_{1}(\theta)\rangle$ that takes the form

[TABLE]

which is a generalized version of $\langle B_{1}\rangle$ obtained by Espinosa Lara & Rieutord (2011) for rigid rotators. In Fig. 10, $\langle B_{1}\rangle$ is shown as a function of $\eta$ for $-0.5\leq\alpha\leq 0.5$, $k=0.5$, 2, 6, and for both functions $\Upsilon(\theta)$. The expression of $\langle B_{1}\rangle$ for $\alpha=0$ was sometimes used in the literature to compare the exponents $b_{1}$ derived from interferometric observations.

7.2.2 Considering the viewing angle

Any parameter describing a physical property of a rotating stellar atmosphere that depends on the co-latitude $\theta$, 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. (88). 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 $\langle{T}_{\rm eff}(i)\rangle$ and $\langle{g}_{\rm eff}(i)\rangle$, respectively.

Let us define $\langle{T}_{\rm eff}(i)\rangle$ and $\langle{g}_{\rm eff}(i)\rangle$ so as to obtain a slope that might be assimilated with an empirically derived gravity-darkening slope. The apparent effective temperature $\langle{T}_{\rm eff}(i)\rangle$ can be written as

[TABLE]

where $\mathcal{S}(i)$ 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)

[TABLE]

where $\hat{\rm\bf n}=$ $-\@vec{g}_{\rm eff}/g_{\rm eff}$ is the unit vector normal to the stellar surface; $\hat{\bm{\iota}}$ is the unit vector directed towards the observer; $\hat{\@vec{e}}_{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. $\mu$ is the directional cosine $\cos(\hat{\rm\bf n},\hat{\bm{\iota}})=$ $\hat{\rm\bf n}\cdot\hat{\bm{\iota}}$; $I(\mu)$ is the bolometric specific intensity of the radiation for which we use a quadratic form

[TABLE]

whose limb-darkening coefficients $\epsilon_{1,2}$ have been interpolated in the tables given by Claret (2000) for each pair of local fundamental parameters [$T_{\rm eff}(\theta),g_{\rm eff}(\theta)$]. Knowing that for each point over the stellar surface we have $|\@vec{F}|=F(\theta)$ $=2\pi\int_{0}^{1}I(\mu)\mu\,d\mu$, the intensity $I_{1}=I(\mu=1)$ can be written as

[TABLE]

which depends on $f(\theta)$ and $\gamma(\theta)$ through $F(\theta)$ given in Eq. (77).

Consistently with the definition of $\langle{T}_{\rm eff}(i)\rangle$, we can use Eq. (77) to estimate $\langle{g}_{\rm eff}(i)\rangle$ using the expression

[TABLE]

that with Eqs. (92) and (94) leads to

[TABLE]

where $\langle{f}\rangle_{i}$ is given by

[TABLE]

For the numerical estimates of the apparent $\langle{T}_{\rm eff}(i)\rangle$ and $\langle{g}_{\rm eff}(i)\rangle$ the knowledge of absolute values of $T_{\rm eff}(\theta)$ and $g_{\rm eff}(\theta)$ are required. To get a rough insight on their amplitude variation, we have assumed model-stars with mass $M=9M_{\odot}$ and fractional age $t/t_{\rm MS}=0.6$ ($t_{\rm 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_{\rm e}/R_{\rm o}$ ($R_{\rm 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[\langle{T}_{\rm eff}(i)\rangle/\langle{T}\rangle]$ and $\ln[\langle{g}_{\rm eff}(i)\rangle/\langle{g}\rangle]$ shown in Fig. 11 are calculated using $\Upsilon(\theta)=\cos^{k}\theta$ and $\sin^{k}\theta$, for $\eta=0.8$, $k=0.5$, 2.0, 6.0, and $-0.5\leq\alpha\leq-0.5$. While Fig. 8 represents the variation of $T_{\rm eff}$ against $g_{\rm 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 $\eta$, it turns out that the resolution among these curves according to the differential parameter $\alpha$ is a strong function of the power $k$ used in $\Upsilon(\theta)$.

From Fig. 11 it seems legitimate to assume relations of the type

[TABLE]

for each $[\Upsilon,k,\alpha]$ so as to represent the dependence of the slope $\flat(i)$ with the stellar inclination, as we did for $b_{1}(\theta)$ with $\theta$, i.e. by taking here $\flat(i)$ and $\langle A(i)\rangle$ as constants between two successive inclinations $i_{j}$ and $i_{j+1}$ in a given discretization of this variable in $0\degr\leq i\leq 90\degr$. The values of $\flat(i)$ obtained for three inclination angles: $i=0^{\circ}$, $45^{\circ}$ and $90^{\circ}$ against the ratio $\eta$ are shown in Fig. 12. Figures 12a,b,c are for $\Upsilon=\cos^{k}\theta$ with $k=0.5$, 2 and 6, while Fig 12d,e,f are for $\Upsilon=\sin^{k}\theta$ and the same exponents $k$. In these figures we note that the $\flat(i)$ curves have a global behaviour that is similar to $\langle{B_{1}}\rangle$ in Fig. 10, except for their spread as a function of $i$.

In Fig. 13 we show the values of $b_{1}^{\rm obs}$ determined with interferometric data (see Domiciano de Souza et al., 2014), where the interpretation of observations was done using Eq. (83). These points have been originally given as a function of the stellar flattening $\varepsilon=1-R_{\rm p}/R_{\rm e}$. We have converted $\varepsilon$ into $\eta$ considering, as the authors themselves did, that stars are rigid rotators, i.e. $R_{\rm e}/R_{\rm p}=1+\eta/2$, which imply that $\eta=2\varepsilon/(1-\varepsilon)$. In this diagram the curves $\flat_{1}(i)$ are reproduced for $\Upsilon(\theta)=\cos^{k}\theta$ with $k=2$ and different values of $\alpha$, in particular those suited for rigid rotators, i.e. $\alpha=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 $\flat_{1}(i=0)$ and $\flat_{1}(i=\pi/2)$ of rigid rotators. Two of them are likely situated in the region of $\alpha>0$, while there is one in the zone of $\alpha<0$, which indicates that other solutions than rigid rotation might also suit these objects.

However, according to the diagrams in Fig. 12 it clearly appears that there is not a unique way to characterize the studied stars with $\beta_{1}$ as rotators, since many solutions can be envisioned in terms of $[\Upsilon,k,\alpha]$. On the other hand, there is a strong concern with the interpretation of the observed values $\flat_{1}^{\rm obs}$, because $\beta_{1}$ and $C_{\rm F}$ were considered as constants in the reduction of data. Due to these ambiguities, it would perhaps be wiser to abandon the parameter $\beta_{1}$ to model fluxes and rather concentrate on determining them using the differential rotation parameter $\alpha$ 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. (8) and (12), 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 $V\!\sin i$ 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(\theta)$, which describes the distribution of the effective temperature as indicated by Eqs. (79), (72), (68) and (60), 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 $(\eta,\alpha,\Upsilon)$ instead of the simple pair $(\eta,\beta_{1})$.

In the approach where differential rotation is considered, a significant uncertainty remains in the choice of the function $\Upsilon(\theta)$. Indeed, up to now, and except for the Sun, no empirical indication exists on the actual nature of $\Upsilon(\theta)$ in stars. The quantity that needs to be inferred first is a parameter related to the average surface gradient $\langle{d\,\Omega/d\,\theta}\rangle$. To this end, the simplest expression for $\Upsilon(\theta)$ could be $2\theta/\pi$. This makes $\alpha$ the unique free parameter characterizing the surface differential rotation, but in this case the calculation of the integral $\mathscr{T}(\theta)$ [Eq. (60)] has some annoying numerical subtleties. However, we have shown in Sect. 2 that for the same ratio $\Omega_{\rm p}/\Omega_{\rm e}$, no significant differences in the prediction or analysis of observed parameters may arise if we use only $\Upsilon(\theta)=\cos^{k}\theta$ or $\sin^{k^{\prime}}\theta$, because reliable correspondences can be established between $(\alpha,k)$ and $(\alpha^{\prime},k^{\prime})$. According to this finding, we can adopt $\Upsilon(\theta)=\cos^{k}\theta$ and show easily that the weighted average gradient over the stellar surface

[TABLE]

is equal to $\langle d(\Omega/\Omega_{o})/d\theta\rangle=$ $2\alpha/\pi$ if $k=2.55$. Although nothing prevents us from making numerical models with $k=2.55$, the results for $T_{\rm eff}(\theta)$, $V\!\sin i$ and other parameters, show that no significant differences occur on the predicted or analyzed parameters if we simply adopt $k=2$. Moreover, the representation $\Upsilon(\theta)=\cos^{2}\theta$ 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 $\Upsilon(\theta,k)$, in particular for $k=2$, analytical expressions of $\mathscr{T}(\theta)$ exist and can help to proceed in a similar way as for rigid rotators with the formulas developed by Espinosa Lara & Rieutord (2011).

Acknowledgements.

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.

Appendix A The functions $\mathscr{T}(\theta)$

The analytical expressions obtained of the function $\mathscr{T}(\theta)$ with $\Upsilon(\theta)=\cos^{k}(\theta)$ by integrating Eq. (56) for the natural numbers $k=1$, 2, 3 and 4 are:

$k=1$

[TABLE]

$k=2$

[TABLE]

$k=3$

[TABLE]

$k=4$

[TABLE]

The corresponding analytical expressions of $\mathscr{T}(\theta)$ with $\Upsilon(\theta)=\sin^{k}(\theta)$ by integrating Eq. (56) for the natural numbers $k=1$ and 2, are:

$k=1$

[TABLE]

$k=2$

[TABLE]

In spite of the formal difference between Eq. (116) and Eq. (135), they are equivalent. They can be converted to each other knowing that $\sin^{2}\theta+\cos^{k}\theta=1$, making the transformation $\alpha$ into $-\alpha/(1+\alpha)$ and considering the sign of $\alpha$.

Appendix B Limit forms of $f(r,\theta)$

Integrating Eq. (60) by parts, we obtain the following expression equivalent to Eq. (58)

[TABLE]

As $\theta\to 0$ and $\theta\to\pi/2$, also $\vartheta\to 0$ and $\vartheta\to\pi/2$, respectively. Moreover, $\theta\to 0$ both $\theta$ and $\vartheta$ are infinitesimals of the same order, so that we can write $\vartheta=p\theta$ where $p$ is constant. Expanding the function $\tan x$ for small values of $x$ and knowing that $\alpha\Upsilon<1$, we derive

[TABLE]

Taking the limit of this expression for $\theta\to 0$, and identifying $r(\theta=0$ with stellar surface it comes that

[TABLE]

which introduced into Eq. (72) leads to

[TABLE]

where $\Omega_{\rm p}$ is the polar angulat velocity, which carries that $\Omega_{\rm p}/\Omega_{\rm e}=1+\alpha$ when $\Upsilon(\theta)=\cos^{k}(\theta)$ and $\Omega_{\rm p}/\Omega_{\rm e}=(1+\alpha)^{-1}$ for $\Upsilon(\theta)=\sin^{k}(\theta)$. We note that the dependence of $f^{\rm pole}$ with $\alpha$ and $k$ is given through the radii ratio $R_{\rm p}/R_{\rm e}$.

To calculate the limit of $f(r,\theta)$ as $\theta\to\pi/2$ we make the change of variables $\vartheta\to\pi/2-\epsilon_{\vartheta}$ and $\theta\to\pi/2-\epsilon_{\theta}$, where $\epsilon_{\vartheta}\to 0$ and $\epsilon_{\theta}\to 0$ are infinitesimals of the same order. This operation implies the following changes of functions

[TABLE]

and similar expressions for the functions dependent on $\vartheta$. Introducing this terms into Eq. (138), it is obtained

[TABLE]

which for $\epsilon_{\theta}\to 0$ reduces to

[TABLE]

where

[TABLE]

and $\Omega_{\rm o}/\Omega_{\rm e}=1$ when $\Upsilon(\theta)=\cos^{k}(\theta)$, or $\Omega_{\rm o}/\Omega_{\rm e}=1/(1+\alpha)$ if $\Upsilon(\theta)=\sin^{k}(\theta)$.

Finally, with Eq. (155) and Eq. (72) we obtain

[TABLE]

which is independent of $\alpha$ and $k$, and where for the stellar surface it is $\lim_{\theta\to\pi/2}r=1$ .

The forms Eq. (144) and Eq. (161) can also be derived using the analytical expressions from Eq. (111) to Eq. (135) by using series expansions for $\ln\tan{x}$, $\arctan{x}$ and $\ln(1\pm{x})$ valid for $x<1$ and relations in Eq. (151).