© K. Bouchaud et al. 2020

1. Introduction

A large fraction of intermediate-mass and massive stars have high rotation rates (Zorec & Royer 2012; Ramírez-Agudelo et al. 2013). Understanding the physics of their interiors requires an understanding of how rotation affects them. Indeed, at high initial masses, stars rotate (on average) rapidly, up to several hundred km s⁻¹ for the fastest (Royer 2009). Such high rotation rates, which sometimes bring stars close to break-up, have a strong impact on their internal structure and evolution. For example, unless they possess a strong magnetic field, and early-type stars rarely do, many of these stars exhibit signs of differential rotation (Reiners 2007). This differential rotation induces meridional circulation and small-scale turbulence that carry matter and angular momentum (Mestel 1953). This effect is commonly called “rotational mixing”, and while it brings fresh hydrogen to the core, increasing the time the star spends in the main sequence (MS) phase (some stars may even skip the Blue Loop in the giant phase Georgy et al. 2013), it also brings elements formed in the core of the star up to the surface, where they can be observed. Abnormal abundances are thus expected in fast rotating stars, compared to the slowly rotating ones (e.g. Heger & Langer 2000; Meynet & Maeder 2000). Yet, some fast rotators do not show the expected abnormal abundances, while some slow rotators do show abnormal abundances, as summarised in Cazorla et al. (2017). Obviously, rotational mixing still requires further investigation.

These examples show that taking rotation into account in stellar interior models is paramount to a more accurate description of stellar evolution. Unfortunately, as rotation has two-dimensional (even three-dimensional) effects, it cannot be easily accounted for by 1D models. For instance, the 1D MESA models (e.g. Paxton et al. 2018) include rotational mixing as a pure diffusive process while advection is known to be essential in the transport of angular momentum (Zahn 1992). The difficulty is increased when one has to deal with data that are directly influenced by the non-spherical nature of rotating stars, like spectro-interferometric observables. For example, Domiciano de Souza et al. (2002) performed a detailed study on how physical models can be used to investigate and constrain some model parameters from the observed interferometric signatures of rotation, namely rotational flattening and gravity darkening (GD). From the observational side, many results were obtained from spectro-interferometry, and van Belle (2012) give a review of some of the different attempts made on several stars. Although important advances were achieved thanks to these and other theoretical and observational studies, they are based on simplified models, not taking into account the 2D internal structure of the rotating stars.

Recent progress in programming techniques and computer power has enabled the creation of fully two dimensional stellar models by the ESTER code (Rieutord et al. 2016). ESTER models indeed predict the differential rotation profile and the associated meridional circulation of an early-type star at a given stage of its MS evolution. The solution given by the code is presently the steady state solution of an isolated rotating star. Time evolution has not been implemented yet. However, by tuning the hydrogen mass fraction in the convective core, a good approximation of an evolved state on the MS can be computed. Espinosa Lara & Rieutord (2013) have shown that the fundamental parameters of three nearby rapidly rotating stars obtained with interferometry could be reproduced fairly well. The three stars are α Leo (M ≃ 4.15 M_⊙), α Lyr (M ≃ 2.2 M_⊙) and α Oph (M ≃ 2.4 M_⊙). However, the fundamental parameters of these stars have been derived from interferometric data using simple stellar models, namely Roche models with uniform rotation. Gravity darkening, a non-uniform flux distribution on the stellar surface as a result of rotation, was modelled using the modified von Zeipel’s law , with the exponent β being adjusted. Hence, several approximations of stellar structure are presently used to interpret the interferometric data. One might therefore wonder what would be the fundamental parameters derived from interferometric data if more realistic 2D models were used. This, in turn, would inform us on the validity domain of 1D models and their approximations.

To test this new way of modelling fast rotators, we chose to focus on the well studied A7V star Altair (α Aquilae, HD 187642) for which numerous data sets are available.

To determine its parameters, we first endowed ESTER models with appropriate stellar atmosphere models. We thus obtained spectra and intensity maps of Altair, the latter being constrained with interferometric data. The results gave new fundamental parameters for Altair based on the most up-to-date two-dimensional models. Interestingly, the new models enable an interpretation of the δ-Scuti type oscillations of Altair detected by Buzasi et al. (2005).

The paper is organised as follows: in Sect. 2, we list the observational data that we use. In Sect. 3, the ESTER code is briefly described. We present the preliminary models that match the previous determinations of Altair’s parameters (Sect. 4). We then focus on the combination of atmospheric models with interior models (ESTER and ω-model) and compute interferometric observables and spectra from the resultant intensity maps, to be compared with the observational data (Sect. 5). The model-fitting method we used and the results obtained are presented in Sects. 6 and 7. Finally, we discuss several points of interest (Sect. 8) and give some prospects for the future that this work offers.

2. Data

2.1. Interferometry

Two sets of near-IR interferometric data from two different VLTI beam-combiners were used to study Altair: 8 observations from PIONIER (H band; 1.65 μm) and 6 observations from GRAVITY (K band; 2.2 μm). The PIONIER observations (listed in Table 1) were obtained from the JMMC¹ web service OiDB² and combine two different observing programmes. The data reduction and calibration process for both programmes were done using the PNDRS pipeline (Le Bouquin et al. 2011). The first part is made of five observations done over two nights in September 2011. Seven channels in the H band are available. The second part is composed of three observations done in one night, in October 2014, for the Exozodi survey (the details of the observation are given in Marion et al. 2014). They were conducted with a compact configuration (AT) and three wavelength channels in the H band. This second part of PIONIER data was used despite the three-year gap with the first one, as the compact configuration will only be sensitive to large scales, such as the overall shape of the star, which does not change significantly in such a short time for Altair.

The second set of data is made of two nights of GRAVITY observations (see Table 2). They were obtained during the first Science Verification Time (SVT) of the instrument (programme ID: 60.A-9164(A)). They were reduced using the GRAVITY pipeline (Lapeyrere et al. 2014), with the star HD 188310 used as calibrator, through reduction recipes made available by the GRAVITY consortium in their python toolkit³. Six squared visibilities and four closure phases were obtained with each observation, for both the p and s polarisation directions. This yields two data sets at high resolution in the K band (R ∼ 4000), and two at low resolution for the fringe tracker. As there was no significant difference between the data sets for the two polarisations, we averaged them. Only the science detector data was used in our analysis.

The uv-plane for both instruments is shown in Fig. 2, and the observables are shown in Fig. 1 along with the theoretical data.

Fig. 1. Interferometric observables (squared visibilities and closure phases) versus theoretical predictions from the best model (see Sect. 7). Each colour indicates a different baseline, with the corresponding telescopes shown in legend.

Fig. 2. uv coverage of PIONIER and GRAVITY observations of Altair.

2.2. Spectroscopy

The spectrum used for the spectroscopic analysis of the star was obtained on October 1st, 2003, using the ELODIE instrument at the Observatoire de Haute-Provence. After merging the orders of the initial echelle spectra, Reiners & Royer (2004) combined five single exposures, with a spectral resolution of 42000. They calibrated and normalised it using HD 118623 and γ Boo as calibrators. The spectrum has a signal-to-noise ratio (S/N) of 228 (average S/N for the five exposures), between λ = 3850 and 6800 Å. As no error bars were computed in the process, we averaged the relative errors found in the files of the five single exposures, available in the ELODIE archives⁴. The part of the spectrum used in our analysis is shown in Fig. 3 along with the theoretical line.

Fig. 3. Observed and theoretical MgII line at 4481 Å of the best ESTER model (see Sect. 7. The parameters of the corresponding model are shown in Table 5).

3. ESTER: Modelling stellar interiors in 2D

3.1. Partial differential equations

Several attempts at modelling stars in two dimensions have been made over the last 50 years or so, but all have failed to describe the large-scale dynamics of the interior, with realistic internal rotation. The aim of the ESTER project has therefore been to address this issue, thus leading to the eponymous code. We shall briefly describe the models here, but a more detailed presentation may be found in Espinosa Lara & Rieutord (2013) and Rieutord et al. (2016).

ESTER models describe the (quasi) steady state of an isolated, non-magnetic early-type star made of a convective core and a radiative envelope. The core is assumed to be isentropic, which is very close to what the mixing-length modelling says.

The four partial differential equations that determine such a model are:

(1)

These are Poisson’s equation (ϕ is the gravitational potential, ρ the density, and G the gravitational constant), the equation of entropy S (T is the temperature, v the velocity, F the heat flux, and ε_* the nuclear heat sources), the momentum equation in an inertial frame (P is the pressure and F_v the viscous force), and the equation of mass conservation.

