© The Authors 2024

1 Introduction

Spectral synthesis is one of the most powerful techniques to characterise a star. Comparing the high-resolution spectra of a given target with synthetic models usually provides very accurate stellar parameters. The spectroscopic analysis must be complemented with a detailed analysis of the spectral energy distribution, built from photometric observations, and, when feasible, with the use of astrometric and interferometric observations.

A vast amount of work has been done in the field of spectral synthesis: an extensive list of the main 1D-LTE codes available can be found in the introduction of the paper by Wheeler et al. (2023). All these codes allow us to compute spectra for a given set of parameters; in particular, single values of the effective temperature, T_eff, and gravity, log g (other inputs, such as metallicity, [M/H], and microturbulence are also required). In some cases – for example, SYNTHE (Kurucz 2014), the codes can be adapted to simulate the spectrum of a rotating star by computing individual surface intensities at different inclinations through the atmosphere, applying the appropriate Doppler shifts, corresponding to the projected rotation speed, υ sin i, to the emergent spectra.

Single values of T_eff and log g in modelling a stellar spectrum imply the underlying limitation of spherical symmetry. This approximation breaks down for rapidly rotating stars¹: in that regime, the star becomes oblate, and all the relevant photospheric variables, in particular the radius, temperature and gravity, become functions of the latitude, making the problem complex both from the theoretical and computational points of view. Examples of the departure from spherical symmetry are the results of the works – all based on interferometric observations – by Bouchaud et al. (2020) on Altair (α aql, A7V), who find that R_eq/R_pole ≃ 1.282; Domiciano de Souza et al. (2014) on Achernar (α Eri, B6Vpe), giving a ratio R_eq/R_pole ≃ 1.352; and Monnier et al. (2012), on Vega (α Lyr, A0V), for which R_eq/R_pole ≃ 1.13 (Model 3 of that paper).

In consequence, the first issue to be tackled before proceeding to the computation of a synthetic spectrum for a rotating star is that of its structure; in particular, finding out the dependence of R, T_eff, and log g_eff with latitude; the effective gravity, g_eff, is defined as the vector sum of the classical gravity and the centrifugal acceleration (see Eq. (4) in Sec. 3). This area of research has a long history, whose starting point can be set in the pioneering works by von Zeipel (1924a,b), who found that in barotropic stars the energy flux is proportional to the local effective gravity, leading to $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\propto g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\beta }$, with β = 0.25. This is the well-known so-called von Zeipel law, which was modified by Lucy (1967), proposing a smooth dependence, β = 0.08, for stars with a convective envelope. β is traditionally called the ‘gravity-darkening exponent’, the term ‘gravity darkening’ encompassing all the phenomena involved when the rotation of the star is considered – polar temperature and gravity larger than the equatorial values, polar radius shorter than the equatorial radius – is the common terminology today. Espinosa Lara & Rieutord (2012) discussed the relevance of gravity darkening and warned about the caveats posed by the dependence of the resulting laws on the stellar atmosphere models chosen.

The review by Rieutord (2006) gives a summary of the advances in modelling rapidly rotating stars in the decades preceding that paper. In recent years, substantial progress has been made. In particular, concerning the work in this paper, we mention ESTER (Evolution STEllaire en Rotation, Espinosa Lara & Rieutord 2013; Rieutord et al. 2016, and references therein). ESTER is the first code computing, in a consistent way, 2D models of fast-rotating stars, including their large-scale flows. A semi-analytical approximation of this code is used in this work.

The main goal of this paper is to provide a set of methods, described in as much detail as possible, to carry out from scratch the structure and computation of the synthetic spectrum for a rotating star. To our knowledge, there is no publicly available code to carry out these combined tasks. In particular, the prescription presented here for computing the stellar structure has the advantage of being valid for any rotation and is not restricted to slow rotators, in contrast with the von Zeipel approximation (see Sec. 3). A good example of the utility of these tools is the work by Lazzarotto et al. (2023), in which the authors combine the use of synthetic spectra and the ESTER model to carry out a photometric determination of the inclination, rotation rate, and mass of rapidly rotating intermediate-mass stars.

The paper is organised as follows: in Sec. 2, we describe how a star with an inclination angle, i, with respect to the line of sight is seen by the observer, and how to carry out the projection onto a 2D plane. In Sec. 3, we describe how to obtain the relevant parameters for a rotating star required to carry out the spectral synthesis. In Sec. 4, we describe how to build the synthetic spectrum of a star where R, T_eff, and log g_eff are functions of the latitude. In Sec. 5, we apply the whole formalism to Vega as a testbed. Sections 6 and 7 include a discussion of the results and some conclusions. Since this paper deals with formalisms of different areas of stellar physics – namely, geometry, structure, spectral synthesis, and limb darkening – we give in each section the basic information and equations, and direct the reader to the appropriate references.

Fig. 1 Geometry of the problem. Left: an oblate star, rotating with an angular speed, Ω, and seen equator-on by the observer. Right: the star is inclined by an angle, i, with respect to the line of sight. Polar coordinates (r, θ, ϕ) are used, where θ is the colatitude (see Appendix A for details).

2 The geometry

Figure 1 shows the geometry that is used throughout the paper. Initially, all the calculations are done considering that the rotation axis is perpendicular to the line of sight, which coincides with the z axis (left). The star is then inclined by an angle of α = π/2 − i around the x axis, where i is the inclination (i = 0, pole-on, i = π/2, equator-on) (right). Polar coordinates (r, θ, ϕ) are used, where θ is the colatitude (0 for latitude π/2, π for latitude − π/2), and ϕ the azimuthal angle. Since the stellar rotation takes place around the y axis, all variables are only functions of colatitude, and are symmetric with respect to the equator.

