Photometric determination of rotation axis inclination, rotation rate, and mass of rapidly rotating intermediate-mass stars

Intermediate-mass stars are often fast rotators, and hence are centrifugally flattened and affected by gravity darkening. To analyse this kind of stars properly, one must turn to 2D models to compute the visible radiative flux and to take the geometrical effect of the star inclination into account. Assuming a given stellar age and chemical composition, we aim to derive the mass and rotation rates of main sequence fast rotating stars, along with their inclination, from photometric quantities. We chose three observables that vary with mass, rotation, and inclination: the infrared flux method temperature T_IRFM, the Str\"omgren c1 index, and a second index c2 built in the same way, but sensitive to the UV side of the Balmer jump. These observables are computed from synthetic spectra produced with the PHOENIX code and rely on a 2D stellar structure from the ESTER code. These quantities are computed for a grid of models in the range 2 to 7~M_Sun, and rotation rates from 30% to 80% of the critical rate. Then, for any triplet (T_IRFM, c1, c2), we try to retrieve the mass, rotation rate, and inclination using a Levenberg-Marquardt scheme, after a selection step to find the most suitable starting models. Hare-and-hound tests showed that our algorithm can recover the mass, rotation rate, and inclination with a good accuracy. The difference between input and retrieved parameters is negligible for models lying on the grid and is less than a few percent otherwise. An application to the real case of Vega showed that the u filter is located in a spectral region where the modelled and observed spectra are discrepant, and led us to define a new filter. Using this new filter and subsequent index, the Vega parameters are also retrieved with satisfactory accuracy. This work opens the possibility to determine the fundamental parameters of rapidly rotating early-type stars from photometric space observations.


Introduction
The influence of rotation on the evolution of stars is a longstanding question for their modelling (Maeder 2009).The major consequence of rotation is the so-called rotational mixing.It affects radiative regions thanks to the forcing of baroclinic flows and their associated instabilities, which generate a small-scale turbulence (e.g., Mathis et al. 2018).As a result, the lifetime on the main sequence is enhanced and surface abundances are modified (Maeder et al. 2014).In addition, the interpretation of observations is also affected, in particular for fast rotating earlytype stars, which is our focus in this article.
According to Zorec & Royer (2012), fast rotation (equatorial velocity higher than 100 km s −1 ) involves 50% of early-type stars.For such stars, 1D (spherically symmetric) models may fail to correctly interpret the observations.To illustrate this point, the example of Altair is emblematic.Altair is an A7V-type star in the solar neighbourhood, 5.14 pc away from the Sun, and whose V sin i is around 240 km s −1 .Because of its proximity, Altair has been observed in many ways; it has been a privileged target for interferometry (van Belle et al. 2001; Domiciano de Souza et al. The whole set of synthetic spectra and their associated data are only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/vizbin/cat/J/A+A/676/A502005; Peterson et al. 2006;Monnier et al. 2007;Bouchaud et al. 2020;Spalding et al. 2022), seismology (Buzasi et al. 2005;Suárez et al. 2005;Le Dizès et al. 2021), and spectroscopy (Reiners & Royer 2004).A long-standing question about Altair is its evolutionary stage, in short its age.Estimates using 1D models were rather inconclusive: Suárez et al. (2005) gave a range of 225-775 Myr, while Domiciano de Souza et al. (2005) gave another range between 1.2 and 1.4 Gyr.However, Peterson et al. (2006) interpreted their interferometric data with a simple Roche 2D model, and deduced from the measured polar radius and luminosity that Altair was barely off the zero-age main sequence (ZAMS).Some years later, with the advent of self-consistent 2D models 1 (Espinosa Lara & Rieutord 2013;Rieutord et al. 2016), Bouchaud et al. (2020) reanalysed available data of Altair from interferometry, spectroscopy, and seismology, and derived a 2D concordance model that points towards 100 Myr for Altair's age, thus indeed barely off the ZAMS.
The foregoing example of Altair is rather extreme since its equatorial velocity exceeds 300 km s −1 (Bouchaud et al. 2020).However, it emphasises the fact that 2D models are necessary to analyse fast rotators.Altair's example interestingly shows that interferometry is well suited to determine the inclination of the rotation axis on the line of sight.This angle, i, is a crucial parameter since it conditions the determination of the true luminosity, the polar and equatorial radii, and the effective temperature distribution over latitude.Fast rotating stars are singled out by gravity darkening that makes their polar regions brighter than their equatorial regions (e.g., Rieutord 2016).Hence, inclination strongly influences the way we see a fast rotating star, and therefore how we interpret its observations (e.g., line profiles, spectral energy distribution).For fast rotating early-type stars, the best way to determine i seems to be interferometry, but the angular diameter should not be too small.Presently, 1.5 milliarcsec (mas) seems to be the minimum size, but in the near future with the forthcoming Space Infrared Telescope for Cosmology and Astrophysics (SPICA) to be implemented on the Center for High Angular Resolution Astronomy (CHARA) array (Mourard et al. 2018), this limit may be reduced to 0.8 mas.Nevertheless, this new limit is still stringent on the number of early-type stars that can be analysed since they should typically be closer than 30 pc from the Sun.If we wish to go farther, we can only rely on spectroscopy and/or photometry, and try to get the relevant information through asteroseismic and line profile analyses, for instance.
Since we are interested in fast rotation, both spectroscopy and seismology face difficulties.The eigenspectrum of a fast rotator is very complex and requires 2D calculations of eigenmodes over a 2D model (e.g., Reese et al. 2021).Until now, determination of inclinations solely from asteroseismic data has been attempted only for late-type stars where the amplitude of stochastically excited waves can reasonably be guessed (Gizon & Solanki 2003;Gehan et al. 2021).For early-type fast rotators where mode excitation usually comes from the κ-mechanism, no such determination has been achieved.However, if eigenmodes can be identified, then the rotation period can be deduced.With this quantity and the projected equatorial velocity V sin i, obtained with a spectroscopic analysis, inclination can be estimated as in Moravveji et al. (2016), among others, for the slowly pulsating B-type star (SPB) star KIC7760680.Actually, since the rotation of this star is not too fast (72 km s −1 ), 1D models are still suitable to perform its seismic analysis.
Very few attempts exist from a purely photometric perspective, while several spectroscopic studies have been done.For instance, Frémat et al. (2005) considered rapidly rotating B-type stars and showed that neglecting gravitational darkening produces systematic errors on surface parameters, such as abundance, effective temperature T eff , and V sin i, emphasising the result of Townsend et al. (2004) on the underestimation of V sin i in Be stars.More specifically, Reiners & Royer (2004) attempted to determine Altair's inclination from a line profile analysis and found i > 68 • on a 1σ level and above 45 • on a 2σ level.Interferometry provides a value of 50.7 • ±1.2 • (Bouchaud et al. 2020) or 57.2 • ± 1.9 • (Monnier et al. 2007).The spectroscopic exercise was applied to Vega by Takeda et al. (2008), who find a value of 7.2 • , actually quite close to the value provided by the concordance model of Monnier et al. (2012), namely i = 6.2 • ± 0.4 • .
In the present work we propose another route to the determination of fundamental parameters of a fast rotating star.The novelty is to use the recent 2D-ESTER models in combination with the PHOENIX atmosphere models to derive a grid of spectral quantities that can be inverted to retrieve the actual physical parameters of the star.We still restrict ourselves to early-type stars of a certain age, and focus on two fundamental parameters of a fast rotating star, namely its mass and flattening, along with the inclination of the rotation axis on the line of sight.For this we compute synthetic spectra and determine three spectral quantities: T IRFM , the temperature determined by the infrared flux method (IRFM; Blackwell et al. 1980); c 1 , the Strömgren photometric index; and c 2 , a new photometric index using the UV part of the spectrum.These three quantities characterise the shape of the spectrum and we show here that it is possible to use them to determine the stellar mass, the centrifugal flattening, and the inclination on the line of sight for intermediate-mass stars with fast enough rotation.
In Sect. 2 we briefly present the numerical tools we used to generate spectra of fast rotating stellar models, then in Sect. 3 we describe the properties of the observable quantities we used.Our inversion algorithm is detailed in Sect.4, and we show how it performs using a hare-and-hounds exercise.Section 5 shows the application of our method to the concrete case of Vega and how we overcome the problems that arise when dealing with this real star.In Sect.6 we discuss questions raised by our method and by the next steps towards a full determination of fundamental parameters of early-type fast rotating stars.Our conclusions follow in Sect.7.