These equations are completed by the expressions of the heat flux, the viscous force, the energy generation, the equation of state and opacities. For these latter quantities OPAL tables are used (Rogers et al. 1996).

3.2. Boundary conditions

The system of equations needs to be completed with appropriate boundary conditions, in particular at the surface: (i) the gravitational potential must match that of a field in vacuum, vanishing at infinity; (ii) the velocity field must match stress-free conditions; and (iii) the temperature must meet the local black body radiation condition.

These conditions are applied on the bounding surface, chosen as the isobar where the pressure is equal to the polar one, defined as

(2)

where τ_s is the optical depth that defines the polar photosphere, g_pole is the polar gravity and κ_pole is the polar opacity (for more details see Espinosa Lara & Rieutord 2013). The surface radius, chosen as the radius of the polar isobar, coincides with the polar photospheric radius, but slightly differs from the actual photosphere radius elsewhere. However, it can be shown (see Appendix A) that as long as the flattening of the star, defined as

(3)

does not exceed 0.2, the difference between the two radii is less than one percent. Interestingly, ε ∼ 0.2 is the value obtained for the models of Altair presented in Sect. 7.

Finally, either the total angular momentum or the surface equatorial velocity V_eq has to be specified.

3.3. Numerical solver

For the numerical part, the star is subdivided into multiple domains in the radial direction. In each domain, a Gauss-Lobatto collocation grid appropriate for spectral methods based on Chebyshev polynomials is used. The advantage of a multi-domain approach is that even for spectral methods, it is possible to increase the resolution where it is needed, in particular near the stellar surface where rapid variations occur. An appropriate set of coordinates, which adapts itself to the centrifugal distortion of the star, is used.

The set of Eq. (1) is then solved over the stellar interior using Newton’s method. This iterative method has a quadratic convergence if the initial solution given in input is not too far from the real one. Usually, a 1D model computed with ESTER is given as input for the computation of a 2D one when the rotation velocity does not exceed 50 to 70 % of the break-up value (the break-up rotation velocity is defined as the velocity at which the centrifugal force and the gravitational force at the equator are equal). Above that rotation rate, an intermediate 2D model must be computed first and given as input.

3.4. Limitations

Apart from the already mentioned assumptions concerning the star, such as the absence of a magnetic field, ESTER models face other limitations. The major one is that convection in surface layers (other than the core) has not yet been implemented, thus stars with a mass below about 1.6 M_⊙ and evolved stars cannot be modelled, as convective layers start to form in these stars just below the surface. Altair is thus a good test for ESTER since it is a well studied star and its mass is close to the lower mass limit manageable by the code. On the high mass side, the limit is around 40 M_⊙ presently (e.g. Gagnier et al. 2019).

4. Altair’s fundamental parameters

Altair is a nearby (5.13 ± 0.02 pc, van Leeuwen 2007) rapidly rotating star that has already been studied extensively in OLBI. van Belle et al. (2001) were the first ones to measure the oblateness of Altair, using data from the Palomar Testbed Interferometer (PTI). They found an oblateness ε (Eq. (3)) of 0.122 ± 0.022, and also derived a projected rotation velocity v sin i, independently of spectroscopic analyses, of 210 ± 13 km s⁻¹. This value falls within the range of velocities determined through different spectral line studies, from 190 km s⁻¹ (Carpenter et al. 1984) to 250 km s⁻¹ (Stoeckley 1968), the latest estimate being 227 ± 11 km s⁻¹ (Reiners & Royer 2004).

Ohishi et al. (2004) then showed evidence of gravity darkening from NPOI data (Navy Prototype Optical Interferometer). Domiciano de Souza et al. (2005) confirmed that Altair was compatible with the von Zeipel relationship between T_eff and g for hot stars (T_eff ∝ g^β, with β = 0.25) after adding VLTI’s VINCI data to the two previously mentioned data sets, analysing the star in different spectral bands. They also gave a broad estimate of the inclination angle i (angle between the polar axis and the line of sight), between 40 and 65°. Peterson et al. (2006) gave a more precise value of 63.9 ± 1.7°.

The latest interferometric study of Altair led to the first “direct” imaging (via image reconstruction techniques), by Monnier et al. (2007), exploiting data from the CHARA array’s instrument MIRC. They found a departure from von Zeipel’s theorem, with a value of β = 0.19 in their models allowing a better fit to the data.

Along with interferometry, other observation techniques give us additional information on Altair. In addition to the projected velocity, Reiners & Royer (2004) determined its inclination to be greater than 68° at a 1σ-level (45° at a 2σ-level). On the asteroseismic side, let us mention the works of Buzasi et al. (2005) and Suárez et al. (2005) who showed with data from the WIRE satellite that Altair exhibits δ Scuti-type pulsations.

Using two-dimensional models from the ESTER code, we aim at further improving the determination of Altair’s fundamental parameters by better modelling the effects of fast rotation. This can be done by comparing the observational data with our models, but it requires replacing the Roche models (simple and fast to create) with the more numerically demanding ESTER models. To test the feasibility of this processing we first derived an ESTER model of Altair by a manual adjustment of its parameters to the parameters found by Monnier et al. (2007). The result given in Table 3 shows that ESTER models can reproduce fairly well the parameters obtained by Monnier et al. (2007). We note that the equatorial velocity is slightly off the error bars and that the derived mass is smaller than the one derived by Monnier et al. Nevertheless, this first attempt encouraged us to take up the challenge of deriving again, ab initio, the fundamental parameters of Altair from the above data using ESTER models.

5. Atmospheric models for ESTER models

The analysis of interferometric and spectroscopic data requires monochromatic intensity maps of the surface of the star from which the complex visibility will be computed. The interferometric observables and the spectrum will then be extracted from it. To obtain these intensity maps, model atmosphere codes coupled with radiative transfer codes must be used.

5.1. Atmosphere models

Among the different existing atmosphere codes, we decided to use models from the PHOENIX code⁵ (Husser et al. 2013), as it is the code that best matches our needs in terms of surface parameters (see Table 3), spectral range, and specific intensity.

Indeed, these spherically symmetric atmospheres were used to produce specific intensities (low resolution, but for different values of the emission angle) and spectra (high resolution), from 500 to 26 000 Å, which suits our study since we are analysing observational data in the visible range (ELODIE), the H band (PIONIER), and the K band (GRAVITY). Models computed for 0.5 < log g < 6.0, 2300 K < T_eff < 12 000 K and −1.0 < [M/H] < 1.0 were selected from the PHOENIX online library.

To compute the intensity maps, as seen from Earth, we need the specific intensities as a function of the emission angle I(μ). These are available in the PHOENIX database, but only with a 1 Å step. This kind of resolution is sufficient to analyse PIONIER and GRAVITY data, but as most lines are blended together by rapid rotation, initial high resolution spectra are needed to obtain an accurate broadened spectrum in the visible, which we need for the analysis of our spectroscopic data. Such high resolution spectra are only available in the database as fluxes, independent of the direction of emission. To compute the desired quantity, we used Claret (2000)’s expression of limb darkening (Eq. (6)), using the a_k-coefficients found in the dedicated Vizier catalogue (Claret 2018, most recent catalogue for the wavelengths of interest). The expression

(4)

can be coupled with

(5)

to yield

(6)

Thus I_λ(μ = 1) can be obtained from the flux F_λ, followed by I(μ), that is the specific intensity emitted in any direction. Let us point out that this method, based on band-integrated limb darkening coefficients, can only be approximate, especially if some thermal convection is present at the surface (Dravins 1982). A more precise method is however premature at this stage of the investigations.

5.2. Linking interior and atmosphere models: ESTERIAS

Once a complete grid of stellar models and a grid of atmosphere models are made, we use them both to create intensity maps. This is done with ESTERIAS (ESTER for Interferometry, Asteroseismology, Spectroscopy), the Python tool we developed to compute interferometric and spectroscopic observables with ESTER.

5.2.1. The stellar surface grid

A grid is used to represent the surface of the star, dividing it in surface elements dS. We chose to construct the grid as “rings” of constant co-latitude θ, each point of a ring being defined by a ϕ value. In order to do this, we impose the number of rings (θ values) from pole to pole, and the number of azimuthal points (ϕ values) of the first ring (from the pole). The radial distance (function of θ) corresponding to every ring (the star being axisymmetric) must then be interpolated from the ESTER model over the θ values chosen for the grid. Indeed, ESTER models only require 20 grid points in latitude to achieve a precise solution. However, this is not enough for an accurate representation of the stellar spectra. Fortunately, the spectral representation with spherical harmonics allows the user to get the value of any parameter at any co-latitude, without losing numerical accuracy. Some tests showed us that the relative difference between a model computed with N_θ = 100 θ-values and parameters interpolated on these 100 co-latitudes from a model made with only N_θ = 10 never exceeded 10⁻⁴.