In Appendix A, all the details concerning the projection of the 3D star onto the 2D plane of the sky, and how some quantities are seen from the point of view of the observer, are given. Information on how to compute all the relevant geometrical variables that arise when dealing with an oblate star is also provided.

3 Parameterisation of the stellar surface

In this work, we follow the model of Espinosa Lara & Rieutord (2011; ER11, hereafter), also called ‘ω-model’. A very detailed discussion of its derivation can be found in the work of Rieutord (2016). Their starting point is the fact that the gravity darkening of rapidly rotating stars is not well described by the von Zeipel (1924a) law, parameterised, as we mentioned before, as $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\propto g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\beta }$, where T_eff and g_eff are the effective temperature and effective gravity, respectively. Their work was triggered by the fact that some interferometric works (see references in ER11) showed that von Zeipel’s approach seems to overestimate the temperature difference between the pole and the equator of the star. The formalism presented in ER11 allows the computation of R, T_eff, and log g_eff of the photosphere of a rotating star, improving the results obtained following von Zeipel’s prescription. Although simple, mainly from a computational point of view, the ER11 model is able to reproduce the results of more complex models; in particular, the above-mentioned ESTER, as can be seen in Fig. 2 of ER11. The model is tailored for stars with radiative envelopes.

The two basic equations of the model are $$\frac{1}{{\omega }^{2}r}+\frac{1}{2}{r}^2{\sin }^2\theta =\frac{1}{{\omega }^2}+\frac{1}{2},$$(1)

that is, the Roche model (see ER11), and $$\cos \vartheta +\ln \tan \frac{\vartheta }{2}=\frac{1}{3}{\omega }^2{r}^3{\cos }^3\theta +\cos \theta +\ln \tan \frac{\theta }{2},$$(2)

where $$\omega =\mathrm{\Omega }\sqrt{\frac{R_{\mathrm{e}\mathrm{q}}^3}{GM}}=\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\mathrm{K}}}$$(3)

is the non-dimensional rotation rate, Ω, Ω_K, and R_eq are the angular velocity, the Keplerian angular velocity, and the radius, respectively, the last two at the equator, r = R/R_eq, is the non-dimensional radial coordinate, and θ, as we mentioned, is the colatitude. ϑ is an auxiliar angular variable (see ER11 for details, also concerning the two singularities in Eq. (2) at θ = 0 and π/2).

Equation (1) provides the values of the photospheric radius, r, as a function of θ; then for each colatitude, r is introduced into Eq. (2) to obtain ϑ. Both equations can be solved by bisection, or by a Newton-Raphson method. Once these two variables, r(θ) and ϑ(θ) are computed, the effective gravity and the effective temperature can be obtained from the following expressions: $${\mathit{g}}_{\mathrm{e}\mathrm{f}\mathrm{f}}=(-\frac{GM}{r^2}+{\mathrm{\Omega }}^{2}r{\sin }^2\theta ){\mathit{u}}_r+({\mathrm{\Omega }}^{2}r\sin \theta \cos \theta ){\mathit{u}}_{\theta }$$(4) $$\begin{array}{rl}{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}& ={\left(\frac{L}{4\pi \sigma GM}\right)}^{1/4}\sqrt{\frac{\tan \vartheta }{\tan \theta }}{g}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{1/4}\\ & ={\left(\frac{L}{4\pi \sigma R_{\mathrm{e}\mathrm{q}}^2}\right)}^{1/4}{(\frac{1}{r^4}+{\omega }^4{r}^2{\sin }^2\theta -\frac{2{\omega }^2{\sin }^2\theta }{r})}^{1/8}\sqrt{\frac{\tan \vartheta }{\tan \theta }}.\end{array}$$(5)

We note that Eqs. (4) and (5) include three quantities, namely, the stellar mass and luminosity, and the equatorial radius, which – in particular M and R_eq – are not usually known with an acceptable degree of accuracy. Even the luminosity, L, is a more subtle parameter to estimate in the case of very flattened stars since the same object, seen pole-on or equator-on, would show to the observer different spectral energy distributions, which would lead to different apparent effective temperatures, and hence luminosities, the reason being that the expression $L=4\pi \sigma R^2{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}^4$ loses its meaning since both temperature and radius are functions of the latitude.

In a practical case, when attempting to model a stellar spectrum by building a grid of models, a reasonable range of masses, consistent with the estimated spectral type of the object, can be used in Eq. (4). As for the luminosity and equatorial radius, the first parenthesis of the second expression of Eq. (5), including the exponent 1/4, is basically the equatorial effective temperature, and therefore can be substituted by an estimation of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{e}\mathrm{q}}$, or alternatively of the polar temperature, $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{e}}$, using Eq. (32) of ER11: $$\frac{T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{e}\mathrm{q}}}{T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\text{pole }}}=\sqrt{\frac{2}{2+{\omega }^2}}{(1-{\omega }^2)}^{1/12}\exp (-\frac{4}{3}\frac{{\omega }^2}{{(2+{\omega }^2)}^3}),$$(6)

building the grid using as inputs a range of values of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{e}\mathrm{q}}$ or $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{e}}$ that are also consistent with an initial estimate of the spectral type, and iterating until both the interferometric results are reproduced, and/or the synthetic spectrum matches the observed one. In the most common situation, in which interferometric observations are not available, the peculiar shape of some spectral lines (see Sec. 5) can give hints about the inclination of the star, and together with an estimate of the projected rotation speed, υ sin i, iterate using a range of temperatures, until an agreement between the observations and the model is reached.

