The understanding of the rotational evolution of early-type stars is deeply related to that of anisotropic mass and angular momentum loss. In this paper, we aim to clarify the rotational evolution of rapidly rotating early-type stars along the main sequence (MS). We have used the 2D ESTER code to compute and evolve isolated rapidly rotating early-type stellar models along the MS, with and without anisotropic mass loss. We show that stars with $Z=0.02$ and masses between $5$ and $7~M_\odot$ reach criticality during the main sequence provided their initial angular velocity is larger than 50% of the Keplerian one. More massive stars are subject to radiation-driven winds and to an associated loss of mass and angular momentum. We find that this angular momentum extraction from the outer layers can prevent massive stars from reaching critical rotation and greatly reduce the degree of criticality at the end of the MS. Our model includes the so-called bi-stability jump of the $\dot{M}-T_{\rm eff}$ relation of 1D-models. This discontinuity now shows up in the latitude variations of the mass-flux surface density, endowing rotating massive stars with either a single-wind regime (no discontinuity) or a two-wind regime (a discontinuity). In the two-winds-regime, mass loss and angular momentum loss are strongly increased at low latitudes inducing a faster slow-down of the rotation. However, predicting the rotational fate of a massive star is difficult, mainly because of the non-linearity of the phenomena involved and their strong dependence on uncertain prescriptions. Moreover, the very existence of the bi-stability jump in mass-loss rate remains to be substantiated by observations.

Gravity darkening (GD) and flattening are important consequences of stellar rotation. The precise characterization of these effects across the HRD is crucial to a deeper understanding of stellar structure and evolution. We seek to characterize such important effects on Sargas, an evolved, fast-rotating, intermediate-mass star, located in a region of the HRD where they have never been directly measured as far as we know. We use our numerical model CHARRON to analyze interferometric (VLTI/PIONIER) and spectroscopic (VLT/UVES) observations through a MCMC model-fitting procedure. The visibilities and closure phases from the PIONIER data are particularly sensitive to rotational flattening and GD. Adopting the Roche approximation, we investigate two GD models: (1) the beta-model (classical von Zeipel's law), and (2) the omega-model. Using this approach we measure several physical parameters of Sargas, namely, equatorial radius, mass, equatorial rotation velocity, mean Teff, inclination and position angle of the rotation axis, and beta. In particular, we show that the measured beta leads to a surface flux distribution equivalent to the one given by the omega-model. Thanks to our results, we also show that Sargas is most probably located in a rare and interesting region of the H-R diagram: within the Hertzsprung gap and over the hot edge of the instability strip. These results show once more the power of optical/IR long-baseline interferometry, combined with high-resolution spectroscopy, to directly measure fast-rotation effects and stellar parameters, in particular GD. As was the case for a few fast rotators previously studied by interferometry, the omega-model provides a physically more profound description of Sargas' GD, without the need of a beta exponent.

Recent interferometric data have been used to constrain the brightness distribution at the surface of nearby stars, in particular the so-called gravity darkening that makes fast rotating stars brighter at their poles than at their equator. However, good models of gravity darkening are missing for stars that posses a convective envelope. In order to better understand how rotation affects the heat transfer in stellar convective envelopes, we focus on the heat flux distribution in latitude at the outer surface of numerical models. We carry out a systematic parameter study of three-dimensional, direct numerical simulations of anelastic convection in rotating spherical shells. Restricting our investigations to hydrodynamical models with a thermal Prandtl number fixed to unity, we consider both thick and thin (solar-like) shells, and vary the stratification over three orders of magnitude. We measure the heat transfer efficiency in terms of the Nusselt number, defined as the output luminosity normalised by the conductive state luminosity. We report diverse Nusselt number profiles in latitude, ranging from brighter (usually at the onset of convection) to darker equator and uniform profiles. We find that the variations of the surface brightness are mainly controlled by the surface value of the local Rossby number: when the Coriolis force dominates the dynamics, the heat flux is weakened in the equatorial region by the zonal wind and enhanced at the poles by convective motions inside the tangent cylinder. In the presence of a strong background density stratification however, as expected in real stars, the increase of the local Rossby number in the outer layers leads to uniformisation of the surface heat flux distribution.

Stars may be understood as self-gravitating masses of a compressible fluid whose radiative cooling is compensated by nuclear reactions or gravitational contraction. The understanding of their time evolution requires the use of detailed models that account for a complex microphysics including that of opacities, equation of state and nuclear reactions. The present stellar models are essentially one-dimensional, namely spherically symmetric. However, the interpretation of recent data like the surface abundances of elements or the distribution of internal rotation have reached the limits of validity of one-dimensional models because of their very simplified representation of large-scale fluid flows. In this article, we describe the ESTER code, which is the first code able to compute in a consistent way a two-dimensional model of a fast rotating star including its large-scale flows. Compared to classical 1D stellar evolution codes, many numerical innovations have been introduced to deal with this complex problem. First, the spectral discretization based on spherical harmonics and Chebyshev polynomials is used to represent the 2D axisymmetric fields. A nonlinear mapping maps the spheroidal star and allows a smooth spectral representation of the fields. The properties of Picard and Newton iterations for solving the nonlinear partial differential equations of the problem are discussed. It turns out that the Picard scheme is efficient on the computation of the simple polytropic stars, but Newton algorithm is unsurpassed when stellar models include complex microphysics. Finally, we discuss the numerical efficiency of our solver of Newton iterations. This linear solver combines the iterative Conjugate Gradient Squared algorithm together with an LU-factorization serving as a preconditionner of the Jacobian matrix.