Once the “radius” and co-latitude of each ring are determined, the number of azimuthal points must be set. To determine the appropriate numbers of points per ring, we first imposed a number of points N_ϕ0 on the first ring. The area dS₀ of these N_ϕ0 surface elements is then computed and is identical for all points on the ring, as dS only depends on θ, dθ and dϕ, see Eq. (7). Finally, N_ϕ is computed for each ring such that dS is as close to dS₀ as possible. This is done using Eq. (A.51) of Rieutord et al. (2016):

(7)

with dS the surface element area, and r and r_θ the radial coordinate and its derivative with respect to the θ coordinate.

Apart from r and r_θ, the relevant stellar surface parameters that must be extracted from the models and interpolated on our θ values are: T_eff, the effective temperature; log g, the logarithm of the effective gravity; and Ω, the angular velocity. The linear velocity v_ϕ = r(θ)Ω(θ) sin θ is then computed.

We are thus left with a set of points, each associated with a set of (r(θ),θ, ϕ) coordinates, and physical parameters (T_eff(θ),log g(θ),Ω(θ)), that define our star surface.

5.2.2. The visible grid

Once the surface grid is determined, it can be used to generate a number of “visible grids”, representing the visible side of the star, as viewed from an earthbound observer, for any value of the inclination i of the star’s rotation axis. For this we associate with each point a parameter μ which is the cosine of the angle α between the normal to the surface and the line of sight (see Fig. 4):

(8)

Fig. 4. Side view of a flattened star. The primed coordinates are the Cartesian coordinates used to describe the plane of the sky, as seen from Earth. (Y′, Z′) are the axes of the sky plane, and X′ the line of sight. The non-primed coordinates are linked to the star, and can be obtained by rotating the primed coordinates around the Y′ axis so that the rotation axis Z is at an inclination angle i from the line of sight. α is the angle between the normal n to the surface at an arbitrary point on the stellar surface and the line of sight. We recall that μ = cos α (cf. Eq. (8)). The star is shaded based on a 3D rendering of its shape, and not on physical effects such as gravity or limb-darkening. Other figures such as Figs. 12 and 13 show such effects.

We then keep only the points for which μ ≥ 0.

Once the “visible grid” is ready, we can compute the geometrical parameters Y′, Z′ and ds_proj, respectively the two Cartesian coordinates projected onto the plane of the sky (see Fig. 4) and the projected surface element, as well as the linear velocity projected onto the line of sight v_proj. We find that:

(9)

At this step in the process, we have a two-dimensional representation of the stellar surface (a 2D projection of a 3D object actually), giving a map of its shape, temperature, gravity, and rotation velocity.

5.2.3. The intensity map

Several operations (mainly interpolations) are necessary to get the final intensity maps from the specific intensity files. First it is necessary to associate the original intensity I_λ from the files, which were computed for fixed values of effective temperature and gravity, with every point on the grid. Then, for each μ associated with each visible grid surface element, the I_λ are interpolated from the original files, or computed with Claret’s formula (Eq. (4)). Finally, the Doppler shift due to the rotation of the star is computed. We thus get one intensity map per wavelength value (monochromatic intensity maps), for which the several observables (e.g. interferometric, spectroscopic) can be computed.

5.2.4. Observable quantities

For the interferometric data, the two components of the sky-projected baselines (B_x, B_y), the squared visibilities and closure phases with their associated errors, and the wavelengths are extracted, and the spatial frequencies are computed. Here, the x and y subscripts refer to the cartesian sky (angular) coordinates, with x positive towards the east and y to the north. A transformation taking into account the position angle of the star (PA, the angle between the northern direction and the axis of rotation projected onto the plane of the sky, counted positive to the east) must be applied on Y′ and Z′ (from Fig. 4) to get x and y (e.g. in Eq. (11)). As the star is made of a series of rings of different sizes, with non-uniform spatial discretisation, a Fast Fourier Transform algorithm cannot be used. We thus computed the discrete Fourier Transform via the classical formulas:

(10)

(11)

d being the distance to the star. The complex visibility, at each λ, is then given by

(12)

with the flux integrated over the visible surface of the star (observed flux). The squared visibilities and closure phases are computed from the complex visibilities with

Ψ_ij being the phase of the complex visibility computed at the projected baseline defined by telescopes i and j, for i ≠ j, and i, j = 1, 2, 3 (three telescopes).

6. Model fitting

We apply a model-fitting approach to find the parameters that best represent the star, using the emcee code (Foreman-Mackey et al. 2013). This Python implementation of the Markov chain Monte Carlo (MCMC) method allows the user to easily tweak the parameters of the fit: the number of values to test at every step, the number of steps, the priors to apply to constrain our fit, etc. The MCMC method consists in drawing a constant number of random “walkers”, i.e sets of values for the free parameters of the problem, within the allowed range for each parameter, and computing the posterior probability density function associated with each set. A new ensemble of values is then randomly drawn, and each walker will be replaced by the new value with some probability, the probability being higher if the new value gives a better fit than the previous one. The process is then repeated for a preset number of steps, or until convergence is met. Indeed, as the process proceeds, the walkers will be more and more closely gathered around the best-fitting solution(s). The randomness of the process enables to detect and avoid local minima, at the condition that the parameter space be sufficiently sampled (by choosing a high enough number of walkers, and the right parameter space domain).

Unfortunately, as some parameters become more “extreme” (mass below 2 M_⊙, low or high metallicity, high rotation velocity, etc.), ESTER will fail to create the model if the input model is not really close to the output solution. A step-by-step approach, creating intermediate models before reaching the desired one, can be used, but not as part of an MCMC model-fitting. Indeed, as there is no way to know beforehand whether ESTER will converge on a solution, this step-by-step approach has to be done manually, changing the increment, or the input parameter to increment, if needed. Furthermore, computing ESTER models is very time-consuming, as one model takes between 30 s and several minutes to produce. This may be fast for the computation of the full two-dimensional structure of a star, but isn’t appropriate for a model-fitting needing the computation of several tens of thousands of models.

We thus decided to first create a grid of stellar models with ESTER, with fixed increments in parameter values. The parameter space covered was chosen by taking into account either the parameter values found in the literature (when available) or various assumptions concerning plausible ranges for other parameters (X_c for example). Then, ESTERIAS was made so that when random values of parameters are drawn by the MCMC routine, the parameters which are relevant for our star surface (effective temperature and gravity, radius, and rotation velocity) are linearly interpolated for these values from the closest models in our grid (selected via Delaunay triangulation).

We present in the next section the different attempts at determining the following parameters of Altair: mass (M), equatorial radius (R_eq), angular velocity (Ω, expressed as a fraction of the Keplerian angular velocity⁶, ), metallicity (Z), hydrogen mass fraction in the core (X_c), inclination (i), position angle (PA), and the metallicity of the atmosphere models ([M/H]). The envelope hydrogen content X is also considered, as it turns out it plays an important role in the convergence towards a physically realistic solution.

7. Results

7.1. Fitting interferometric and spectroscopic data

First, our attempts to determine Altair’s parameters by fitting interferometric and spectroscopic data are described. Additional constraints obtained through other methods are presented in the following sections.

7.1.1. ω-model

Initially, an attempt was made to constrain simultaneously all the parameters previously mentioned, using ESTER models and only interferometric data. The χ² was computed in a way that the two sets of data (PIONIER and GRAVITY) had the same weight. Convergence could not be achieved, as there were several correlations between parameters, such as M, Z and X_c, or Ω and i, and constraining all seven parameters with only interferometry was impossible. To disentangle this, we decided to study the parameters in smaller groups. As ESTER models need all the parameters as input, and are quite heavy to use, they were replaced by the so-called “ω-model” (Espinosa Lara & Rieutord 2011; Rieutord 2016) for this first step. This gravity-darkening model provides all of the relevant surface parameters (T_eff, g_eff, r, etc.) without going through the difficult process of computing the whole internal structure of the star (a Roche-model is used to compute the effective gravity). It stems from simple assumptions, namely energy is conserved, there are no energy sources in the envelope, and the radiative flux is anti-parallel to the effective gravity. Thus, the resulting models, while being in almost-perfect agreement with ESTER models (see Fig. 5 from Rieutord 2016 and Fig. 3 from Espinosa Lara & Rieutord 2011), do not allow us to adjust internal parameters such as the metallicity or hydrogen fraction in the core. They are still suitable for studying the shape of the star and the distribution of flux on the surface though, making it perfect for the analysis of the interferometric data, which is mostly sensitive to these two signatures of rotation. Such an analysis was already done by Domiciano de Souza et al. (2018), who succeeded in reproducing the interferometric observables of the star Sargas, an evolved 5.1 M_⊙ rapidly rotating star, with the use of the ω-model.