We note that, according to the first expression in Eq. (5), von Zeipel’s law, $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\propto g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{1/4}$, is recovered for slow rotations since as ω decreases, ϑ/θ → 1 (see Eq. (2)). As a final remark, we also note that in previous works modelling interferometric data (Aufdenberg et al. 2006; Monnier et al. 2012) the methods involving the calculation of the stellar radius, temperature, and gravity make explicit use of the gravity-darkening law $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\propto g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\beta }$, whereas the ER11 formalism used in this work allows us to avoid the whole discussion of what the appropriate value of the exponent, β, is.

4 Spectral synthesis

The spectral synthesis of a star whose relevant photospheric variables are functions of the latitude is, from a computational point of view, substantially more difficult than the classical single-temperature, single-gravity modelling; however, it is conceptually fairly intuitive and can be carried out by following these steps:

1. Once the structure is computed according to the prescription described in Sec. 3, the star is divided into cells delimited by the intersection of a mesh of parallels and meridians with separations of Δθ = Δϕ = 1°; that implies 180 × 360 = 64 800 cells, of which half are visible to the observer. We point out that this discretisation of the stellar surface is uneven, in the sense that the areas of cells near the polar regions are smaller than those of cells near the equator. A discretisation keeping the surface area of all cells constant (see e.g. Bouchaud et al. 2020; Lazzarotto et al. 2023) leads to exactly the same results as the ones presented in this work. A finer mesh does not result in any improvement or refinement of the output spectrum. No numerical noise appears in the results from any of the discretisations.

2. Since all variables, in particular T_eff and g_eff, are only functions of latitude, 90 synthetic spectra, corresponding to the cells with colatitudes between 0 and π/2, are computed for the corresponding values of temperature and gravity. The angular speed, metallicity, and microturbulence are fixed. These synthetic spectra, which contain the fluxes in erg cm⁻² s⁻¹ Å⁻¹, are not rotationally broadened.

3. After the star is rotated by an angle, α = π/2 − i, around the x-axis, each individual cell is seen by the observer with a projected area, (ΔA)_p, a radial velocity, $v_z^{\mathrm{\prime }}$, and an angle, γ, between the line of sight and the normal to the surface element, the latter being relevant for the correction for limb darkening, C_ld(λ), to be applied (see Appendix B for details of the computation of the limb-darkening coefficients). Taking into account all these factors, the total flux at a given wavelength, λ_j, is $$F\left({\lambda }_j\right)={\displaystyle \sum _{i=1}^{N_{\text{cells }}}}{F}_i({\lambda }_j,v_{i,\mathrm{z}}^{\mathrm{\prime }})(\mathrm{\Delta }A{)}_{i,\mathrm{p}}{C}_{\mathrm{l}\mathrm{d}}({\gamma }_i,{\lambda }_j),$$(7)

   where $F_i({\lambda }_j,v_{\mathrm{z}}^{\mathrm{\prime }})$ is the flux at λ_j of the synthetic spectrum computed for the particular values of T_i and log g_i of that cell, redshifted or blueshifted, according to the value of the radial velocity, $v_{z,i}^{\mathrm{\prime }}$, of the cell (see Eqs. (A.3)–(A.5) and (B.3)).

The individual synthetic spectra are computed using the codes ATLAS and SYNTHE (Kurucz 2014) and the models containing the elemental abundances and the stratification of the stellar atmospheres as a function of temperature, gravity, metallicity, and microturbulence velocity (Castelli & Kurucz 2003). The ATLAS code allows us to compute a model atmosphere for any value of temperature and gravity from a close model already computed in the Castelli & Kurucz grids. The spectral synthesis is carried out using SYNTHE, with a resolution of λ/Δλ = 100 000 at 450 nm (0.0045 nm pixel⁻¹). The GNU-linux version of the codes by Sbordone (2005) is used².

A grid of 3668 synthetic models with T between 7000 and 20 000 K (step 100 K), log g = 3.0, 3.5, 4.0, 4.5, and metallicities of [M/H] = −2.5, −2.0, −1.5, −1.0, −0.5, 0.0, and +0.5 was computed beforehand; the microturbulent velocity is 2 km s⁻¹. For all the cells at a colatitude, θ_i, with parameters (T_i, log g_i), the corresponding synthetic spectrum is computed by linear interpolation between the four closest neighbouring models in the grid, those bracketing at a time the temperature and the gravity of the cell; that is, the four models in the grid – for a given metallicity − (T_j, log g_k), (T_j+1, log g_k), (T_j, log g_k+1) and (T_j+1, log g_k+1), have to fulfil T_j < T_i ≤ T_j+1, log g_j < log g_i ≤ log g_j+1.

The interpolation is easily carried out in this way: first, two constants are defined as $$C_{\mathrm{T}}=\frac{T_{j+1}-T_i}{T_{j+1}-T_j}{C}_{\mathrm{g}}=\frac{\log g_{j+1}-\log g_i}{\log g_{j+1}-\log g_j},$$(8)

then two intermediate models are computed: $$\begin{array}{rl}& (T_i,\log g_j)=C_{\mathrm{T}}\cdot (T_j,\log g_j)+(1-C_{\mathrm{T}})\cdot (T_{j+1},\log g_j)\\ & (T_i,\log g_{j+1})=C_{\mathrm{T}}\cdot (T_j,\log g_{j+1})+(1-C_{\mathrm{T}})\cdot (T_{j+1},\log g_{j+1})\end{array}$$