Computation of the observed flux of fast rotating stars
For the kind of stars we are interested in, the geometry departs enough from the spherical case so that a single pair of (T eff , log g) does not make sense anymore to characterise the atmosphere.
The emergent flux as seen by an observer consists of the sum of the contributions of all the visible parts of the star, each with its own surface effective temperature and gravity.Therefore, for the very same star, the observed flux depends on the inclination of the rotation axis with respect to the line of sight.We thus need to calculate the spectrum as seen from an observer before computing the fluxes in various photometric bands in order to infer photometric indices.For this we use the 2D structure of the star, which provides the local effective temperature and gravity of each surface element of the visible part.These two quantities are used as inputs for the computation of a model atmosphere, which in turn is used for the synthesis of the spectrum emerging from this surface element.

Two-dimensional internal structure code ESTER
The 2D internal structure is computed with Evolution STEllaire en Rotation (ESTER), a stellar evolution code dedicated to the modelling of rapidly rotating stars (Rieutord & Espinosa Lara 2009, 2013a,b;Espinosa Lara & Rieutord 2013;Rieutord et al. 2016).ESTER relies on the following assumptions: a steady state structure that includes the baroclinic flows in radiative regions (differential rotation and meridional circulation), and hydrogen burning in the convective core.Evolution is presently only possible along the main sequence.The equations are solved in a non-spherical geometry that follows the centrifugal distortion of the star.Solutions are obtained using spectral elements and a discretisation of the star in 'onion' shells, where the radial and the horizontal directions are respectively discretised using a Gauss-Lobatto and a Gauss-Legendre collocation grid.We thus have access to the latitudinal variations in all the fundamental parameters such as the effective temperature T eff (θ), the surface ω n θ ≤ 0.3 4 0.3 < ω ≤ 0.7 8 0.7 < ω ≤ 0.9 10 >0.9 ≥12 gravity g(θ), and the radius R(θ) thanks to the Legendre polynomial interpolation, where θ is the colatitude.Thanks to the spectral discretisation of the θ grid, accurate models can be computed with a small number of grid points.The spectral convergence of ESTER models in the angular dependence of T eff is shown in Fig. 1, where ω is the ratio of the actual equatorial rotation rate to the Keplerian rate at the equator.This figure shows that we need fewer than 22 spherical harmonics (or 11 grid points2 ) to ensure a truncation error less than 10 −3 at all rotation rates less than ω = 0.9.In Table 1 we summarise the latitudinal resolution we used to minimise both spectral error and computing time.