The parameters we are constraining with interferometric data and the ω-model are the equatorial radius, R_eq, the angular velocity, Ω (again, here, as a fraction of the Keplerian equatorial velocity), the inclination, i, and the position angle, PA. The mass, M, and luminosity, L, must also be given as input, but as they only affect the scale of T_eff and g_eff, not their distribution, the interferometric observables won’t be sensitive to those. We used M = 1.8 M_⊙ and L = 10.6 L_⊙, both taken from Peterson et al. (2006). Altair’s distance was fixed to 5.13 pc, from the HIPPARCOS parallax of 194.95 ± 0.57 milliarcsecond (mas). Clear convergence was obtained for R_eq and PA (see Fig. 5). For Ω and i, convergence is still achieved, but an effect due to the resolution of the intensity maps is clearly visible. This is discussed in Sect. 8. The envelope of the peaks in the histograms of Fig. 5 is of Gaussian shape, with a relatively small width (judging from the ±1σ values).

Fig. 5. Corner plot showing the convergence of the MCMC model-fitting for Req, Ω, i, and PA, using interferometric observations (PIONIER + GRAVITY) and the ω-model. The run was made with 200 steps and 100 walkers. Only the last 50% steps are displayed here as convergence was met within the first 100 steps (“burn-in” phase). These results are listed in Table 5.

In the above MCMC model-fitting, we used model atmospheres with a fixed value of the metallicity, namely [M/H] = 0.3. However, one may wonder what impact this metallicity can have on the fitting process. Accordingly, we computed the χ² value of the fit to the interferometric data with only Ω and i as free parameters, for different values of the metallicity of the atmosphere models. We recall that the metallicity of the PHOENIX atmospheres is defined as

(13)

where n is the number density of elements in the star, and M the sum of all metals. As seen in Fig. 6, showing the i and Ω coordinates of the χ² minima, the location of the best fit (in terms of (Ω, i) coordinates) hardly changes. This means that the fitting of the interferometric data is nearly independent of this parameter. Spectroscopy is the go-to method to determine the metallicity of a star’s atmosphere. Yet, only an extensive analysis of Altair’s spectrum would allow us to get an accurate estimate of its composition. Indeed, at such a high rotation rate, neighbouring metallic lines are blended together. Thus, the integrated spectrum will not only depend on the lines’ depth, but also on the list of lines included in the computation of the opacities and on the relative abundances of all elements in the atmosphere. These effects may be responsible for a mismatch between theoretical and observed spectra. Such work needs a dedicated study and is beyond the scope of this paper. We settled on using the line of MgII at 4481 Å to get an appropriate atmosphere metallicity for the next steps. This line is strong and isolated enough that blending is not too much of an issue, and was deemed by Royer (2009) as a valid candidate for v sin i measurement in an A7-type star, which will give us one more constraint on this parameter. The line depth will give us an estimate on [M/H], keeping in mind that the abundance ratios in PHOENIX atmospheres are solar-like. We computed a high resolution spectrum in the range of the MgII line, and found that a value of [M/H] = +0.45 gave the best concordance with the observed spectrum (the integrated theoretical line is identical to the one obtained with our best ESTER model, shown in Fig. 3). We also computed v sin i, and found ∼238 km s⁻¹, at the limits of the range 227 ± 11 km s⁻¹ found by Reiners & Royer (2004).

Fig. 6. Positions in the (i, Ω) plane of the χ2 minima when fitting interferometric data with Ω and i as free parameters, using M = 1.80 M⊙, and several values of [M/H] for the atmosphere models. [M/H] goes from −0.5 to 0.9 with a 0.1 step. The grey rectangle shows the area enclosed within the error on i and Ω shown in Fig. 5.

7.1.2. ESTER models

We now want to determine the mass M, the metallicity Z, and hydrogen mass fractions in the envelope and core, X and X_c respectively. ESTER models are used, with i, PA and Ω fixed to the previously obtained values. Interferometric and spectroscopic data are analysed simultaneously.

The metallicity of the stellar model must be addressed here. Ideally, this metallicity should match that of the model atmospheres unless some physical phenomena, such as diffusion, were to lead to a different composition in the surface layers. In the present case, we treat the two metallicities as independent parameters as they are subject to different constraints. Indeed, the atmospheric metallicity is especially sensitive to spectroscopic constraints as shown above. In contrast, the stellar model bulk metallicity is subject to both interferometric and spectroscopic constraints due to its impact on the stellar structure. Furthermore, it is correlated with M, X, and X_c. Indeed, when the mass decreases, the size of the star also decreases; but increasing Z increases its size, thereby compensating the effects of the mass decrease. Decreasing X_c (for a fixed X) mimics a more evolved star and, therefore, also increases its size. The dependence of stellar size on X is more complex. In order to circumvent these difficulties, we decided to separate the parameters, and search for Z and X_c simultaneously, for different values of the mass and envelope hydrogen mass fraction. We chose two values of X, namely X = 0.700 (Grevesse & Noels 1993) and X = 0.739 (Asplund et al. 2005). Figure 7 shows the resulting χ² maps. For both values of X, X_c is such that 0.5 <  X_c/X <  1.0. We started with X = 0.700, covering masses from 1.70 to 1.80 M_⊙, expecting Altair to be in this range, between our first estimate of Sect. 4 and Peterson et al. (2006)’s value. Then, as an asteroseismic study of these models favoured higher masses (see Sect. 7.3), we opted for masses between 1.75 and 1.85 M_⊙ for X = 0.739. A mass of 1.90 M_⊙ or higher would lead to a higher hydrogen mass fraction in the core than in the envelope, which wouldn’t make sense with regard to current theories of stellar formation and evolution. Figure 8 shows cuts at fixed Z values (upper row), and along the valley with low χ² values (lower row), for X = 0.739. The figure would be the same for X = 0.700, with χ² minima shifted towards higher values of X_c/X.

Fig. 7. χ2 maps of the fitting of ESTER models to both interferometric and spectroscopic data. The upper and lower rows correspond to X = 0.739 and 0.700, respectively. For each mass and X value, only Z and Xc vary, while every other parameter is fixed to the values previously obtained (see Fig. 5). The decimal logarithm of the reduced χ2 is shown in colour. The red dot in each subpanel corresponds to the minimum χ2 value.

Fig. 8. χ2 cuts at fixed Z values (upper row), and along the valley with low χ2 values (lower row), for X = 0.739. The points in the lower row correspond to the minima in the upper row. As can be seen, a definite minimum appears even within the valley with low χ2 values. Z values are shown in legend, and only the values for which the optimal Xc/X is in the range [0.5−1.0] are shown in the lower row.

For a fixed value of X, the loci (in (Z, X_c) coordinates) of the χ² minima as a function of the mass (red points in Fig. 7) show almost linear relations for Z = f(M) and X_c = f(M), but a degree 2 polynomial seems to fit them better. We thus obtain two degree 2 polynomials as a function of M for Z, and two for X_c. We then do a linear interpolation on the coefficients of these polynomials between X = 0.700 and X = 0.739, and obtain a single expression for Z, as a function of both M and X, and likewise for X_c:

(14)

and

(15)

We may now check whether these relations allow us to retrieve the parameters we had found in Sect. 4 by matching the temperature and radius of Monnier et al. (2007). Extrapolating these relations to a lower mass of M = 1.65 M_⊙ and an X value of 0.700, we get Z ∼ 0.008, and X_c ∼ 0.30. This is somewhat different from our first estimate (see Fig. 3). This is not surprising, as our results using both the ω-model and ESTER give a slightly smaller size and broader surface temperature range than Monnier et al. (2007) (this is discussed in Sect. 8).

If we can now obtain a value for the mass and hydrogen mass fraction through other means, we will get Z and X_c at the same time. In what follows, we investigate what constraints can be placed on the mass.

7.2. Spectral energy distribution