which are combined to give the final one (T_i, log g_i) for the ith cell: $$(T_i,\log g_i)=C_{\mathrm{g}}\cdot (T_i,\log g_j)+(1-C_{\mathrm{g}})\cdot (T_i,\log g_{j+1}).$$(9)

5 Vega as a testbed

Vega (α Lyr, HD 172167, HIP 91262, HR 7001) is one of the most extensively studied stars. It is well known that it is used as standard for the calibration of several photometric systems (Bessell 2005) and that it is surrounded by a debris disc, discovered by Aumann et al. (1984), which triggered an intensive study at infrared wavelengths (Sibthorpe et al. 2010, and references therein). However, this object turned out not to be the perfect standard, showing anomalies in its luminosity (Petrie 1964; Millward & Walker 1985), some peculiarly shaped absorption lines (Gulliver et al. 1991), and its radius (Hanbury Brown et al. 1967; Ciardi et al. 2001) in comparison with other A0 V stars.

Concerning this paper, our interest focuses on the fact that all these anomalies are now explained by the fact that Vega is a rapidly rotating star being seen almost pole-on; that is, with a small inclination angle, as has been shown by a number of interferometric studies (Aufdenberg et al. 2006; Peterson et al. 2006; Monnier et al. 2012) and spectroscopic analyses (Gulliver et al. 1991; Hill et al. 2004, 2010; Takeda et al. 2008; Yoon et al. 2010; Takeda 2021). Table 1 in Takeda (2021) gives a summary of the values of the projected equatorial velocity, υ_eq sin i, the equatorial velocity, υ_eq, the inclination angle, i, the polar and equatorial radii, R_pole and R_eq, and the rotation period, P, according to different analyses.

To check the reliability of the methods described in this paper, we use some of the results of the above-mentioned works to check whether the structural model (Sec. 5.1), and then the spectral synthesis model (Sec. 5.2), can reproduce the observed properties.

5.1 Structure and stellar parameters

Table 1 shows in the upper part the input parameters of the model; the inclination, stellar mass, and equatorial radius have been taken from Monnier et al. (2012) (their Model 3, the ‘concordance model’). Since the model also requires the polar temperature and ω as inputs, T_pole was explored within the uncertainty interval given by Monnier et al. to match the luminosity, and ω was fixed to match the value of υ sin i derived by Takeda (2008). The lower part of the table (Col. 2) shows the results of the formalisms described in Sects. 2 and 3; Col. 3 shows, for comparison, some parameters derived from interferometric (Model 3 by Monnier et al. 2012) and spectroscopic analyses (Takeda 2021).

The results derived in this work are in general consistent with those from previous modellings, although some discrepancies are apparent: the temperature drop from pole to equator, 1098 K in our case, is in agreement with that by Monnier et al. (2012) (1160 K); both are substantially smaller than those by Peterson et al. (2006; > 2400 K), and Aufdenberg et al. (2006; 2250 K).

Some details about the calculations: the stellar luminosity was computed by adding for all cells the quantity $\sigma (\mathrm{\Delta }A{)}_i{T}_{\mathrm{e}\mathrm{f}\mathrm{f},i}^4$, where (ΔA)_i and T_eff,i are the surface area and the effective temperature of the ith cell; and the average effective temperature was estimated from the expression $L=4\pi R_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}^2\sigma T_{\mathrm{e}\mathrm{f}\mathrm{f},\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}^4$, where R_aver is an average of the radius, computed in the interval of colatitudes, [0, π/2].

Figure 2 shows the stellar radius normalised to R_eq (black) and the temperature (red), plotted against the colatitude for the northern hemisphere (the results for the southern hemisphere are symmetrical). Figure 3 shows two colour plots showing the temperature and effective gravity profiles for Vega according to the results of our modelling.

Fig. 2 Results for Vega from the model computed using the input parameters shown in Table 1 and the formalism described in Sec. 3. Radius, normalised to the equatorial radius (black), and temperature (red) are plotted against the colatitude, θ.

5.2 Spectral synthesis

In this section, we check whether the results derived in Sec. 5.1 together with the formalism described in Sec. 4 are able to reproduce the peculiar shape of some features of the spectrum of Vega. The high-resolution, high-signal-to-noise spectrum atlas of Vega from Takeda et al. (2007) has been used throughout. We do not intend here to make an analysis of the elemental abundances, which has already been carried out by other authors (see e.g. Ilijic et al. 1998; Qiu et al. 2001; Takeda 2008; Takeda et al. 2008), but rather to make sure that the whole set of procedures described in the previous sections is able to reproduce the peculiar shapes of the absorption lines of the spectrum of Vega.

Figure 4 shows the profiles of 30 lines, from different species, both neutral and ionised. The observed profiles and the results of our modelling are plotted in black and red, respectively. In cyan, the profiles resulting from a single-temperature, single-gravity synthetic spectrum computed with the average parameters, $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}$ and log $g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\text{aver }}$, listed in Table 1, are also plotted. For the sake of a better graphical display of the whole set of lines, both the observed and the synthetic profiles have been re-scaled and normalised, placing the continuum at intensity 1.0 and the bottom of the profiles at intensity 0.9, whereas the single-temperature, single-gravity profiles have been scaled so that they fit the wings of the observed absorptions.