Spectrum synthesis
Synthesising a spectrum requires a model atmosphere.Unlike the 1D case where effective temperature and gravity are constant over the stellar surface, we have to compute several atmosphere models following the local surface flux and gravity provided by the internal structure at each point of the visible surface of the star.Each atmosphere model is used to calculate a contribution to the observed spectrum.This approach is used in Bouchaud et al. (2020) with the PHOENIX spectra from the Göttingen database (Husser et al. 2013).However, using precomputed spectra requires applying a weighting to account for the limb darkening effect (Claret 2000).In the present study we use the PHOENIX code (Hauschildt et al. 1999) to compute all the model atmospheres and corresponding spectra from scratch, which provides us with the intensity spectra for different directions, the limb darkening effect being included inevitably.To retrieve the specific intensity for any direction, we used a Legendre polynomial defined on a four-point Gauss-quadrature grid in µ = cos γ, where γ is the angle between the normal to the surface and the direction of interest.We checked that such a small Gauss grid was able to ensure a numerical precision better than 5% everywhere (actually better than 1% in most cases).
The model atmospheres and the spectra are calculated with the conservative microturbulence value of 2 km s −1 and the same set of abundances as the stellar structure, namely the chemical composition of Asplund et al. (2009) with the meteoritic values of Lodders & Palme (2009) for the refractory elements, as suggested by Serenelli (2010).Convection is treated with the mixing-length theory using a ratio of 1.25 for mixing length to pressure scale height.Since the relative radiative flux difference between a plane-parallel and a spherical model is below 0.1 % in the wavelength interval we use (see Fig. 2), we performed all the calculations with the plane-parallel geometry to reduce the computing time without impairing the results.Local thermodynamic equilibrium (LTE) is assumed for the model atmosphere and the spectrum synthesis.

Geometric grid of the visible stellar surface
The method we adopted to discretise the stellar surface is close to that described in Bouchaud et al. (2020) and ensures that all the surface elements share roughly the same area.The star is first sliced along a meridian into N θ latitudinal strips of constant angular width ∆θ = π/2N θ , where N θ is the number of atmosphere strips and n θ is used to converge the internal structure, and are independent.Keeping the size of each surface element roughly constant implies that the number of longitudinal elements varies from one latitudinal strip to the other.The typical size ∆S 0 of a surface element is set in the strip surrounding the pole (of mean colatitude θ 0 ) so as to have N φ 0 = N φ (θ 0 ) surface elements.The parameter N φ 0 is chosen to have a spectrum as smooth as possible and we show below how its value is determined.To have the Doppler shifts symmetric compared to the rest frame, the surface elements are located symmetrically with respect to the meridian that contains the line of sight.Then, for any other strip of mean colatitude θ k , the number of longitudinal A50, page 3 of 12 points N φ (θ k ) is set to have ∆S (θ k ) as close as possible to ∆S 0 .The quantity ∆S (θ k ) is defined as where r is the distance to the star centre, r θ = ∂r/∂θ, and ∆φ = 2π/N φ (θ k ) is the angular step in longitude.The visible surface is then obtained by selecting elements with where i is the angle between the line of sight and the rotation axis.

Spectrum in the observer's frame
At this stage, for each visible surface element located at (θ k , φ j ), the monochromatic contribution ∆I λ (θ k , φ j ) to the observed spectrum at wavelength λ is defined as where is the specific intensity along the line of sight of a model atmosphere computed with (T eff (θ k ), g(θ k )) and ∆S proj (θ k ) is the surface element ∆S (θ k ) located at φ j and projected on the sky.The observed spectrum is then the sum of each of these contributions, Doppler-shifted according to the motion along the line of sight of the surface element in question.
Figure 3 shows a close-up of the resulting synthetic spectrum for a 5 M star with ω = 0.6 around the MgII λ4481 line, computed for three different inclinations of the rotation axis, various sets of (N θ , N φ 0 ), and a 1 pm wavelength step.The spectra are normalised to the continuum level, which was computed in the same way, but without the line opacities.The wavy aspect of the spectra, especially visible at high inclinations, is due to the discretisation of the stellar surface.This numerical noise becomes negligible by setting (N θ , N φ 0 ) to (100,20) or higher, as shown in Fig. 3.

Selected photometric quantities
The flux received by an observer depends on the effective temperature and effective gravity latitudinal distributions together with the inclination of the rotation axis on the line of sight.However, Espinosa Lara & Rieutord (2011) showed that the flux and the effective gravity are tightly related in the radiative envelope of early-type stars.Thus, the flux and the effective gravity distributions are both governed by the shape of the star.Hence, if we consider main sequence stars at a given evolutionary stage, we may expect that the energy flux towards an observer basically depends on two fundamental parameters of the star, its mass M and its centrifugal flattening ε, and on the angle i between its rotation axis and the line of sight.Instead of the centrifugal flattening ε, we use ω, the ratio of the equatorial angular velocity to the Keplerian angular velocity at equator.In the Roche model ω and ε are simply related by ε = ω 2 /(2 + ω 2 ), which is very close to the relation verified by the more realistic ESTER models.
These three quantities are expected to influence the shape of the star's spectrum, even more so when the star rotates rapidly.We therefore investigate the possibility that two colour indices together with a temperature expressing the surface heat flux may be able to constrain i, ω, and M, if the age is given.
To this end, we chose to focus on the ZAMS, with solar metallicity.We built a grid of models with masses from 2 to 7 M and rotational rates ω between 0.3 and 0.8.The rotation axis inclinations ranged from 0 to 90 • with a 5 • step.The spectra were computed from 10 nm to 9000 nm with a 0.1 nm step.