One solution to determine the mass is to compare the spectral energy distribution (SED) of the model with the observed SED of Altair. The SED, which is a measure of the energy received on Earth from the star at different wavelengths, increases as the mass increases (as it mainly depends on the effective temperature of the star, and its radius). But this is only true when the other physical parameters of the star are kept constant. The correlation between M, Z, X, and X_c makes it so that the models which follow the previously found relation have a similar temperature range and distribution (this includes the size and shape of the model star), despite having a different mass. Their SED should thus be fairly similar, the higher-mass models nonetheless having a slightly lower mean temperature (∼270 K decrease from M = 1.70 M_⊙ to M = 1.80 M_⊙, for X = 0.700, and ∼170 K decrease from M = 1.75 M_⊙ to M = 1.85 M_⊙, for X = 0.739). To test this hypothesis, we computed the SED of three models for X = 0.700 (corresponding to four different masses), and 3 models for X = 0.739, with Z and X_c following the relations previously found. The result is shown in Fig. 9. The observational data used is available on Altair’s page of the Vizier photometry viewer⁷. Apart from the blue part of the curve, where the effect of the temperature difference is visible, the SED are nearly identical.

Fig. 9. Spectral energy distribution produced in the direction of the observer (the inclination of the star being the one shown in Table 5), for different ESTER models with X = 0.700 and 0.739. is the mean surface temperature of the models. The blue dots correspond to the observed SED data from the Vizier website.

Comparing the projected rotation velocities of the models did not help, as the v sin i found for these 4 ESTER models ranges from 229 to 238 km s⁻¹ with increasing mass, still within the error bars of Reiners & Royer (2004).

7.3. The word of asteroseismology

Since Altair is a δ Scuti, an alternate approach to constraining its mass is to model its observed pulsation spectrum. Indeed, acoustic pulsation modes (p modes) potentially provide a tight constraint on the mean density (e.g. Reese et al. 2008; García Hernández et al. 2015), which when combined with the volume, provides the mass. δ Scuti type pulsations were discovered in Altair thanks to the star tracker on the WIRE satellite (Buzasi et al. 2005). The pulsation frequencies and amplitudes are reproduced in Table 4 for convenience, with the first seven being arranged in order of increasing frequency, while keeping the original indices of Buzasi et al. 2005, who ordered the modes by decreasing amplitude. Suárez et al. (2005) subsequently attempted to interpret Altair’s pulsations using a second order perturbative approach to model the effects of rotation on the pulsations. In the present study, we use a full 2D approach, both for stellar structure and pulsations, which is necessary for a star rotating as rapidly as Altair (e.g. Reese et al. 2006).

One of the difficulties with which Suárez et al. (2005) were confronted and which still remains in the current study is the fact that the mode identification, i.e. the correspondence between observed and theoretically calculated modes (as characterised through a set of quantum numbers), is unknown. Rather than carrying out an exhaustive search for best fitting pulsation spectra using a χ² minimisation, we prefer to make a few simplifying assumptions and see to what conclusions they lead. First of all, the frequencies f₁, f₂, f₇, f₃, and f₆ form a regular pattern with subsequent frequencies being spaced by roughly 5 or 2.5 cycles/day. Furthermore, their amplitudes alternate between higher and lower values. A compelling interpretation is to assume that these are a sequence of island acoustic modes with the same and m values and consecutive values, with a mode missing between f₁ and f₂ (see e.g. Lignières & Georgeot 2009; Pasek et al. 2012 for an explanation on island modes and associated quantum numbers). The alternating amplitudes could partially be explained by the fact that modes with even values are symmetric with respect to the equator whereas those with odd values are antisymmetric. Indeed, for a given mode normalisation, this would lead to different apparent pulsation amplitudes (or mode visibilities) when integrating the intensity fluctuations over the visible disk because of differing degrees of cancellation.

Figure 10 shows a comparison between observed and theoretical frequencies for six models that satisfy the interferometric and spectroscopic constraints. For the theoretical modes, we used which was the simplest assumption to make and leads to selecting the most visible modes. In terms of spherical quantum numbers, these would be the rotating equivalents to ℓ = 0 and ℓ = 1 modes. Furthermore, we calculated disk-integrated mode visibilities using the approach given in Reese et al. (2013), the inclination obtained from interferometry (i.e. i = 50.65°), and a photometric band deduced from Fig. 8 of Buzasi (2004). The modes were normalised such that the maximal Lagrangian displacement times the square of the frequency is kept constant. In Reese et al. (2013), this was found to keep approximately constant visibilities for island modes with the same (ℓ,m) values and to penalise gravity modes which tend to have a low surface amplitude. As can be seen, the alternating amplitudes are correctly reproduced in most cases, at least qualitatively, by the alternating mode visibilities, the most visible modes being symmetric with respect to the equator, i.e. their value is even (or ℓ = 0).

Fig. 10. Observed and theoretical pulsations for Altair. The vertical grey lines that span the plots are the observed pulsations (modes f1 to f7). Their thicknesses are proportional to the observed amplitudes. The blue vertical segments are theoretical pulsations for six different ESTER models which satisfy Eqs. (14) and (15). The theoretical modes are successive m = 0, (or equivalently ℓ = 0 and 1) modes apart from the first mode (see text for details). The lengths of these segments are proportional to the disk-integrated mode visibilities. Left plot: theoretical frequencies whereas the right plot shows the same set of frequencies after applying a suitable homology scaling (e.g. Kippenhahn & Weigert 1990), thus leading to the scaled masses indicated on the figure (assuming the equatorial radius is fixed).

A cluster of modes is shown around 15 c/d. This is where the fundamental mode is expected. We note that Suárez et al. (2005) also found solutions for which f₁ was the fundamental. In our case, it was difficult to precisely identify which of the numerically calculated modes corresponds to the fundamental mode since a large number of gravity modes (g modes) also appear in the same frequency range thus leading to multiple mode interactions and an unclear mode geometry. Indeed, rotation causes acoustic modes, which are located in the envelope, to decrease in frequency due to the increase in equatorial radius, whereas the gravity modes, which are located more deeply in the star, are less affected. This causes the two mode domains to overlap, in particular where the fundamental mode is located. The modes shown in Fig. 10 have large visibilities and a relatively simple geometry in the outer portion of the star. This interaction between g modes and the fundamental mode might explain why there are three observed pulsations in this region. In any case, one can easily exclude the possibility that f₁, f₄, and f₅ are a rotational multiplet. Indeed, the rotation rate as based on the above best-fitting models is around 2.9 c/d which is much larger than the frequency separations between these modes.

Overall, the model frequencies are somewhat too low (see left panel of Fig. 10). This is an indication that the mean densities of the models are too low, provided the mode identification is correct. We therefore scaled the models using a homologous transformation to see what masses (assuming the equatorial radius is fixed as given by interferometry) would lead to a better match with the observations (specifically, we fitted the three upper pulsation frequencies). The scaled frequencies and masses are indicated in the right panel of Fig. 10. As can be seen, the scaled frequencies match the observations fairly well, and the scaled masses converge towards 1.86 − 1.89 M_⊙. For X = 0.700, due to the correlation between M and X_c (see Eq. (15)), such a mass leads to X_c/X slightly larger than unity, i.e. a core with slightly more hydrogen than the envelope. From a physical point of view, such a solution must of course be rejected. For the models at X = 0.739, however, X_c/X remains below unity. We therefore searched for the best fitting model with X = 0.739, thus obtaining M = 1.863 M_⊙. The corresponding values for Z and X_c are shown in Table 5. A full theoretical pulsation spectrum is displayed in Fig. 11 where the mode visibilities are calculated using the same normalisation as above. As can be seen, the islands modes (highlighted in red) are a good match to the observed frequencies, especially at higher frequencies. At lower frequencies they seem to be more strongly affected by avoided crossings, thus potentially explaining the offsets between these and the observed pulsations. At high frequencies, the island modes have the highest visilibities, whereas at low frequencies a large number of gravito-inertial modes with non-negligible visibilities are present.

Fig. 11. Theoretical pulsation spectrum for the M = 1.863 M⊙, X = 0.739 model. The dark grey lines that vertically span the plot are the observed pulsations (modes f1 to f7). Their thicknesses are proportional to the observed amplitudes. The light grey lines of different heights are the theoretical pulsations computed from the model. The red line segments correspond to identified island (or mixed gravito-island) modes, assumed to correspond to the observed modes. Meridional cross-sections of these 6 modes, labelled with the letters “a” to “f”, are shown in Fig. B.1. We note that no island pulsation mode was clearly identified around 18 c/d, probably as a result of an avoided crossing.