The first five profiles of the left panel, from top to bottom, show rounded shapes, typical of lines broadened with a classical rotation profile; however, the remaining 25 profiles have a variety of shapes, the most extreme cases being those that are almost rectangular, such as Fe I 402.19, 406.80, 449.45 nm and Ba II 455.40, 493.41 nm, or those showing two deeper components at velocities close to ±υ_eq sin i, like Ca I 445.48 nm. It is apparent that in all cases the agreement between the shape of the observed profiles and the results from the formalism described in Sec. 4 is remarkable. To our knowledge, the only work that also accurately reproduces the peculiar profiles of the Vega spectrum is that from Takeda et al. (2008).

In order to understand the peculiar profile of some lines, we show as an example an analysis of two nearby lines; namely, Fe II 445.16 nm, which has a rounded shape, and the above-mentioned Ca I 445.48 nm. The left panel of Fig. 5 shows the colour plot of the radial velocity of each surface element for the Vega model. It is well known that the loci of equal radial velocities, in the case of solid-body rotation, as seen by an observer, are lines parallel to the rotation axis projection; that is, all points in the disc with x=constant (see e.g. Gray 1992, chapter 18) have the same value of radial velocity. Since the star is seen almost pole-on, the regions close to the borders of the disc – the limb, which coincides with the equator – and in particular those with the highest projected rotation speeds, are the ones with the lowest temperatures and gravities. Therefore, the ionisation balance of some species differs from regions with low temperatures and gravities to regions near the pole (see Fig. 3). The vertical purple lines superimposed on the colour plot of the star in the left panel of Fig. 5 delimit 16 strips. Their widths have been computed to fulfil the condition that each one contributes to the full synthetic spectrum with the same amount of flux in the continuum near the lines.

The right panel of Fig. 5 shows the stellar lines (black), with the corresponding model profiles (red) framed in green (Fe II) and blue (Ca I) boxes, respectively, and the equivalent widths (EWs) of the contributions to the profiles from each one of the 16 strips, using the same colour code. The different ranges of EWs for both lines are apparent: whereas the values for the Fe II values are more even, leading to a rounded profile, the outer strips dominate the absorption of the Ca I line, producing its peculiar shape.

6 Discussion

6.1 The ω-model versus the von Zeipel approach

Despite the fact that the ω-model, and the corresponding spectral synthesis, work well for Vega, it is interesting to point out that Takeda (2008), using the Roche model and the von Zeipel value of the gravity darkening exponent, β = 0.25, also found a good agreement between modelling and observations. Monnier et al. (2012) showed that the value of β that best fit the observations was 0.231 ±0.028, in agreement with the von Zeipel value; therefore, both the ω-model and the von Zeipel approximation give accurate results for Vega.

To show the real potential of the ω-model, it is interesting to explore the case of a much faster rotator for which the β exponent is much less than 0.25. The case of Achernar is a good one to check, since a value of β = 0.166 must be used to reproduce the observed results (Domiciano de Souza et al. 2014). The ω-model, using as inputs – all extracted from Domiciano de Souza et al.’s paper – M/M_⊙ = 6.1, T_pole = 17124 K, R_eq/R_⊙ = 9.17, and ω = 0.838, gives T_eq = 12700 K, in excellent agreement with the best fit of the CHARRON RVZ model to the VLTI/PIONIER H band observations, which gives T_eq = 12673 K. In contrast, the von Zeipel model (β = 0.25) gives T_eq = 10 880 K, almost 1800 K off the value derived from observations, which is nicely reproduced by the ω-model (see Domiciano de Souza et al. 2014, for further details). Obviously, that deviation would have a large impact on reproducing the observed spectrum, via spectral synthesis.

Fig. 3 2D structure of the temperature and effective gravity of Vega as seen by an observer. Parallels and meridians separated by 10 degrees are also shown. The values in the axes are distances scaled to Vega’s equatorial radius.

6.2 Contexts in which this work can be useful

Once the whole formalism of the ω-model plus the spectral synthesis have been put together and successfully tested with the paradigmatic case of Vega, and the reassuring case of Achernar mentioned in the previous subsection, it is interesting to point out explicitly in which contexts all this can be useful. With this purpose, we have computed several models whose details are given in Table 2. The models share as fixed inputs some of the parameters of the structural model of Vega (see Table 1); namely, the stellar mass, polar temperature, equatorial radius, and metallicity. Models in the upper half of Table 2 have been computed with a fixed inclination, i = 6.2°, the value obtained for Vega, and five values of ω, resulting in five values of υ sin i; namely 10, 15, 20, 21.4 (Vega), and 25 km s⁻¹. Models in the lower half of Table 2 have been computed with a fixed value of the projected υ sin i, 21.4 km s⁻¹, which is again the value we obtained for the Vega model, and five values for the inclination; namely, 5, 6.2 (Vega), 10, 15, and 20 degrees.

Figure 6 shows, as an illustrative example, the synthetic spectra of the two sets of models in a short wavelength interval between 445.0 and 446.0 nm. The upper and lower panels in the figure correspond to the models in the upper and lower parts of Table 2. The spectra have been normalised to the intensity at 445.2 nm and contain five lines: Fe II 445.16, Ca I 445.48, Fe II 445.58, Ti II 445.66, and Fe II 445.91 nm. The colour codes of the spectra – red, black, cyan, purple, and orange – correspond to decreasing values of υ sin i (upper panel) and increasing values of inclinations (lower panel), the model for Vega being plotted in black.

What is interesting in this plot is how sensitive the profiles are to changes in inclination and υ sin i, a conclusion that can be extended to the full spectral range. In particular, the profiles of the lines with peculiar shapes, as in the cases of Ca I 445.48, Fe II 445.58, and Fe II 445.91 nm, change very dramatically as the inclination decreases. Very interesting, too, is the comparison between the behaviour of the normal rounded-shape Fe II 445.16 nm line and the peculiar Ca I 445.48 nm profile in the upper panel: whereas the Fe II line behaves as one would expect as the value of υ sin i increases, the shape and depth of the Ca I line changes drastically.