The infrared flux method
Since we need an effective temperature representative of the surface heat flux of the star, we resort to that derived with the wellknown IRFM, which was first devised by Blackwell & Shallis (1977) to yield angular diameters along with effective temperatures.As a follow-up, Blackwell et al. (1980) introduced the ratio R = F/F λ IR of the bolometric flux F to the infrared monochromatic flux F λ IR to find effective temperatures.They computed theoretical values of this ratio R th for a set of T eff and several wavelengths λ IR .Then, for an observed value of R, the IRFM temperature T IRFM was retrieved by interpolation in the table linking R th and T eff .
We computed values of R th with PHOENIX atmosphere models with 8000 K ≤ T eff ≤ 25 000 K and log g = 4.0, which are typical values for the kind of stars we study.The bolometric flux is obtained by an integration between the above-mentioned boundaries, the relative difference with larger intervals being less than 10 −5 .The results are shown in Fig. 4 along with the original values of Blackwell et al. (1980) for log g = 4.For similar temperatures and given λ IR , our values of R th and those of Blackwell et al. are almost superimposed, showing that the relation between T eff and R th depends very weakly on the model atmosphere and the set of abundances, as already stated by Blackwell et al. (1980).
We then investigate the sensitivity of T IRFM to measurement errors of the flux.For that we introduce a 1% variation to the bolometric flux and monitor the resulting variations in T IRFM on our grid of ESTER/Phoenix models with various rotation rates.As shown in Fig. 5, the induced variation in T IRFM in absolute value is 70 K at most.A similar 1% variation introduced to the A50, page 4 of 12 Fig. 5. Variation in T IRFM induced by a 1% variation in the bolometric flux for our grid of ESTER/Phoenix models.For each mass, several rotation rates ω were considered, ranging from 0.3 to 0.7 with a 0.1 step.flux at a selected IR wavelength has almost no effect on T IRFM .Hence, the error on T IRFM essentially comes from that on the bolometric flux.
We next investigate the sensitivity of T IRFM to the geometry of the star, which is determined by the pair (ω, i).To this end, we consider the ratio namely the IRFM temperature scaled by that of the associated non-rotating model of the same mass.Quite nicely, the ratio Q is almost independent of the mass if M ≥ 3 M , as shown by Fig. 6.Hence, the IRFM temperature of a centrifugally distorted star viewed under an inclination angle i is that of the corresponding non-rotating star corrected by a 'shape factor' Q 1 actually describing the impact of gravity darkening on the flux received by the observer.The subsequent approximate relation main sequence evolutionary tracks computed with non-rotating 1D models.This independence of the Q-factor with respect to mass seems to disappear below 3 M (e.g., 2 M models, red crosses in Fig. 6).The dependence is not strong but notable, especially at high angular velocity and high inclination, for an as-yetunknown reason.

The Strömgren index c 1
The next spectral quantity that we consider is the Strömgren index c 1 .Strömgren (1963)  We use the Kurucz (1993) program to compute synthetic c 1 indices from our spectra.To allow direct comparisons with observed photometry, this program takes into account the transmission of the filters, the transmission of the Earth's atmosphere, and the response of the detector.The zero point is set by equating  the c 1 value, which we obtained with the spectrum of our Vega model (see Table 2), to the observed value of Crawford & Barnes (1970).
The variations in c 1 with T eff and log g of spherical models are illustrated in Figs. 2 and 3 of Moon & Dworetsky (1985).For stars hotter than about 9500 K, at given log g, c 1 is an indicator of T eff , which increases with decreasing c 1 .A small dependence versus log g also exists.In the case of a fast rotating star, log g and T eff decrease from pole to equator and both variations increase c 1 .The overall effect with respect to the rotation axis inclination for such stars is plotted in the top panel of Fig. 8 for a 5 M model: the lowest value of c 1 is met when the star is seen pole on.
For temperatures less than 8500 K, this effect reverses: c 1 decreases from pole to equator, as illustrated in the bottom panel of Fig. 8 with a 2 M model.In this case the temperature effect overcomes the log g effect.As can be guessed, between these two temperatures, the dependence of c 1 with inclination is weaker and can even be non-monotonic.

Hubble Space Telescope colour index c 2
To probe the shape of our spectra at wavelengths shorter than the Balmer jump, we built a new index c 2 in a similar way to c 1 .The c 1 index can be interpreted as a measure of the curvature of the spectrum at the low-frequency end of the Balmer jump.Thus, we built c 2 to measure the curvature of the spectrum at the high-frequency end of the Balmer jump.Thus, we need two filters on the blue side of the Balmer jump.For this new index we three filters from the Hubble Space Telescope (HST) filter set, namely HSP_VIS_F419N_B, ACS_HRC_F250W, and HSP_UV1_F220W_A, respectively referred to as F419, F250, and F220 in the following.Other choices are possible as long as the curvature of the spectrum is monitored.Thus, we define c 2 as The transmissions of these three filters are shown in Fig. 9.In Fig. 10 we show the variations in c 2 with inclination and rotation rate for a 2 M and a 5 M model.These variations are actually opposite to those of c 1 , which is a good property for inversion matters.

M, ω, and i determination from photometric observables
In the previous section we computed the three functions on a grid of triplets (M, ω, i) for ZAMS models.Since the computation of the three observables is indeed computationally very demanding, if we need to evaluate these three functions at any mass, rotation or inclination, we are forced to resort to interpolations by using tri-cubic splines.The drawback is that these interpolations introduce some numerical noise, which perturbs somewhat the process to retrieve (M, ω, i) from (T IRFM , c 1 , c 2 ) as we show below.