7.4. Final result

Based on the results presented in the previous sections, we give in Table 5 the parameters of our best model along with the values obtained by Monnier et al. (2007). These results and their uncertainties are discussed in the next section.

Surface maps of several relevant parameters of this model are shown in Figs. 12–16. The temperature map and the associated monochromatic intensity (in the continuum at one arbitrary wavelength in the H band) are shown in Figs. 12 and 13. Figures 14–16 show the internal and surface rotation profiles of the model star, while Fig. 17 shows its surface velocity projected onto the line of sight. Interestingly, the slight curvature in the lines of constant projected velocity is caused by the differential rotation. The fit of the interferometric observables was already displayed in Fig. 1, while Fig. 3 shows the fit to the spectrum.

Fig. 12. Surface map of the effective temperature of our best ESTER model (parameters in Table 5). The dashed line marks the equator. The values are in Kelvin.

Fig. 13. Monochromatic intensity map of our best ESTER model (parameters in Table 5), at 1.5 μm in the H band. The values are in erg s−1 cm−2 cm−1 srad−1.

Fig. 14. Meridional cut of our best ESTER model (parameters in Table 5). The colours represent the angular rotation rate.

Fig. 15. Surface map of the rotation rate of our best ESTER model (parameters in Table 5). The colour represents the rotation period (in hours).

Fig. 16. Surface rotation rate as a function of the co-latitude, for our best ESTER model (parameters in Table 5).

Fig. 17. Surface map of the velocity of our best ESTER model (parameters in Table 5) projected onto the line of sight. The values are in km s−1.

8. Discussion

8.1. Obtained parameters

We discuss here the parameter values obtained for our best ESTER model, shown in Table 5, and how they compare with previous estimates.

Radius. Through our analysis of interferometric data, we found for the equatorial and polar radii, R_eq and R_pole, smaller values than both Domiciano de Souza et al. (2005) and Monnier et al. (2007; hereafter D05 and M07, respectively), outside of their error bars. Our polar radius especially, being 4% smaller than that of M07, induces a higher flattening of the star in our best model (∼0.22 against ∼0.19). This translates into a smaller apparent angular diameter in the polar direction of mas, compared to ∼3.30 mas for D05 (M07 only gave an equivalent angular radius to the polar radius of their model, not the actual observed apparent angular radius in the polar direction). One can make the assumption that this effect comes directly from the difference in intensity distribution over the stellar surface compared with previous studies, which stems from the more realistic physics in our models (of great importance here is the modelling of GD). This could also explain our slightly broader range of temperatures, while our mean temperature of 7594 K agrees well with Erspamer & North (2003)’s value of 7550 K.

To compute the error ΔR_eq on the equatorial radius, we proceeded as follows: as the distance Earth–Altair was fixed in our MCMC run, the resulting ±1σ error on R_eq (given by our MCMC run) corresponds in fact to an error on the angular radius ⊘_eq. The ±0.001 R_⊙ error on R_eq thus corresponds to an error on ⊘_eq of ±0.001 mas. Coupled with the 0.3% error on the distance given by Hipparcos (quadratically adding them), we get a relative error on the linear radius of 0.3%. This error is high enough to cover for the combined effects of the wavelength calibration accuracy of PIONIER (∼0.4%) and GRAVITY (a few ∼0.01%). The error on the polar radius directly comes from its definition in the ω-model, as a function of R_eq and Ω.

Position angle. Also based on the interferometry, the position angle was accurately constrained to a value of 301.13 ± 0.34°, which is slightly above the value of 298.2° obtained by M07. D05 found 298 ± 17° from the analysis of NPOI closure phases and squared visibilities from PTI and VINCI (Table 3, BMIRCP column), agreeing nicely with our result.

Ω, i, and v sin i. Our v sin i agrees very well with that of M07, with 243 km s⁻¹ against 240 km s⁻¹ for them. Yet, the individual values of the inclination angle i and rotation velocity Ω are heavily dependent on the brightness distribution over the surface of the star. As stated above, this distribution differs in our models from those of previous studies, leading to a lower inclination and higher rotation velocity, with the inclination still within the 2σ estimate of Reiners & Royer (2004), that is i >  45°. Errors on i, Ω, and PA are the ±1σ tolerance on the MCMC model-fitting with the ω-model.

Mass. Our mass determination (M = 1.863 M_⊙) is significantly higher than previous determinations, with D05 citing Malagnini & Morossi (1990)’s estimate of 1.80 M_⊙, and M07 using Peterson et al. (2006)’s value of 1.791 M_⊙. Both masses were obtained by looking for a 1D non-rotating Geneva model (Schaller et al. 1992) which would reproduce the position of the star in the HR diagram for the former, and the corrected luminosity and polar radius estimates of the latter. Finding a ∼4% difference in mass with their estimates, through the seismic study of a 2D stellar model which reproduces well the interferometric data, is not surprising. The error on the mass was determined through asteroseimology. As was done in Sect. 7.3, we computed the scaled masses needed to match the observed pulsation frequencies for models corresponding to the upper and lower boundaries for Z and X_c, for all M and X values shown in Fig. 7. The resulting scaled masses were all between 1.84 and 1.89 M_⊙, that is at most 0.03 M_⊙ away from the value of M = 1.863 M_⊙, which is the mass that comes out when studying the pulsation frequencies of models which verify relations (14) and (15), as previously stated. This is the value we adopted for the error on the mass.

Temperature. The error on the temperature is trickier. Computing the temperature of the models at M = 1.86 ± 0.03 M_⊙ with Z and X_c following relations (14) and (15), for both X = 0.700 and 0.739, we get errors on the equatorial and polar temperature of about 40 K and 55 K respectively. Yet, this only accounts for the uncertainty on the mass. We can estimate an uncertainty on Z and X_c by multiplying their standard deviation (computed as the square root of the diagonal elements of the covariance matrix) by the square root of the minimum of the reduced χ², as is often done to account for the dispersion of the data. While not ideal, this method gives an idea of the uncertainties on the model parameters as a result of the dispersion of the observations with respect to the fit. This gives us a ∼0.008 uncertainty on Z, and ∼0.185 on X_c. We may only compare the retained model (see Table 5) with the one corresponding to the lower bounds on Z and X_c, as the upper bound on X_c leads to X_c/X >  1. Doing that, we get a difference of about 640 K at the equator, and 810 K at the pole. This again highlights the need for an accurate determination of Altair’s surface composition, as more constraints on Z and X would greatly help reduce this 9.4% error on the temperature.

Z, X, and X_c. This considerable uncertainty on Z, X, and X_c led us not to give the error on these parameters in Table 5. However, our results all point towards the fact that, whatever Altair’s composition is, as long as it is in the range of compositions found for stars in the vicinity of the Sun (where Altair is located), Altair is a young star, close to the zero-age main sequence (ZAMS). This conclusion actually agrees with the latest estimate of Altair’s age made by Peterson et al. (2006) who also place Altair close to the ZAMS. We used the CESAM code (Morel & Lebreton 2008) to get an estimate of the time it would take for a 1.86 M_⊙ star with an initial hydrogen content X = 0.739 to evolve from X_c = 0.739 to X_c = 0.710. The resulting age is ∼93 Ma. Since this age was evaluated with a 1D model taken outside its range of validity, as far as rotation is concerned, we estimate the age of Altair to be around 100 Myr.

[M/H]. We note that the metallicity obtained for the best ESTER model does not match that of the model atmosphere used to reproduce the observed spectrum. Indeed, the correlation between M and Z (see Eq. (14)) leads to Z = 0.019 for X = 0.739. Knowing that

(16)

we get a corresponding [M/H] = 0.19 using the solar composition from Asplund et al. (2005). This is lower than the metallicity used in the model atmosphere (i.e. [M/H] = 0.45). As this mismatch is not entirely satisfactory, we carried out preliminary calculations with an atmospheric metallicity matching the bulk metallicity. This led to a slightly worse fit to the spectroscopic data and incompatible results with the interferometric constraints. Remembering that the metallicity for the atmosphere was obtained by fitting the spectroscopic data for a single absorption line of MgII only, more work is needed to properly determine the atmospheric abundances of Altair. It may indeed be possible that Altair’s atmosphere is more metallic than its interior, as a possible consequence of a recent accretion of metals from a residual disc or planetoids. This idea is comforted by the fact that Nuñez et al. (2017) found an extended, weak IR excess (a few % in K band) for Altair, which suggests the presence of a tenuous circumstellar material (possibly from a debris disc) within a few AU of the star.