The model plotted as a dotted grey line in the lower panel has been computed for i = 5°, ω = 0.632, assuming the von Zeipel approximation; this model must be compared with the one plotted in red, computed with the same parameters, but under the assumptions of the ω-model. That value of ω would be associated with a β exponent ~0.197, quite far from β = 0.25. As can be inferred from the values of T_eq/T_pole for both models in Table 2, the equator is almost 400 K cooler when the von Zeipel approximation is used, which results in deeper lines, in particular those with peculiar profiles, leading to erroneous determinations of abundances. This is a good example of the influence of the β exponent on the line shapes and intensities.

All this shows the potential of a detailed spectral analysis to find structural and physical parameters and inclinations of this kind of stars. Good examples are the works by Takeda et al. (2008) and Takeda (2020, 2021) disentangling Sirius A and Vega’s properties using spectral line profiles, or Fourier analysis.

A quantitative analysis of the usefulness of the proposed formalisms is relevant. Regarding the inclination of the star, it is apparent that very clear changes are observed in certain line profiles as i moves in the range between 0 and ~ 20 degrees; at larger inclinations, most of the lines are insensitive to this parameter. It is easy to prove that the probability of a star having an inclination between i and i + Δi is P(i, i + Δi) = sin i Δi, and therefore the probability of finding a star, among a large set of objects with an inclination in the interval of [0,20] degrees is ~ 0.0603. As a first impression, one might consider the whole modelling effort to be disproportionate considering that the number is small; however, a query to the Gaia DR3 catalogue asking how many stars with spectral types between A0 and A9 − for which the methods presented here would be useful − with parallaxes, ϖ, with relative errors of Δϖ/ϖ < 0.20, are ~ 146 000 (ϖ ≥ 2 mas) and ~ 231200 (ϖ ≥ 1 mas). In other words, in a sphere with a radius of 1 kpc, we would find around ~ 14 000 stars in that range of spectral types with inclinations less than 20 degrees. The constraint of −0.037 ≤ BP − RP ≤ +0.377 to bracket the interval A0-A9 has been used³.

Fig. 4 Plotted in black, profiles of 30 absorption features of the spectrum of Vega from different species. Superimposed on the observed lines are the results of the spectral synthesis described in Sect. 4, and the spectral single-T, single-log g synthesis carried out with the values of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}$ and log $g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\text{aver }}$, listed in Table 1, plotted in red and light blue, respectively (see the text for further details).

Fig. 5 Left: colour plot of the projected radial velocity of each surface element of the star. The purple lines delimit 16 regions contributing the same amount of flux to the continuum in the region near the Fe II 445.16 and Ca I 445.48 nm lines. Right top: observed (black) and model (red) profiles of the Fe II (green frame) and Ca I (blue frame) lines. Right bottom: equivalent widths of the contributions to the lines from each of the individual 16 strips shown in the left panel.

7 Conclusions

In this paper, we provide a combined method to compute the structure of rapidly rotating stars and build their synthetic spectra. A summary of the main features of the whole formalism follows:

1. The ω-model by Espinosa Lara & Rieutord (2011) has been implemented to compute the relevant parameters of the photosphere of rapidly rotating stars – namely, the radius, R, effective temperature, T_eff, and effective gravity, g_eff – as a function of the colatitude, θ. The method, relatively simple from a computational point of view, is able to reproduce the results of more complex models. One of the big advantages of this formalism is that it avoids the discussion, and hence the subsequent computation, or ad hoc assignment, of the appropriate gravity darkening exponent, β. The model is applicable to stars with radiative envelopes (Sect. 3). In those situations in which some of the approximations inherent in the ω-model – that is, mass concentrated near the centre of the star and rigid rotation (Roche model) – are no longer valid, the original ESTER model should be used.

2. A detailed method of how to compute the synthetic spectrum of a rapidly rotating star, at any inclination angle, i, with respect to the line of sight, is presented. The model makes use of the suite of codes, ATLAS and SYNTHE (Kurucz 2014), and the grid of model atmospheres by Castelli & Kurucz (2003; Sec. 4 and Appendices A and B).

3. The combined methods summarised in items 1 and 2 above were applied to the particular case of Vega, obtaining results regarding both the structure and the synthetic spectrum that are compatible with previous works. The fitting of the spectral lines was remarkable, both for those with normal, rounded shapes and those with peculiar profiles (Sec. 5).

4. In addition, Appendix A describes in detail how to treat, from a strict geometrical point of view, all the relevant variables when a rotating star is seen with a given inclination with respect to the line of sight.

Although this work has focused on the spectral synthesis of rapid rotators, the tools provided in this paper can be useful in other contexts:

1. To locate the position of a star in colour-magnitude diagrams, since a star deformed by rapid rotation appears brighter and hotter when it is observed near pole-on (see e.g. Pérez Hernández et al. 1999; Bastian & de Mink 2009; Girardi et al. 2019).

2. To find and estimate the inclination of the rotation axis with respect to the line of sight in those cases without making use, in the first instance, of interferometric measurements (Takeda et al. 2008). This would be useful to search for potential pre-main-sequence stars of spectral types earlier than F hosting Jupiter-like planets. In particular, there is indirect evidence that Herbig Ae/Be stars with low metallicities could be good candidates to host such giant planets (Kama et al. 2015; Guzmán-Díaz et al. 2023). The method presented here would allow a detailed metallicity analysis, and a subsequent filtering of targets according to their inclination, which is suitable in the case of low-inclination systems of potential interferometric and/or direct imaging studies.