Initialisation
As for any resolution of a non-linear problem, we need an initial guess not too far from the target solution to start the solver.Hence, we first need a reasonable guess of M.
Since our models have a fixed metallicity and age, their effective temperature, and thus T IRFM , is essentially controlled by their mass, as illustrated by Fig. 11.Clearly (ω, i) impacts T IRFM , and the relation between T IRFM and M can only be approximate.We derived the following fit which is shown in Fig. 11 as a black line.This relation provides an approximate mass with a maximum error of 15% for the worst case (i.e., when M = 2 M , ω = 0.3, i = 90 • ).Hence, this relationship determines a set of possible values of M, while ω and sin i are still unconstrained.
To limit the range of potential solutions, we first consider two series of models: pole-on models (i = 0 • ) and equator-on models (i = 90 • ).For these two inclinations, we determine the range of ω for which the three observed values of (T IRFM , c 1 , c 2 ) are between the extreme values imposed at i = 0 • and i = 90 • : Here the superscript 'obs' stands for the observed (i.e., target) values.The result of this selection is illustrated in Fig. 12, where we show the values of ω) that meet the above requirements, for an observed ZAMS model of 5 M at ω = 0.5 and i = 65 • .8), colour-coded by the scattering between the inclinations obtained from each observable, as defined by the left-hand side of Eq. ( 10).Green crosses stand for initialisation models for which Eq. ( 10) is verified with ε = 0.1.The big red cross shows the location of our observed star.Models complying with system (8) are still numerous, as shown in Fig. 12.To further restrict the possibilities, we consider each possible solution (M k , ω k ) of this set (i.e., each dot in Fig. 12) and solve for the inclination for each observable.In other words, we solve independently the three equations with the Newton-Raphson algorithm, where f T IRFM , f c 1 , and f c 2 are the tri-cubic splines.We thus obtain three values for the star's inclination for each possible pair (M k , ω k ).Obviously, a suitable pair is one for which the three inclinations are equal.Due to the analytic fits, we cannot expect a strict equality and we select the acceptable pairs (M k , ω k ) by imposing that there is an inclination i 0 such that where i(T IRFM ) is the inclination which meets the T IRFM constraint (9a) for (M k , ω k ), and mutatis mutandis for i(c 1 ) and i(c 2 ).
If such an inclination i 0 exists for a small enough value of ε, then a model matching the observation exists as well.Numerical experiments showed that ε should not be larger than 10 −1 for the algorithm to converge to the target solution.This initial guess (M k , ω k , i 0 ) is then used to derive the final solution with a Levenberg-Marquardt algorithm.We show in Fig. 12 an example of the set of models that match Eq. ( 10) (green crosses).Some difficulties arise, however, and come from the fact that the three functions of Eq. ( 6) are only known thanks to interpolating splines, which reproduce the values of our grid of models with limited accuracy out of the grid points.Hence, some differences can appear between the 'true' values of T IRFM , c 1 , and c 2 and those computed using the fits for model that are not on the grid points, so that we may not be able to retrieve the exact solution (M, ω, i) in these cases.(top), ω (middle), and i (bottom) for models with (M, ω) = (4.8,0.42), (5.5, 0.62), (6.6, 0.72), (6.6, 0.32), and (3.4,0.72; blue, orange, green, red, and purple lines resp.), with M in M .

Robustness of our resolution
Since the interpolating splines cannot give us an exact inversion, it is worth evaluating the error generated by the inversion To this end, we apply the inversion to all the observed quantities sampled on our gridpoints.In most cases the algorithm succeeds in retrieving the parameters of the model, but it fails sometimes, especially in the neighbourhood of the grid boundaries.In Fig. 13 we show maps of errors in the (ω, i) plane for three masses of the grid.In most cases the parameters of the models are obtained with a negligible error.
The next test is to invert observables associated with models that are not on the gridpoints.In Fig. 14 we show the error measured on five models that are representative of the general trend.Rieutord (2013) that match the observational constraints (interferometric and spectroscopic) derived by Monnier et al. (2012).
Typically, the error strongly depends on the mass range: it is a few percent when M ≥ 3.5 M , but may exceed 10% for lower masses.This sensitivity to interpolation errors is amplified by the weak variations in indices c 1 and c 2 for these masses.Inclination errors may reach ±15 • in the worst cases, namely for rotators viewed pole-on.Even though these errors may seem important, they can easily be reduced with the use of a finer grid, and focusing on a narrower range of parameters in the search process.