8.2. Surface convection at the equator

We mentioned in Sect. 3 that ESTER models do not include surface convective layers. This is not a problem for computing the structure of intermediate-mass and massive stars, since the surface convective layers are very thin and convection carries a small fraction of the flux. However, we can still have an idea of the shapes of convective layers from the distribution of the squared Brunt-Väisälä frequency. In Fig. 18 we show this distribution, which reveals two thermally unstable layers, corresponding to the hydrogen ionisation peak (reaching the surface) and helium first ionisation (below the surface). Convective layers clearly thicken at the equator, as can also be noticed in Espinosa Lara & Rieutord (2013)’s model of Vega. This thickening may be related to the X-ray coronal emission of Altair, which seems to be concentrated around the equator (Robrade & Schmitt 2009). Most probably, this activity is also related to Altair’s observed UV emission (Redfield & Linsky 2005). Finally, surface thermal convection likely alters the broadening of the lines, through micro and macro-turbulence, which in turn influences our fitting of Altair’s spectrum.

Fig. 18. Squared Brunt-Väisälä frequency N2 of the final model, as a function of the radius r and colatitude θ. The thermally unstable layers (regions where N2 <  0), are colour-coded. These regions are bounded by solid lines, with a dashed line showing the minimum of N2 in this particular region.

8.3. Resolution of the intensity maps

The resolution of the grid used to compute the intensity maps turns out to be important. Indeed, we first decided to use only N_θ = 25 sectors in latitude, from pole to pole, to speed up the process. Even with such a low number of points on the visible surface of the star (876), the largest surface element is about 0.02 squared milliarcsecond, whereas the resolution achieved with the longest baseline of the GRAVITY configuration used, i.e. 130 m at ∼2 μm, is about 3 mas. No difference could be seen on the interferometric observables, and even though a spectrum done with a “low resolution” grid (N_θ = 25) was more noisy than the one done at “high resolution” (N_θ = 100), when smoothed (e.g. with a moving average), both spectra were similar (between 0.1 and 0.3% difference). Despite all this, convergence was not reached for Ω and i with N_θ = 25, even though preferred values seemed to come out of the fitting process. On the other hand, using N_θ = 100 allowed convergence towards a solution for all four parameters (see Fig. 5). Unfortunately, for Ω and i, the solution seems to be multi-modal, as can be seen in the histograms of Ω and i in Fig. 5. This effect is more important at low resolution, which suggests that a higher resolution should solve the problem. We thus run the MCMC code with N_θ = 60, 80, 100, 120, obtaining 2.1, 1.6, 1.3, and 1.0 degree modulations on the inclination, respectively. This comforts us in the idea that the resolution in latitude is the main effect at play here, but as the most probable values are identical for all resolutions, we settled on a resolution of N_θ = 100. This gives a high precision on all parameters (inclination included, considering the small dispersion of the peaks) while requiring reasonable numerical resources.

8.4. Influence of the geometry of atmosphere models

As observations became more and more accurate, limb-darkening laws derived from plane-parallel model atmospheres provided a less satisfactory fit to the observed data (Fields et al. 2003; Barros et al. 2012). Sphericity has thus been used extensively, in recent years, in atmosphere codes, to better take into account the actual shape of a (non-rotating) star (spectra made from MARCS atmosphere models are available in both plane-parallel and spherical geometries, as is the case for PHOENIX models). Yet, when rotation comes into play, the question of the curvature radius of the model atmospheres arises. Indeed, as the rotation velocity of the star increases, so does its flattening. This gives a number of different curvature radii of the surface, as a function of latitude (the star is considered, at least in ESTER, axisymmetric). This effect is further complicated by the fact that the curvature radius is not limited to a single value at each latitude on the star, but depends on the direction one considers. Expressions of the limits on the values of the curvature radius at any point on the surface can be obtained from Harris (2006). Their Eq. (18) shows that the minimal and maximal values for such a radius correspond to those computed along the latitudinal and azimuthal directions. Adapting their Eq. (15) to suit ESTER’s spheroidal coordinates, we find that the curvature radius in the θ direction is:

(17)

where r_θθ is the second derivative of r with respect to θ. Likewise, the radius of curvature in the ϕ direction is:

(18)

This allows us to compute the extrema of the curvature radius from pole to pole, and compare it to the actual radius of the star (Fig. 19). The figure confirms that (i) the curvature radii along the θ and ϕ directions are equal at the pole (as expected, since the polar axis is the axis of symmetry of the star), and (ii) the curvature radius along the azimuth is equal to the actual radius at the equator (also expected, as the equatorial plane is a plane of symmetry of the star).

Fig. 19. Curvature radii along the latitudinal () and azimuthal () directions, and the actual radius of the model star r(θ). The vertical line marks the equator.

If we take our ESTER model that best describes Altair, we find a curvature radius at the pole of R_c = 2.08 R_⊙, and radii of = 0.90 R_⊙ and = 2.01 R_⊙ at the equator. As the temperature and gravity distributions at the surface of the star are not uniform, the radii of the model atmospheres associated with the different points on the grid vary, from R_⊙, to R_⊙. We see that, at the equator, the effective radius chosen in the atmosphere models is off the upper limit by about 1 R_⊙. This and the fact itself that no point on the surface of the star (except for the poles) can be associated with a unique curvature radius make us wonder whether the need for atmospheres computed with spheroidal coordinates arises when dealing with rapid rotators.

It has been shown that the importance of the geometry of the atmosphere depends primarily on the extension of the atmosphere compared to the radius of the photosphere (see Plez 1990 for example). Baschek et al. (1991) formulated this extension ΔR/R_photosphere as being proportional to the pressure scale height, H_P, thus inversely proportional to the effective gravity g. Effectively, Heiter & Eriksson (2006) determined that for log g >  3.0, the estimated abundances of a number of elements are similar in both the plane-parallel and spherical geometries. Neilson & Lester (2013) found discrepancies in limb-darkening and apparent diameters, even for stars with compact atmospheres (log g ≥ 4.0), concluding that even in the case of “main sequence stars with large gravities and small atmospheric extensions”, one should use spherical model atmospheres as they are more physically representative when trying to measure precise angular diameters and fundamental parameters of stars from optical interferometry. Yet, this effect only amounts to 1% at maximum in the K-band, way below our uncertainty on the radius of the star. Furthermore, they also find that the difference in gravity darkening between both geometries is negligible for this kind of stars.

A star like Altair has an effective gravity at the surface log g >  3.0, with our best model showing 3.8 ≲ log g ≲ 4.3 from equator to pole. We should then be safe from any significant effect of the geometry of the atmospheres we used, and a spheroidal code for computing model atmospheres would not improve the accuracy of our results in a way that would justify the time and effort put into it. If surface convection and time evolution are one day successfully implemented into ESTER, then modelling evolved stars with extended atmospheres would become feasible, and this matter would have to be resolved before a work such as this one should be considered.

9. Conclusions

We conducted a multi-technique analysis of the star Altair (HD 187642) using interferometry, spectroscopy, and seismology. For the first time, this kind of analysis was performed by comparing observational data with intensity maps computed from a full two-dimensional model of the stellar interior (ESTER) and surface (atmosphere models).

PHOENIX atmosphere models were used when fitting interferometric observables, and were furthermore supplemented by Claret’s 4-coefficients limb-darkening law when interpreting high resolution visible spectra from the ELODIE instrument. For the stellar structure, the ω-model (i.e. the GD model from Espinosa Lara & Rieutord 2011 applied to a Roche model) first allowed us to constrain Altair’s equatorial radius, position angle, rotation velocity and inclination with interferometry, while the spectrum provided a preliminary value for the metallicity of the atmosphere. Full 2D rotating stellar models from the ESTER code were subsequently used to determine the star’s mass, metallicity (treated as a separate parameter from the atmospheric metallicity), and hydrogen content (both core and envelope). The correlations between these four parameters prevented convergence towards a unique solution, and the fitting of Altair’s SED hardly helped. The analysis of Altair’s pulsations, however, provided a solution. The observed frequencies clearly point to the upper end of the mass range and more specifically to M = 1.86 M_⊙. However, due to the correlations between M, Z, X, and X_c, such a mass would lead to X <  X_c if X = 0.700, thus pointing towards higher values of X, such as 0.739 (i.e. the solar hydrogen content based on Asplund et al. (2005), as might be expected for stars in the solar neighbourhood). Even for such a value of X, the value of X_c is only 4 % below that of X, thus indicating that Altair is young. Even then, the solution is not unique due to the correlations between M, Z, X and X_c thus pointing to the need for a full spectroscopic study, based on multiple absorption lines. Such a study may also help resolve the difference found between the atmospheric metallicity (based on a single line) and the bulk metallicity.