3. The role of the inclination is particularly important in modelling accretion processes for young objects of intermediate mass. In the scenario of magnetospheric accretion, the shape and intensity of the spectral lines are strongly dependent on the assumed inclination (Muzerolle et al. 2004; Mendigutía et al. 2011). Regarding the alternative scenario of boundary layer continuum models, the dependence on the inclination is also critical (see e.g. Fig. 5 in Mendigutía 2020).

Fig. 6 Models of a small section of the optical spectrum showing the variations in the spectral lines with υeq sin i for a fixed inclination (upper panel), and with the inclination, i, for a fixed value of υeq sin i. The parameters of the models are given in Table 2. We note that the profile of the Ca I 445.48 nm line, which appears in Figs. 4 and 5, and that of the Fe II 445.16 nm, which also appears in Fig. 5, are scaled in those figures to an intensity at the bottom of ~0.9 for plotting purposes. In this figure, all three lines show their actual relative intensities.

Acknowledgements

The author is very grateful to the referee, Prof. Michel Rieutord, and his colleagues, Alain Hui-Bon-Hoa and Axel Lazzarotto, for providing very useful comments, suggestions and references that, no doubt about, have improved the contents and scope of the paper. This research has been funded by grants AYA2014-55840-P, PGC2018-101950-B-I00 and PID2021-127289-NB-I00 from MCIN/AEI/10.13039/501100011033/ and FEDER. The author is grateful to Francisco Espinosa-Lara for useful discussions on the ER11 formalism, Antonio Claret for some guiding for the computation of the limb-darkening coefficients, and Almudena Alonso-Herrero, Olga Balsalobre-Ruza, Carlos Eiroa, Jorge Lillo-Box, Ignacio Mendigutía, Enrique Solano and Eva Villaver for their help and comments to several sections of this paper. Special thanks also to Antonio Parras and Sergio Suárez for their work keeping up and running the computing centre.

Appendix A The geometry of the problem

According to the notation in Fig. 1, a point on the stellar surface with coordinates $$\mathit{r}=\{\begin{array}{l}x=r\sin \theta \sin \varphi \\ y=r\cos \theta \\ z=r\sin \theta \cos \varphi \end{array}$$

when the star is seen equator-on, is transformed after a counterclockwise rotation around the x axis by an angle, α, is done, into r′ = (x′, y′, z′) by applying to r the matrix R_x(α); namely, $${\mathit{R}}_x(\alpha )=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & -\sin \alpha \\ 0& \sin \alpha & \cos \alpha \end{array}\right)$$ $${\mathit{r}}^{\mathrm{\prime }}=(x^{\mathrm{\prime }},y^{\mathrm{\prime }},z^{\mathrm{\prime }})={\mathit{R}}_x(\alpha )\mathit{r}.$$(A.1)

The observer would see all the points on the stellar surface fulfilling the easy constraint, z′ = y sin α + z cos α > 0. In the case of an spherical object, the unitary vector, u, attached to any point has the direction of the normal to the surface; however, in an oblate object this is not the case, as can be seen in Fig. A.1. The normal to the surface at a given point with colatitude θ, is inclined at an angle of ξ = π/2 − θ + η with respect to the line of sight before proceeding to apply the rotation by an angle, α. The computation of η is a fairly straightforward geometrical problem, as is illustrated in the figure, in which Δθ has obviously been plotted out of scale: tan η ≃ (r₂ − r₁)/(r Δθ).

Fig. A.1 Sketch showing the geometrical problem involved in the computation of the angles, η and ξ, as functions of the colatitude and the radius of a given surface element.

A differential surface area element at a latitude, θ, can be written as ΔA ≃ (r Δθ/ cos η) · (r sin θ Δϕ), and its associated unitary vector normal to the surface, before the star is rotated, has the following expression: $${\mathit{u}}_A=(\cos \xi \sin \varphi ,\sin \xi ,\cos \xi \cos \varphi ).$$(A.2)

Therefore, after applying the rotation, the projected area as seen by the observer would be ΔA multiplied by the z-component of R_x(α) · u_A; namely, $$(\mathrm{\Delta }A{)}_{\mathrm{p}}=\mathrm{\Delta }A(\sin \alpha \sin \xi +\cos \alpha \cos \xi \cos \varphi ).$$(A.3)