Preliminaries
Vega (α Lyrae) is a traditional standard of photometry (Bohlin et al. 2019), and as such has been observed in many ways, especially with interferometers in recent years (Aufdenberg et al. 2006;Peterson et al. 2006;Monnier et al. 2012).Using observations with the CHARA-array and MIRC beam combiner operating at 1.65 µm, Monnier et al. (2012) derived a concordance model that matches interferometric and spectroscopic constraints.Since the aim of our work is to derive fundamental parameters from spectroscopic indices, we naturally test our method on Vega.
As a first step, we derived a 2D-ESTER model that matches the most robust observational constraints.From the concordance model of Monnier et al. (2012), we assume an inclination of Vega rotation axis of i = 6.2 • , which implies an equatorial velocity of 196 km s −1 using the spectroscopically derived value of V sin i = 21.2 km s −1 (Takeda et al. 2008).We also adopt the polar and equatorial radii given by Monnier et al. (2012) and the derived polar and equatorial effective temperatures.We match these observational data with a 2D-ESTER model of mass M = 2.374 M with metallicity Z = 0.0093 (Yoon et al. 2008), a hydrogen mass fraction of X env = 0.7546 in the envelope, and X core = 0.2045 in the convective core, as actually found by Espinosa Lara & Rieutord (2013).Hereafter, we use instead x c = X core /X env = 0.271 to denote the evolutionary stage of the model.All the parameters of our model are summarised in Table 2.We note that our mass is higher than that of the concordance model of Monnier et al. (2012).It may be a consequence of the Roche approximation used by these authors.This difference is however not critical for our test.Comparison between the spectrum from CALSPEC and from our model for Vega.Top panel: spectra of Vega from CALSPEC and from our model (orange and blue lines, respectively; left scale) and energy transmission curves of the filters used to compute the colour indices (right scale).The Strömgren filters are shown as dashed lines, while the HST filters as solid lines, the colours having the same meaning as in Figs.7 and 9, respectively, except that the u filter is replaced here by u .Bottom panel: difference between the flux of our Vega model and that of CALSPEC (black solid line).
The next step is to compute the synthetic spectrum for this model using an inclination of 6.2 • .We used the abundances derived from Royer et al. (2014), which comprises a large sample of chemical elements.The abundance differences with studies that consider gravity darkening (e.g., Yoon et al. 2010;Takeda et al. 2018) are not significant, except for carbon.Nevertheless, the global agreement between our synthetic spectrum and the observed spectrum is quite good.We show in Fig. 15 the spectral energy distribution (SED) of our model together with its difference with the calibrated CALSPEC spectrum derived from observations (e.g., Bohlin et al. 2014;Bohlin & Lockwood 2022).From this synthetic spectrum we derived the three spectral observables (T IRFM , c 1 , c 2 ), which we display in Table 3 along with their values derived from the CALSPEC spectrum.We note that the observed and model values of T IRFM and c 2 agree within a 1% error bar, while the c 1 values show a mismatch of around 6%.We trace the origin of this discrepancy to the flux difference that appears in the spectral band of the Strömgren u-filter, namely for wavelengths between 300 and 360 nm (see Fig. 15).Actually, this difference between the observed and synthetic spectrum has already been noted (e.g., Castelli & Kurucz 1994;Bohlin & Gilliland 2004), and is explained by the difficulty faced by stellar atmosphere models to deal with the complex atomic physics needed to compute this spectral region.This mismatch between observed and synthetic c 1 is a serious obstacle for our inversion procedure since forcing our model to match the observed c 1 would necessarily lead to incorrect stellar model parameters.
To circumvent this difficulty we chose to replace the Strömgren u-filter with a more appropriate one that does not sample this difficult region of the spectrum.Namely, we consider a new filter u , which has the same shape as the u-filter, but is centred at a shorter wavelength λ c .
To determine the best wavelength for this, we computed the index c 1 = (u − v) − (v − b) scanning the spectrum from λ c = 170 nm to 350 nm and computing the difference between c 1 of our Vega model and that of the CALSPEC spectrum.The result is visualised in Fig. 16, where we see that the best match occurs at λ c = 215 nm.The corresponding c 1 values are reported in the last column of Table 3.The difference between CALSPEC and our model is now much less than a percent.

Inversion of Vega data
The challenge is now to use the three observed indices T IRFM , c 1 , and c 2 in Table 3 obtained with the CALSPEC spectrum to retrieve the fundamental parameters M and ω of Vega along with the inclination of its rotation axis on the line of sight.In the present case we do not start from scratch and assume that the evolutionary state of Vega is known.Specifically, we consider a set of ESTER models computed with the metallicity and core A50, page 9 of 12 hydrogen mass fraction of Vega.The model atmospheres and spectra were calculated accordingly.
From the value of T IRFM we deduce that the mass should be between 2 M and 3 M .We thus scan this mass range to derive the analogue of relation ( 7), but with x c = 0.271.From the observed CALSPEC value of T IRFM we find that the approximate mass of Vega is 2.46 M with a typical uncertainty of 15%.This suggests scanning a mass range of [2.1, 2.8] M .Since we have no idea of the rotation rate, except the indication from the V sin i that the equatorial velocity is higher than 21.2 km s −1 , we decided to scan the whole interval of ω ∈ [0.3, 0.8] and to inspect the resulting inclinations.The subsequent set of initial guesses is illustrated in Fig. 17.This diagram already shows that Vega must be a rapid rotator since inclination-consistent solutions (green crosses) point to ω larger than 0.5.Owing to the low value of V sin i we can already deduce that Vega is seen nearly pole-on.
As illustrated in Fig. 18, the selection steps resulted in 42 initial guesses, which converged towards two solutions: 4 to a model with M = 2.374 M , ω = 0.511, and i = 0 • , the others to a solution with M = 2.467 M , ω = 0.663, and i = 32.1 • , all with negligible uncertainties.The existence of these two solutions is due to some degeneracy in the observables.We traced back this degeneracy to the limited precision of our calculations caused by the interpolations between gridpoints.However, the accuracy of the observable determination is not infinite, which is also a source of degeneracy of the solutions.However, the solution with M = 2.467 M yields V sin i = 137 km s −1 , which is easily ruled out by the observed value of V sin i = 21.2 km s −1 (Takeda et al. 2008).If we take the M = 2.374 M model and combine it with the observed V sin i, we get i = 6 • as expected.
The example of Vega shows that our method can reveal the high rotation rate of a star, but may lead to trouble by degeneracy in some parameter range.