This work is the first to combine such a diverse set of constraints in modelling a rapidly rotating star. It highlights the importance of using sophisticated 2D stellar models in interpreting interferometric data and is one of the very few studies which provides a plausible mode identification for acoustic pulsations in a rapidly rotating star. Accordingly, it is also an important step in validating ESTER models from an observational point of view. As such, it paves the way for future studies of other promising targets in a part of the HR diagram which up to now has proven challenging.

Acknowledgments

This research made use of the Jean-Marie Mariotti Center (JMMC) service OiDB (http://oidb.jmmc.fr). This research made use of the SIMBAD and VIZIER databases (CDS, Strasbourg, France) and NASA’s Astrophysics Data System. KB, AD, MR, and DRR acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ESRR (ANR-16-CE31-0007-01), which made this work possible. MR and DRR also acknowledge the International Space Science Institute (ISSI) for supporting the SoFAR international team (http://www.issi.unibe.ch/teams/sofar/). We would also like to thank the referee for his useful insight and comments.

References

Appendix A: Radius of ESTER models

The pressure variation between the polar isobar and the photosphere (defined as where optical depth is unity) at the equator is

(A.1)

assuming that κ_pole ≃ κ_eq. As Δg = Ω²R, we have

(A.2)

For a Roche model, this means

(A.3)

where ε is the flattening coefficient. Then,

(A.4)

Here ℓ is the mean free path of photons. For a 2 M_⊙ model, ℓ = 1/(ρκ)∼2.2 × 10⁹ cm, and R ∼ 1.1 × 10¹¹ cm. If the flattening is maximal (ε = 1/3), ΔR/R ∼ 0.02, if ε = 0.2, ΔR/R ∼ 0.01 and if ε = 0.1, then ΔR/R ∼ 0.004.

Appendix B: Pulsation modes of best-fitting model

Figure B.1 shows the meridional cross-sections of six island modes from the best-fitting model. A pseudo-logarithmic colour scale⁸ is used to bring out faint details. As can be seen, all of these modes are mixed with gravity modes, or even a rosette mode for the first one.

Fig. B.1. Meridional cross-sections of island (or mixed gravito-island) modes in best-fitting model (M = 1.863 M⊙, X = 0.739). The Lagrangian pressure perturbation, normalised by the square root of equilibrium pressure, is shown. A pseudo-logarithmic colour scale is used to bring out faint details.

All Tables

All Figures

	Fig. 1. Interferometric observables (squared visibilities and closure phases) versus theoretical predictions from the best model (see Sect. 7). Each colour indicates a different baseline, with the corresponding telescopes shown in legend.
In the text

	Fig. 2. uv coverage of PIONIER and GRAVITY observations of Altair.
In the text

	Fig. 3. Observed and theoretical MgII line at 4481 Å of the best ESTER model (see Sect. 7. The parameters of the corresponding model are shown in Table 5).
In the text

Fig. 4. Side view of a flattened star. The primed coordinates are the Cartesian coordinates used to describe the plane of the sky, as seen from Earth. (Y′, Z′) are the axes of the sky plane, and X′ the line of sight. The non-primed coordinates are linked to the star, and can be obtained by rotating the primed coordinates around the Y′ axis so that the rotation axis Z is at an inclination angle i from the line of sight. α is the angle between the normal n to the surface at an arbitrary point on the stellar surface and the line of sight. We recall that μ = cos α (cf. Eq. (8)). The star is shaded based on a 3D rendering of its shape, and not on physical effects such as gravity or limb-darkening. Other figures such as Figs. 12 and 13 show such effects.

In the text

	Fig. 5. Corner plot showing the convergence of the MCMC model-fitting for Req, Ω, i, and PA, using interferometric observations (PIONIER + GRAVITY) and the ω-model. The run was made with 200 steps and 100 walkers. Only the last 50% steps are displayed here as convergence was met within the first 100 steps (“burn-in” phase). These results are listed in Table 5.
In the text

	Fig. 6. Positions in the (i, Ω) plane of the χ2 minima when fitting interferometric data with Ω and i as free parameters, using M = 1.80 M⊙, and several values of [M/H] for the atmosphere models. [M/H] goes from −0.5 to 0.9 with a 0.1 step. The grey rectangle shows the area enclosed within the error on i and Ω shown in Fig. 5.
In the text

Fig. 7. χ2 maps of the fitting of ESTER models to both interferometric and spectroscopic data. The upper and lower rows correspond to X = 0.739 and 0.700, respectively. For each mass and X value, only Z and Xc vary, while every other parameter is fixed to the values previously obtained (see Fig. 5). The decimal logarithm of the reduced χ2 is shown in colour. The red dot in each subpanel corresponds to the minimum χ2 value.

In the text

Fig. 8. χ2 cuts at fixed Z values (upper row), and along the valley with low χ2 values (lower row), for X = 0.739. The points in the lower row correspond to the minima in the upper row. As can be seen, a definite minimum appears even within the valley with low χ2 values. Z values are shown in legend, and only the values for which the optimal Xc/X is in the range [0.5−1.0] are shown in the lower row.

In the text

	Fig. 9. Spectral energy distribution produced in the direction of the observer (the inclination of the star being the one shown in Table 5), for different ESTER models with X = 0.700 and 0.739. is the mean surface temperature of the models. The blue dots correspond to the observed SED data from the Vizier website.
In the text

Fig. 10. Observed and theoretical pulsations for Altair. The vertical grey lines that span the plots are the observed pulsations (modes f1 to f7). Their thicknesses are proportional to the observed amplitudes. The blue vertical segments are theoretical pulsations for six different ESTER models which satisfy Eqs. (14) and (15). The theoretical modes are successive m = 0, (or equivalently ℓ = 0 and 1) modes apart from the first mode (see text for details). The lengths of these segments are proportional to the disk-integrated mode visibilities. Left plot: theoretical frequencies whereas the right plot shows the same set of frequencies after applying a suitable homology scaling (e.g. Kippenhahn & Weigert 1990), thus leading to the scaled masses indicated on the figure (assuming the equatorial radius is fixed).

In the text

Fig. 11. Theoretical pulsation spectrum for the M = 1.863 M⊙, X = 0.739 model. The dark grey lines that vertically span the plot are the observed pulsations (modes f1 to f7). Their thicknesses are proportional to the observed amplitudes. The light grey lines of different heights are the theoretical pulsations computed from the model. The red line segments correspond to identified island (or mixed gravito-island) modes, assumed to correspond to the observed modes. Meridional cross-sections of these 6 modes, labelled with the letters “a” to “f”, are shown in Fig. B.1. We note that no island pulsation mode was clearly identified around 18 c/d, probably as a result of an avoided crossing.

In the text

	Fig. 12. Surface map of the effective temperature of our best ESTER model (parameters in Table 5). The dashed line marks the equator. The values are in Kelvin.
In the text

	Fig. 13. Monochromatic intensity map of our best ESTER model (parameters in Table 5), at 1.5 μm in the H band. The values are in erg s−1 cm−2 cm−1 srad−1.
In the text

	Fig. 14. Meridional cut of our best ESTER model (parameters in Table 5). The colours represent the angular rotation rate.
In the text

	Fig. 15. Surface map of the rotation rate of our best ESTER model (parameters in Table 5). The colour represents the rotation period (in hours).
In the text

	Fig. 16. Surface rotation rate as a function of the co-latitude, for our best ESTER model (parameters in Table 5).
In the text

	Fig. 17. Surface map of the velocity of our best ESTER model (parameters in Table 5) projected onto the line of sight. The values are in km s−1.
In the text

	Fig. 18. Squared Brunt-Väisälä frequency N2 of the final model, as a function of the radius r and colatitude θ. The thermally unstable layers (regions where N2 <  0), are colour-coded. These regions are bounded by solid lines, with a dashed line showing the minimum of N2 in this particular region.
In the text

	Fig. 19. Curvature radii along the latitudinal () and azimuthal () directions, and the actual radius of the model star r(θ). The vertical line marks the equator.
In the text

	Fig. B.1. Meridional cross-sections of island (or mixed gravito-island) modes in best-fitting model (M = 1.863 M⊙, X = 0.739). The Lagrangian pressure perturbation, normalised by the square root of equilibrium pressure, is shown. A pseudo-logarithmic colour scale is used to bring out faint details.
In the text