Concerning the rotation speed of each surface element, υ, it has only components x and z, as can be seen in Fig. 1 (left): $$\mathit{v}=\{\begin{array}{l}{v}_x=\mathrm{\Omega }r\sin \theta \cos \varphi \\ v_y=0\\ v_z=-\mathrm{\Omega }r\sin \theta \sin \varphi \end{array}\text{. }$$

After applying the rotation, the component of the velocity in the line of sight is $${\nu }_z^{\mathrm{\prime }}=-\mathrm{\Omega }r\cos \alpha \sin \theta \sin \varphi \text{. }$$(A.4)

Finally, the knowledge of the angle γ, between the normal to a given surface element and the line of sight, in the rotated system, must be known in order to apply the correction for limb darkening to the synthetic spectra arising from that surface element.

Fig. A.2 Sketch showing the the geometry of a surface element, characterised by the unitary vector, ${\mathit{u}}_{\mathrm{A}}^{\mathbf{\prime }}$, whose z component forms an angle γ, with the line of sight. μ = cos γ is required to compute the limb-darkening coefficient (see Appendix B)

In the usual notation for that angle, μ = cos γ, and according to Fig. A.2, it can be written as $\mu =\cos \left(u_{\mathrm{A}z}^{\mathrm{\prime }}/\left|{\mathit{u}}_{\mathrm{A}}^{\mathbf{\prime }}\right|\right)$, and since $\left|{\mathit{u}}_{\mathrm{A}}^{\mathbf{\prime }}\right|$ = 1, the value of cos γ is just the z component of R_x(α) · u_A; namely, $$\mu =\cos \gamma =\sin \alpha \sin \xi +\cos \alpha \cos \xi \cos \varphi .$$(A.5)

Appendix B The limb-darkening coefficients

This appendix shows how the limb-darkening coefficients (LDCs hereafter) and the limb-darkening correction C_ld(λ) (see Eq. (7) are computed. The work and notation by Claret & Bloemen (2011a) are followed in this section. The LCDs a_k, k = 1, 4 are defined in such a way that the most general law is adjusted by the following expression: $$\frac{I(\mu )}{I(1)}=1-{\displaystyle \sum _{k=1}^4}{a}_k(1-{\mu }^{k/2})\text{. }$$(B.1)

I(1) is the intensity at the centre of the disc, and μ = cos γ (see Eq (A.5)), where γ is the angle between the normal to the surface area element and the line of sight. Equation (B.1) can be defined for the intensity in a given passband, although in our case we are interested in a monochromatic estimate of that quantity for each of the wavelengths covered by the synthetic spectra.

Claret & Bloemen (2011a) computed, among others, the LCDs for the photometric Johnson-Cousins UBVRI filters (Table 17 available at Claret & Bloemen 2011b). In order to estimate the limb-darkening correction at each wavelength, λ, we linearly interpolate the LCDs in this way: $$a_k(\lambda )=a_k\left({\lambda }_{\mathrm{F}1}\right)+\left[\frac{a_k\left({\lambda }_{\mathrm{F}2}\right)-a_k\left({\lambda }_{\mathrm{F}1}\right)}{{\lambda }_{\mathrm{F}2}-{\lambda }_{\mathrm{F}1}}\right](\lambda -{\lambda }_{\mathrm{F}1}),$$(B.2)

where F1 and F2 stand for ‘Filter 1’ and ‘Filter 2’ and λ_F1, λ_F2 are the effective wavelengths of the filters adjacent to the wavelength, λ, under consideration; that is, λ_F1 < λ ≤ λ_F2. The wavelengths assigned to the UBVRI filters for our computations are 360, 440, 550, 690, and 950 nm, respectively. Following that notation, the correction for limb-darkening applied to the fluxes at a given wavelength is $$C_{\mathrm{l}\mathrm{d}}(\lambda )=1-{\displaystyle \sum _{k=1}^4}{a}_k(\lambda )(1-{\mu }^{k/2}).$$(B.3)

References

All Tables

All Figures

	Fig. 1 Geometry of the problem. Left: an oblate star, rotating with an angular speed, Ω, and seen equator-on by the observer. Right: the star is inclined by an angle, i, with respect to the line of sight. Polar coordinates (r, θ, ϕ) are used, where θ is the colatitude (see Appendix A for details).
In the text

	Fig. 2 Results for Vega from the model computed using the input parameters shown in Table 1 and the formalism described in Sec. 3. Radius, normalised to the equatorial radius (black), and temperature (red) are plotted against the colatitude, θ.
In the text

	Fig. 3 2D structure of the temperature and effective gravity of Vega as seen by an observer. Parallels and meridians separated by 10 degrees are also shown. The values in the axes are distances scaled to Vega’s equatorial radius.
In the text

Fig. 4 Plotted in black, profiles of 30 absorption features of the spectrum of Vega from different species. Superimposed on the observed lines are the results of the spectral synthesis described in Sect. 4, and the spectral single-T, single-log g synthesis carried out with the values of $T_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}}$ and log $g_{\mathrm{e}\mathrm{f}\mathrm{f}}^{\text{aver }}$, listed in Table 1, plotted in red and light blue, respectively (see the text for further details).

In the text

Fig. 5 Left: colour plot of the projected radial velocity of each surface element of the star. The purple lines delimit 16 regions contributing the same amount of flux to the continuum in the region near the Fe II 445.16 and Ca I 445.48 nm lines. Right top: observed (black) and model (red) profiles of the Fe II (green frame) and Ca I (blue frame) lines. Right bottom: equivalent widths of the contributions to the lines from each of the individual 16 strips shown in the left panel.

In the text

Fig. 6 Models of a small section of the optical spectrum showing the variations in the spectral lines with υeq sin i for a fixed inclination (upper panel), and with the inclination, i, for a fixed value of υeq sin i. The parameters of the models are given in Table 2. We note that the profile of the Ca I 445.48 nm line, which appears in Figs. 4 and 5, and that of the Fe II 445.16 nm, which also appears in Fig. 5, are scaled in those figures to an intensity at the bottom of ~0.9 for plotting purposes. In this figure, all three lines show their actual relative intensities.

In the text

	Fig. A.1 Sketch showing the geometrical problem involved in the computation of the angles, η and ξ, as functions of the colatitude and the radius of a given surface element.
In the text

	Fig. A.2 Sketch showing the the geometry of a surface element, characterised by the unitary vector, ${\mathit{u}}_{\mathrm{A}}^{\mathbf{\prime }}$, whose z component forms an angle γ, with the line of sight. μ = cos γ is required to compute the limb-darkening coefficient (see Appendix B)
In the text