Discussion
These results have raised several questions.Within the assumptions we used, namely the fixed stage of evolution of our models, we can wonder about the impact of uncertainties that naturally come with measurements.For instance, if the three observables are derived with a one percent precision, we expect a precision of a few percent, but polar and equatorial orientations tend to increase the influence of uncertainties.A slow rotation (ω 0.3) producing a weak gravity darkening also amplifies errors on modelling, and this is why we did not consider such cases.Lastly, the example of Vega has shown that some degeneracy may prevent us from determining the parameters unambiguously.In this case our method shows that Vega is a rapid rotator with ω ≥ 0.5 and offers two possible solutions, but we cannot decide between the two.We have to resort to the value of V sin i or to the rotation period to lift the degeneracy.Scanning our grid of models and systematically testing inversions showed that a degeneracy like the previous one especially affects the mass range [2, 3.5] M .The use of other spectral quantities or additional photometric indices may help circumvent such difficulties.This option still needs to be fully explored.
The next question is obviously the determination of the evolutionary stage of the star from the spectroscopic data.From the location in the Hertzsprung-Russell diagram, giants can be separated from main sequence stars.Among the latter, as a star expands during its main sequence evolution, the radius would be the main indicator.To evaluate it, a measurement of the surface brightness and the ensuing angular diameter (Di Benedetto 2005) together with a reliable parallax may help determine the relevant starting models for our inversion method.A50, page 10 of 12 According to Fig. 6, the Q factor varies little with mass at a given age.We wanted to know how this factor may change with age or metallicity.For this we computed Q factors for models with solar metallicity at the evolutionary stage of Vega (i.e., x c = 0.271) on the one hand, and models with Vega's composition on the other (ZAMS and x c = 0.271).At ZAMS, the ratio of Q between solar and Vega metallicity models is very close to unity with a dispersion close to 10 −3 .Hence, Q depends very little on metallicity.The same seems to be true with respect to the stage of evolution when we compared the Q values of ZAMS models with models at the stage of evolution of Vega, both computed with Vega's abundances.Here too, we note that the ratio is around unity and the scatter is a few thousandths.This latter result is interesting as it allows us to use tabulated values of Q to simulate the effect of rotation on isochrones of stellar clusters.Comparisons between simulations and data would lead to new constraints on cluster ages.

Conclusion
In the present work we studied the effect of gravity darkening on three spectral quantities and examined the possibility to retrieve two fundamental parameters of a fast rotating star and the inclination of its rotation axis on the line of sight.To test this idea, we used the 2D-ESTER model of Espinosa Lara & Rieutord (2013) combined with the stellar atmosphere model PHOENIX (Hauschildt et al. 1999) to compute synthetic spectra for a grid of models representing rapidly rotating stars of intermediate mass.Hence, our grid scans the mass range 2 M -7 M and rotation rates from 30% to 80% of the critical angular velocity.Synthetic spectra were computed for any inclination of the rotation axis on the line of sight.
With this grid of models and the associated grid of spectra, we computed the temperature derived from the infrared flux method T IRFM , the Strömgren index c 1 = (u − v) − (v − b), and a new index c 2 using HST filters in the UV domain.The indices c 1 and c 2 are sensitive to the shape of the spectrum near the Balmer jump, while T IRFM reflects the global shape of the spectrum, all three quantities being sensitive to the gravity darkening.Our goal is to derive the mass M, the rotation rate ω, and the inclination i of the rotation axis from the triplet (T IRFM , c 1 , c 2 ), restricting our computations to ZAMS models of solar metallicity in this preliminary exploration.Since our grid is rather coarse and the dependence of (M, ω, i) versus (T IRFM , c 1 , c 2 ) is non-linear, we have to resort to interpolations to solve the inversion.This introduces some numerical noise that limits the precision of the inversion, and in a few cases prevents the computation of the right solution.This drawback can be circumvented somewhat by a refinement of the grids, but at higher computational cost.
We first tested the method with a hare-and-hound exercise using our grid of the ZAMS model.As the results were satisfactory (Figs. 13 and 14), we moved to a test with real data, considering the observed spectra of Vega from the CALSPEC database 3 .After designing a 2D model of Vega that meets the data from interferometric and spectroscopic observations, we computed the synthetic spectrum of this star and derived the three observables (T IRFM , c 1 , c 2 ).We also computed these three values using the calibrated spectrum from CALSPEC and found that c 1 differed from the synthetic value by 6%, while T IRFM and c 2 agreed within a percent.The discrepancy between the observed c 1 and the synthetic index has been traced back to the known discrepancy between the synthetic spectrum and the observed one at wavelengths just below the Balmer jump (e.g., Hauschildt et al. 1999;Bohlin & Gilliland 2004).Since such a mismatch prevents our algorithm from working properly, we replaced the Strömgren u-filter, centred at 350 nm, by another filter of the same profile, but centred at 215 nm.At this wavelength the new index c 1 has the same value in both the observed and the synthetic data.To test our inversion method we computed a grid of 2D-ESTER models with masses ranging from 2 to 3 solar masses and rotation rates from 30% to 80% of the critical rate, at an age (defined by the core hydrogen content) coming from other studies.Browsing within this grid we could immediately conclude that Vega is a fast rotator.However, we encountered a difficulty since our algorithm produced two solutions, revealing a degeneracy in the mass range of Vega, namely between ∼2 M and ∼3.5 M .In the case of Vega, this degeneracy is immediately lifted with the use of the V sin i value, which is compatible with only one of the models.
These first results are quite encouraging, but our method can still be improved.An obvious way is to reduce the interpolation errors by a finer grid in mass and rotation rates, but additional photometric indices may also be helpful.Moreover, we think that the method can be extended to age determination within the main sequence for stars of intermediate mass with a fast rotation, provided that additional information sensitive to the stellar radius is included.
The limits of our approach can be summarised as follows: 1.The gravity darkening should be strong enough, and it seems that this corresponds to a rotation rate above 30% of the critical rate.2. The model should also be precise enough; at the moment 2D-ESTER models are reliable in the mass range 1.8 M -8 M .At lower masses surface convection has an impact on the atmosphere while at higher masses, radiative mass loss also perturbs the atmosphere.Both effects are presently ignored in the standard ESTER models.In addition, the stellar atmosphere model also needs to be closer to reality.The difficulties met at wavelengths just below the Balmer jump prevent us from using the Strömgren u-filter.3. The third limit comes from the fact that we largely use the UV part of the spectrum, which can only be measured from space.4. To compute the temperature with the IRFM approach, we need to know accurately the bolometric flux, which is not usually an easy task.The accuracy of our inversion scheme depends on the precision of the photometry.We implicitely assumed that measurements are accurate at a one percent level.More precise observations will improve the sensitivity to gravity darkening and make the inversion easier.
The next steps in the use of gravity darkening to constrain the fundamental parameters of stars will be to first test other stars for which interferometric data exist together with precise spectrophotometric data.Interferometric measurements are indeed a way of control when the radius is left as a parameter to be determined.Stars like α Oph (M ∼ 2.2 M ), α Leo (M ∼ 4.1 M ), and ζ Aql (M ∼ 2.5 M ) are all rapid rotators seen equator-on with masses appropriate for ESTER models (see Espinosa Lara & Rieutord 2013;Howarth et al. 2023 for their mass determination).They offer the opportunity to test our algorithm with another interesting inclination.
Our method heavily relies on UV data and at the moment these data are from the HST or IUE satellite, but new space A50, page 11 of 12 missions are planned, such as the Ultraviolet Explorer (UVEX; Kulkarni et al. 2021) and the Large Ultraviolet Optical Infrared Surveyor (LUVOIR; The LUVOIR Team 2019).Star clusters are of particular interest because their members are assumed to share the same age and the same initial composition, setting a strong constraint on these parameters.Our method can be used to determine, from simple photometry, the presence of fast rotation among stars in a cluster and its actual impact on their evolutionary stage.The HAYDN mission (Miglio et al. 2021) could bring valuable information about the spin axis distribution among clusters.

Fig. 2 .
Fig.2.Relative radiative flux difference between a plane-parallel and a spherical PHOENIX model with T eff = 20 000 K, log g = 4, and a radius of 2 R , which is the smallest curvature radius for a M = 3 M , ω = 0.7 stellar model.

Fig. 4 .
Fig.4.Values of R th at different wavelengths for log g = 4.The present data are plotted with blue crosses and those ofBlackwell et al. (1980) with red dots.For each wavelength, the effective temperature increases with increasing R th .

Fig. 6 .Fig. 7 .
Fig. 6.Values of Q for the models of our grid (crosses) vs. inclination.Each dashed line corresponds to their average value over masses (M ∈ [3−7] M ) for a given rotation rate, ranging from ω = 0.3 (top) to ω = 0.8 (bottom).The red crosses show the 2 M models.
defined a photometric system consisting of four intermediate-width bands (u, v, b, y).The c 1 index (in magnitudes) is the colour difference (u−v)−(v−b).The transmissions of the filters entering the definition of c 1 are shown in Fig. 7.

Fig. 8 .
Fig. 8. Variations in the c 1 index with the inclination i for several values of ω in the range [0.3;0.7] with a step of 0.1, for a 5 M model (upper panel) and for a 2 M model (lower panel).The dashed curves represent the corresponding fits.

Fig. 11 .
Fig. 11.Models of our grid in the T IRFM -mass plane where each dot represents a model of given (ω, i).The black solid curve shows the fit of the average temperature for each mass.

Fig. 12 .
Fig. 12. Example of a filtered grid (see text) for a 5 M model with ω = 0.5 and i = 65 • .The dots represent models that meet conditions (8), colour-coded by the scattering between the inclinations obtained from each observable, as defined by the left-hand side of Eq. (10).Green crosses stand for initialisation models for which Eq. (10) is verified with ε = 0.1.The big red cross shows the location of our observed star.

Fig. 13 .Fig. 14 .
Fig.13.Error maps in the ω, i plane, for stars with masses equal to 3, 5, and 7 M (top, middle, and bottom panel, respectively).The colour stands for the absolute error between the retrieved and the original models defined as (M sol − M obs ) 2 + (ω sol − ω obs ) 2 + (i sol − i obs ) 2 , where M is in M and i is in radians.The subscripts 'sol' and 'obs' stand for the solution and the original model, respectively.The yellow cells correspond to differences greater than 10 −2 , where the resolution possibly failed.
Fig. 15.Comparison between the spectrum from CALSPEC and from our model for Vega.Top panel: spectra of Vega from CALSPEC and from our model (orange and blue lines, respectively; left scale) and energy transmission curves of the filters used to compute the colour indices (right scale).The Strömgren filters are shown as dashed lines, while the HST filters as solid lines, the colours having the same meaning as in Figs.7 and 9, respectively, except that the u filter is replaced here by u .Bottom panel: difference between the flux of our Vega model and that of CALSPEC (black solid line).

Fig. 16 .
Fig. 16.Influence of the central wavelength of the u filter on the c 1 index.Top panel: difference between the c 1 index determined with our model and that of the CALSPEC spectrum as a function of the central wavelength of the u filter.Bottom panel: spectra of Vega from CALSPEC and from our model (orange and blue lines, respectively).

Table 1 .
Number of latitudinal grid points n θ that ensures a 0.1% precision on T eff and log g for various values of ω.

Table 2 .
Parameters of the ESTER Vega model of Espinosa Lara &