© The Authors 2024

1 Introduction

Spectral synthesis of fast-rotating stars is challenging given the many physical processes occurring within the star, from the stellar core to the photosphere. The life cycle of a massive star is partly determined by physical processes related to its rotation velocity. Convective mixing and mass loss are severely impacted if the star rotates at high rotational rates (Maeder & Meynet 2010). In recent decades, significant effort has been put into the simulation of spectra, synthetic magnitudes, and interferometric profiles of fast-rotating stars, and particularly the inclusion of gravity darkening effects, which are the effects of radiative flux redistribution due to surface gravity changes induced by centrifugal forces (von Zeipel 1924a,b; Collins & Harrington 1966; Lucy 1967; Harrington & Collins 1968; Stoeckley 1968; Townsend et al. 2004; Frémat et al. 2005; Espinosa Lara & Rieutord 2011).

In their pioneering works, von Zeipel (1924a,b) and Chandrasekhar (1933) established the foundations of fast stellar rotation, in particular its connection with the stellar structure and evolution. In its canonical form, the gravity darkening theory assumes the gas pressure to be a function of density, known as the barotropic hypothesis, and the hypothesis of rigid rotation and radiative transport. It was soon realised that these assumptions were conflicting, but the debate persisted until recent interferometric observations (Rieutord 2006).

In the classical picture, the gravity force is related to the gravitational potential gradient, itself ruled by the mass distribution through Poisson’s equation. The angular velocity of a rigid rotator is constant. The balance of the gravity force against the centrifugal one gives the total resultant force acting in every surface element. Poisson’s equation must take all these contributions into account. Assuming the validity of the hydrostatic equation, the gas pressure gradient and the gravitational potential should be parallel. The effective gravity is orthogonal to all isobaric surface elements (Tassoul 2007). Assuming the mean molecular weight to be constant along the equipotentials, all variables in the ideal gas state equation are also constants in these paths, and the temperature gradient can be written in terms of the radiative flux as a function of the equipotentials. The radiative flux divergence balances the nuclear energy generation rate, leading to a problem known as von Zeipel’s paradox. The angular velocity cannot directly determine the nuclear energy generation rate, and so a rotating star cannot be in radiative thermal equilibrium. Due to rotation in regions outside the core, the radiative flux divergence is non-zero and changes its sign, causing some places to cool and others to heat. This leads to slow-motion flows known as the Eddington-Sweet currents, and implies that the classical von Zeipel’s theorem overestimates the temperature contrast between the equator and the stellar poles.

Espinosa Lara & Rieutord (2011) presented a novel approach to the classical von Zeipel’s law. In their new prescription, these latter authors suppose the radiative flux to be anti-parallel to the effective gravity everywhere, and the absence of radiative sources outside the stellar core causes the flux divergence to be null. This results in a smoother latitudinal surface temperature profile compared to the classical approach.

From a computational viewpoint, many attempts have been made in recent years to implement gravity-darkening synthesis codes. Many of these codes were given names, others not. Some older versions simulated gravity darkening through Gaussian profiles instead of using the theoretical spectra of stellar atmospheres. George Collins developed perhaps one of the oldest and most pioneering codes in the mid-1960s (Collins 1965), but did not refer to it by a name. Among the more recent codes, we should mention BRUCE (Townsend 1997), ROTSPEC - written by R. H. Townsend (Vinicius et al. 2007), FASTROT (Frémat et al. 2005), Takeda’s model (Takeda et al. 2008, Takeda et al. 2017), and SUPEROTATE (Dall & Sbordone 2011), among others. The original Espinosa-Lara code does not have a specific name and is integrated within the ESTER project. Also, it is worth mentioning an unnamed program used by Lovekin et al. (2006) and Gillich et al. (2008) to perform the spectral synthesis of fast-rotating stars using the temperature Τ (θ) and surface gravity log g(0) latitudinal profiles as input, which were generated with the ROTORC stellar evolution code (Deupree 1990, 1995).

This work presents and describes our code ZPEKTR, which stands for ‘von ZeiPEl’s code for gravity darKening specTRal synthesis’. Levenhagen et al. (2011) originally developed this code in C ANSI language, and after that, the name ZPEKTR, pronounced ‘spectre’, was chosen to avoid homonyms with other synthesis codes (Levenhagen 2014). Interpolation routines used within the code come from the Python SciPy library (Virtanen et al. 2020), while we take numerical functions from NumPy (Harris et al. 2020). Plot procedures are imported from Mat-plotLib (Hunter 2007). In its current version, the code assumes both the classical von Zeipel and the new Espinosa-Lara model. It performs the spectral synthesis for a fast-rotating star, assuming as input spectra a hybrid grid of synthetic LTE/non-LTE models calculated with TLUSTY/SYNSPEC codes (Hubeny 1988; Hubeny & Lanz 1995). For cooler models below 10 000 K, we assume input models by Kurucz (1979) and synthesised with SYNSPEC (Hubeny & Lanz 1995). Additionally, the most recent versions of ZPEKTR are written in Python (Levenhagen et al. 2021; Solar et al. 2022) and also in Julia; the latest version of Julia generally works faster than common programming languages, mostly due to optimisations of array broadcast fusion instead of classical nested loops.

2 Classical von Zeipel model

Assuming the star to be a rigid rotator with an inner radiative transport mechanism, the norm vector of the radiative flux scales to the temperature gradient. This norm should also scale directly with the pressure gradient in order for the hydrostatic equilibrium to keep holding. As stated by the Stefan-Boltzmann law, the temperature field should be proportional to the radiative flux, which implies that the temperature be proportional to the local gravity: $$T_{\text{eff}}(\text{Ω},\theta )=\text{Γ}{\left(g_{\text{eff}}(\text{Ω},\theta )\right)}^{\beta },$$(1)

where Γ is the Eddington factor, which is dependent on the stellar luminosity, mass, and rotation rate (Maeder & Stahler 2009), and β = 0.25 represents the canonical gravity-darkening exponent for stars with radiative transport. The β exponent depends on several physical processes and has been a subject of debate over the years (Lucy 1967; Claret 2012; Espinosa Lara & Rieutord 2012). Recent observations showed that β can assume values of between 0.19 and 0.25 for typical B stars (van Belle 2012; McGill et al. 2013).

Assuming an effective gravitational potential Φ(r, θ) under the Roche approximation, the stellar surface is found from the equipotential surfaces (Collins & Harrington 1966; Tassoul 2007): $$\text{Φ}(r,\theta )=-\frac{GM}{r}-\frac{1}{2}{\text{Ω}}^2{\mathrm{\hspace{0.17em}sin}}^2\theta ,$$(2)

where G is the universal gravity constant, M stands for the stellar mass, and Ω is the angular velocity. Particularly at the poles, the effective potential reads (Townsend 1997): $$\text{Φ}=-\frac{GM}{R_{\text{p}}}=-g_{\text{p}}{R}_{\text{p}}.$$(3)

Applying Eq. (3) in Eq. (2), we get a cubic expression that delivers the shape of the stellar surface. If the star rotates with a rotation rate ω = Ω/Ω_c, the polar radius can be evaluated as (Collins & Harrington 1966): $$R_{\text{p}}(\omega )=R_{\text{p}}(0)\left(1-0.0540{\omega }^2-0.0547{\omega }^4\right),$$(4)

where the critical angular velocity Ω_c is given by (Collins 1965): $${\text{Ω}}_{\text{c}}={\left(\frac{8GM}{27R_{\text{p}}^3}\right)}^{1/2}.$$(5)

The ratio of the temperature of a non-rotating star T_s to the polar temperature T_p of a rotating star can be regarded as a function ξ of the stellar rotation rate ω: $$\frac{T_{\text{s}}}{T_{\text{p}}}=\xi {(}^{\omega },$$(6)

where T_s is the effective temperature of a non-rotating star and ξ(ω) is given as even powers of the fractional angular velocity ω (Cranmer 1996): $$\begin{array}{l}\xi (\omega )=1-0.19696{\omega }^2-0.094292{\omega }^4+0.33812{\omega }^6\hfill \\ -1.3066{\omega }^8+1.8286{\omega }^{10}-0.92714{\omega }^{12}.\hfill \end{array}$$(7)

At the stellar equator, the following constraints are remarkable: $$T_{\text{e}}=T_{\text{p}}{\left(\frac{g_e}{g_{\text{p}}}\right)}^{0.25},$$(8) $$g_{\text{e}}=g_{\text{p}}{\left(\frac{R_{\text{p}}}{R_{\text{e}}}\right)}^2-\frac{8}{27}{\omega }^2\frac{R_{\text{e}}}{R_{\text{p}}},$$(9)

and $$V_{\text{e}}=\omega \frac{R_{\text{e}}}{R_{\text{p}}}\sqrt{\frac{8}{27}\frac{GM}{R_{\text{p}}}},$$(10)

where T_e is the temperature, g_e is the surface gravity, and V_e is the rotation velocity. In any location of the stellar equipotential surface, the radius can be accomplished through Newton-Raphson iterations of the following relation: $$R(\theta ,\omega )=\frac{3R_{\text{p}}}{\omega \sin \theta }\cos \left(\frac{\pi +\arccos (\omega \sin \theta )}{3}\right).$$(11)

Rendering a rectangular mesh in spherical coordinates (e_r,e_θ,e_ϕ), for each surface element on the star, we can define the infinitesimal surface vector dA in terms of the radial profile R: $$\text{d}A=\left(\frac{\text{d}R}{\text{d}\theta }\times \frac{\text{d}R}{\text{d}\varphi }\right)\text{\hspace{0.17em}d}\theta \text{d}\varphi ,$$(12)

with $$\frac{\text{d}R}{\text{d}\theta }=R_{\theta }{e}_{\theta }+\frac{\text{d}{R}_{\theta }}{\text{d}\theta }{e}_r,$$(13)

and $$\frac{\text{d}R}{\text{d}\varphi }=R_{\theta }\sin \theta e_{\varphi },$$(14)

where R_θ is the θ component of the vector R. Substituting Eqs. (13) and (14) into Eq. (12) results in the area of each surface element: $$\text{d}A=R_{\theta }^2\sin \theta \text{d}\theta \text{d}\varphi \left(e_r-\frac{1}{R_{\theta }}\frac{\text{d}{R}_{\theta }}{\text{d}\theta }{e}_{\theta }\right),$$(15)

and its integration is performed only for those elements where the cosine of the local surface norm is positive; that is, for visible elements to the observer. This norm vector can be calculated by summing the centrifugal and gravitational forces, because this sum is orthogonal to the equipotential surface (Lovekin et al. 2006). As fast-rotating stars have emerging fluxes that reflect several different physical conditions on the stellar surface, it is important to choose some representative values that correspond to the parameters determined by fitting classical LTE/non-LTE synthetic spectra. Collins & Harrington (1966) chose, for instance, the weighted average of temperatures considering the intensity at the Banner’s discontinuity λ 3646 Å, but in principle, any other wavelength reference can be selected. It appears reasonable to adopt the expected value of the physical parameter (temperature or gravity) weighted by the intensity at the centre of the wavelength interval as a general measure of the entire photospheric physical state: $$<T_{\text{eff}}>=\frac{{\displaystyle \underset{\sum }{\int }T}(\theta ,\omega )\widehat{I}(\theta ,\omega ,\mu )\text{d}A}{{\displaystyle \underset{\sum }{\int }\widehat{I}}(\theta ,\omega ,\mu )\text{d}A},$$(16)

and $$<\log g_{\text{eff}}>=\frac{{\displaystyle \underset{\sum }{\int }{\text{log}}_g}(\theta ,\omega )\widehat{I}(\theta ,\omega ,\mu )\text{d}A}{{\displaystyle \underset{\sum }{\int }\widehat{I}}(\theta ,\omega ,\mu )\text{d}A},$$(17)

where the integrals are calculated only over all visible (μ > 0) regions on the stellar surface Σ and Î(θ, ω, μ) is the specific intensity of the continuum spectrum, which is taken, for instance, at the centre of the chosen wavelength interval.

3 Espinosa-Lara model

The Espinosa Lara & Rieutord (2011) model represents an improved gravity-darkening method based on the idea that the flux is antiparallel to the effective gravity, which allows the flux/gravity ratio to be preserved for all latitudes in stars rotating at least as fast as 90% of the break-up angular velocity. The flux expression is given by $$F=-\frac{L}{4\pi GM}{H}_{\omega }(\tilde{r},\theta )g_{\text{eff}}=-h(r,\theta )g_{\text{eff}},$$(18)

where Η_ω is a function that scales with the flux, $\tilde{r}$ is the non-dimensional radius, and ω is the non-dimensional rotation rate. Assuming no energy is generated outside the stellar core (Espinosa Lara & Rieutord 2011), $$\nabla \cdot F=g_{\text{eff}}\cdot \nabla h(r,\theta )+h(r,\theta )\nabla \cdot g_{\text{eff}}=0.$$(19)

Considering ∇ · g_eff = 2Ω² within the Roche approximation, Espinosa Lara & Rieutord (2011) obtained a new gravity-darkening expression by solving the following equation with the method of characteristics: $$\tau =\frac{1}{3}{\omega }^2{\tilde{r}}^3{\cos }^3\theta +\cos \theta +\text{ln\hspace{0.17em}tan}\frac{\theta }{2},$$(20)

finding the following solutions: $$\psi (\tau )={\tan }^2\vartheta ,$$(21)

and $$H_{\omega }=\frac{{\tan }^2\vartheta }{{\tan }^2\theta }.$$(22)

This expression defines the radiative flux for all surface points except in the singularities located at the poles and equator, whose solutions are (Espinosa Lara & Rieutord 2011) $$H_{\omega }(\theta =0)={\text{e}}^{\frac{2}{3}{\omega }^2{\tilde{r}}^3},$$(23)

and $$H_{\omega }\left(\theta =\frac{\pi }{2}\right)={\left(1-{\omega }^2{\tilde{r}}^3\right)}^{-\frac{2}{3}},$$(24)

where Eq. (22) gives a deviation measure from the classical case. The temperature distribution is given by $$\frac{T_{\text{eff }}(\text{equator})}{T_{\text{eff}}\left(\text{pole}\right)}={\left(1-{\omega }^2\right)}^{\frac{1}{12}}\sqrt{\frac{2}{2+{\omega }^2}}\exp \left(-\frac{4}{3}\frac{{\omega }^2}{{\left(2+{\omega }^2\right)}^3}\right)\text{. }$$(25)

The numerical integration of specific intensities is performed in the same way as discussed in the previous section, and also allows the inclusion of limb-darkening effects. Both the classical von Zeipel model and Espinosa-Lara’s new model agree in the low-angular-velocity regime, ω < 0.5, but strongly differ above it. The following input parameters are needed to run the Espinosa-Lara routine: stellar mass, polar radius, equatorial radius, polar temperature, inclination angle i, and the starting and ending wavelengths expressed in Å.

Fig. 1 Model spectra in the temperature range of B stars, evaluated with TLUSTY+SYNSPEC, showing a few illustrative models between Teff = 11 000 K and 33 000 K for log g = 4.0 evaluated with solar abundances, without rotation or microturbulence. By default, output fluxes evaluated with SYNSPEC are expressed as Eddington fluxes. Notable spectroscopic features of hydrogen and other elements are marked.

4 Plane-parallel input model spectra

The flux spectrum of a rotating star can be thought of as the combination of all specific intensities leaving the photosphere. Assuming the star can be divided into a regular square mesh, where each surface element has a particular set of parameters, namely local temperature T_eff, surface gravity log g, and Doppler velocity V (assuming ad hoc a solar chemical composition), then the overall flux distribution on the stellar surface will be given by the numerical integration of all visible specific intensities. For slow-rotating B stars, the calculations can be performed from a set of atmosphere structure models within the range of 11 000 K up to 33 000 K (see Fig. 1). For simulations of late/early-type B and O stars with Ω ≥ 0.7, sets of input spectra far beyond these values are usually required. Typically, a star rotating at ω ~ 0.99 has a polar-to-equator temperature difference that can reach around 8000 K, and in the case of late-type B stars, their equatorial zones can achieve the physical conditions of an A star while the rest of the star has the physical properties of a B-type star. In order to simulate these conditions, we built a grid with effective temperatures ranging from 4000 K to 50 000 K and with log g from 2.0 to 5.0. The low-temperature spectral set of between 4000 K and 10 000 K is synthesised from model-atmosphere structures by Kurucz (1979) with SYNSPEC v.54 (Hubeny & Lanz 1995).

Each surface element has its specific intensity interpolated from an extensive grid of local thermodynamical equilibrium (LTE) and/or non-local thermodynamical equilibrium (non-LTE) models that we computed with SYNSPEC and TLUSTY codes (Hubeny 1988; Hubeny & Lanz 1995). Temperatures higher than 15 000 K (up to 50 000 K) have atmosphere structures in full non-LTE computed with TLUSTY v.208 and are synthesised with SYNSPEC v.54 using the BSTAR2006 and OSTAR2002 atmosphere structure grids. Models between 11 000 K and 15 000 K are de novo computed in LTE using TLUSTY and explicitly taking into account the electronic transitions.

During the synthesis procedure, the code accepts any preformatted spectroscopic grid, both in flux units and normalised-to-continuum spectra. In cases where the user wants to input SYNSPEC/TLUSTY specific intensities, the following limb-darkening laws can be used (Levenhagen 2014): respectively, linear (Gray 1976), quadratic (Kopal 1950), square-root (Diaz-Cordoves & Gimenez 1992), or logarithmic (Klinglesmith & Sobieski 1970): $$I(\mu )=I(1)[1-\epsilon (1-\mu )],$$(26) $$I(\mu )=I(1)\left[1-\epsilon (1-\mu )-\varpi {(1-\mu )}^2\right]\text{, }$$(27) $$I(\mu )=I(1)[1-\epsilon (1-\mu )-\varpi (1-\sqrt{\mu })]\text{, }$$(28)

and $$I(\mu )=I(1)[1-\epsilon (1-\mu )-\varpi \mu \ln \mu ]\text{, }$$(29)

where ε and ϖ are limb-darkening coefficients dependent on the wavelength, temperature, and gravity. These coefficients are interpolated from the tables from (Wade & Rucinski 1985) for linear and quadratic laws. For the other laws, we interpolate tables using the JKTEBOP code (Southworth et al. 2004) to match the nodes of our stellar atmosphere grid. The code outputs are the emergent line spectrum and continuum expressed in Eddington flux units evaluated through angular integration of all visible surface elements. Also, an output normalised to the continuum spectrum is provided.

5 Colatitudinal structure

In what follows, we adopt the same notation as employed by Frémat et al. (2005), and refer to the input ‘deformless’ set of parameters as parent non-rotating counterpart (PNRC) parameters, which are meant to be the parameters of a non-deformed star with the same mass as the observed star. We refer to the set of parameters that we measure in the star as apparent parameters; these are usually obtained by fitting classical stellar atmosphere models to the observed star. Once the static star begins to spin, the average values of its surface parameters deviate from the original PNRC ones. Figure 2 shows the variation of surface temperatures for both the classical Von Zeipel and Espinosa-Lara models. In the classical case, we have widespread temperature and surface gravity colatitudinal distributions compared to the Espinosa-Lara models.

Considering a fast-rotating star, due to the geometric deformation induced by rotation, the equatorial radius gradually increases as the rotation rate approaches its critical value (ω_c). Individual hydrogen, helium, and metal line transitions begin to broaden by the Doppler effect, giving rise to blends of neighbouring lines, many of them formed at different depths in the stellar atmosphere. However, the increasing rotation rates ω imply steeper latitudinal temperature gradients on the stellar surface, and the combination of specific intensities coming from zones with very different temperatures can lead to asymmetries or features in certain line profiles.

To highlight these rotational effects on the line profiles, we present the synthesis of some He I lines commonly observed in the optical domain. In Fig. 3, we consider four helium lines calculated with ZPEKTR assuming the classical von Zeipel (represented with continuous solid lines) and Espinosa-Lara models (dashed lines) for a typical B5 star. Late-type fast rotators up to spectral subtype B5 with rotation rates of ω > 0.95 are limiting cases of interest because their geometrical shapes are oblate enough to allow extreme longitudinal temperature distributions, with hot poles at physical conditions of a Β star and reaching temperature values typical of A-type stars at the equatorial zones (between 8000 and 10 000 K), where He I lines no longer exist because of a change in ionisation regime. Examples of this phenomenon can be found in some fast-rotating late-type stars, such as ω Carinae (B8 Hie) and ο Andromedae (B6 Illpe). These objects display central emission peaks or central reversals that are attributed to the circumstellar environment (Rivinius et al. 1999). Even with the influence of circumstellar emissions, it is not excluded that the origin of central emission within helium lines could be, in these cases, partly due to the influence of gravity darkening.

In these models, we assume a PNRC temperature of Τ = 15 000 Κ, stellar mass 4.4 Μ/Μ_⊙, PNRC radius 6.18 R/R_⊙, aspect angle i = 35°, and five different rotation rates ω = Ω/Ω_c. All the He I spectral lines considered suffer from line blendings with other chemical species, but the He I 6678 Å is less contaminated by lines of other chemical elements with appreciable equivalent widths than the other helium lines, which would be responsible for its rather regular shape through all rotation rates.

Observing the He I 4388 Å line profile in Fig. 3, the classical von Zeipel model only shows a central reversal for ω = 0.985; it does not appear for lower rotation rates. The rest of the lines do not exhibit central reversals. On the other hand, for the Espinosa-Lara models (dashed lines), we do not find central reversals at any rotation rates. One explanation for central reversals in the He I 4388 Å line profile could be non-LTE line-formation processes in a low-gas-density atmosphere. Indeed, the He I 4388 Å line profile arises from a singlet transition 2¹P⁰ – 5¹D with an energy difference between these levels of higher than for the singlets He I 4922 Å (2¹P – 4¹D) and He I 6678 Å (2¹P – 3¹D), which could in principle favour the non-LTE processes in this line. Another scenario explaining the occurrence of central reversals involves the latitudinal temperature distributions in Fig. 2. The classical von Zeipel model considers, in this case, an equatorial temperature near 10 000 K, where no He I lines are formed. In this scenario, the central reversal would be formed by summing the flux contributions coming from hotter latitudes with cooler ones (in the absence of He I) in the same wavelength range. The equivalent Espinosa-Lara model, assuming the same geometrical deformation (i.e. polar and equatorial radii) and considering the same rotation rates in Fig. 2, for T (PNRC) = 15 000 K, shows an equatorial temperature of Τ ≃ 19 000 Κ and does not show central reversals. In this case, He I lines are formed in all latitudes. Also, for low rotation rates, the profiles are deeper at the line core than in the classical von Zeipel model.

It is worth remembering that the highest velocity and lowest flux are found in the equatorial regions. As such, these regions have the smallest contribution to the emergent spectrum. This fact implies that classical plane-parallel models cannot be used to estimate the projected rotation velocities of fast-rotating stars and result in underestimation of ν sin i.

Fig. 2 Colatitudinal variation of surface temperatures for the classical von Zeipel models with β = 0.25 (left panels) and the Espinosa-Lara models (right panels). Top panels show models assuming PNRC Teff = 24000 K, log g = 4.0, and ω = Ω/Ωc = 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, and 0.99. The bottom panels show models with PNRC Teff = 15 000 K, log g = 4.0, and the same rotation rates.

Fig. 3 Illustration of some classical von Zeipel (full lines) and Espinosa-Lara models (dashed lines) for He I 4471 + Mg II 4481, He I 4388, He I 4922, and He I 6678 Å. All models are built assuming a PNRC temperature of Τ = 15 000 Κ, M = 4.4 Μ⊙, R = 6.17 R⊙, and i = 35°, and showing the rotation rate ω = Ω/Ωc for three values, ω = 0.687, 0.879, and 0.985.

6 Comparison between rectangular and triangular meshes

As pointed out in a previous work (Levenhagen 2014), the synthesis procedure starts with the division of the stellar surface into small surface elements, each one assigned its own set of physical conditions. Many works on the subject customarily adopt the rectangular mesh in spherical coordinates; for example, the BRUCE code (Townsend et al. 2004; Townsend 1997). A possible issue with this grid is that the tiles or elements in the stellar surface change their dimensions with latitude. Near the stellar equator, the zone elements are large and decrease their area towards the poles. As the flux that comes through each tile is weighted by the area, the contribution of the fluxes from the equatorial zones is more important than the fluxes coming from the poles. As the rotation rate increases, the tiles near the equator suffer geometrical deformation, increasing more than at the poles.

An alternative procedure is to divide the stellar surface into triangular zones with approximately the same size in all stellar regions from the poles to the equator. The use of triangular tiles has been explored in the SUPEROTATE code (Dall & Sbordone 2011) to model the high oblate spectrum of Altair. The ZPEKTR code operates with both schemes and calculates the triangular tiles using the Delaunay triangulation scheme (Delaunay et al. 1934; see Fig. 4).

We compute models with ZPEKTR for a B5 star with a PNRC temperature of Τ = 15000 Κ, stellar mass 4.4 Μ/Μ_⊙, PNRC radius 6.18 R/R_⊙, inclination angle i = 35°, and rotation rates ω = Ω/Ω_c = 0.687, 0.792, 0.879, 0.944, and 0.985 for both the rectangular and triangular meshes. The syntheses are performed on a large wavelength range, from 4000 Å to 5000 Å.

The absolute differences between the rectangular and the triangular models can be seen in Fig. 5. The largest differences occur at the cores of the Banner lines, slightly decreasing from the Hδ towards Hβ but reaching at most 0.7% of continuum units. For helium lines, these differences are lower, at most 0.6%, while for other chemical species, we have even lower discrepancies of less than 0.4%. Given these results, we think it is safe to use the rectangular grid if we incorporate a sufficient number of square elements to secure the resolution of the line profile.

The number of tiles on the stellar surface to be considered during the synthesis is an important parameter and influences both the processing time and the line shapes. We performed a test on the mesh resolution, choosing a Β5 star model with T(PNRC) = 15 000 Κ again, stellar mass 4.4 Μ/Μ_⊙, PNRC radius 6.18 R/R_⊙, inclination angle i = 35°, rotation rate ω = Ω/Ω_c = 0.9, and four different numbers of tiles (N), from Ν = 300 up to 320000. We set a reference model spectrum with Ν = 500 000 tiles and built up the residuals for each model, as shown in Fig. 6. The residuals built with Ν = 300 tiles show differences in the Balmer line cores of around 0.8% with respect to the continuum and of up to 0.6% in the core of helium lines.

Increasing to N = 180000 tiles, we see a substantial improvement in the overall spectrum, with the residuals in Banner lines decreasing to 0.025% and less than 0.01% for helium lines. At this resolution level and beyond, the shapes and cores of the hydrogen and helium line profiles do not change perceptibly. Beyond Ν = 320 000 tiles, the residuals drop below 0.01% for Balmer lines and are lower than 0.005% in helium lines.

Fig. 4 ZPEKTR grid scheme, assuming Delaunay triangulation on one hemisphere of a spheroid.

Fig. 5 Top panel: Model spectra built with ZPEKTR using the triangular mesh (blue) and the rectangular mesh (orange). The two spectra are almost coincident, and minimal differences are detected only at the line cores of strong lines. Bottom panel: The absolute differences between these two spectra are shown in the bottom plot and are expressed in continuum percentages.

Fig. 6 Model spectra of a B5 star considering four levels of resolution, from Ν = 300 (first panel) up to Ν = 320 000 tiles (bottom panel). Again, the absolute differences with respect to a reference model with Ν = 500 000 tiles are expressed in continuum percentages. It can be seen that, for models with more than Ν = 180 000 tiles, the residuals are lower than 0.025% within the line cores.

7 Application to a small sample of Be stars

In order to illustrate the application of this code, we performed the fittings of classical von Zeipel models to a small sample of 31 Be stars observed with the PUCHEROS spectrograph (Vanzi et al. 2012). Further detailed comparisons between the classical prescription and Espinosa Lara’s approach to larger samples will be treated separately in a forthcoming publication.

The observed spectra are available at the Be Stars Observation Survey (BeSOS), which is a large catalogue with a resolution of R = 17 000 (Arcos et al. 2018), spanning a wavelength range 4200–7300 Å up to magnitude V = 9.

Assuming the classical von Zeipel formalism with ZPEKTR, we built up a large model grid with 15 561 models spanning from 4000 Å up to 8000 Å and Δλ = 0.1 Å. The ranges of parameters for the model grid are 1 ≤ Μ/Μ_⊙ ≤ 14 with ΔΜ/Μ_⊙ = 0.5, 2.5 ≤ R/R_⊙ ≤ 11 with ΔR = 0.5, 12000K ≤ T_eff ≤ 30 000K with ΔT = 100 K, 0.5 ≤ Ω/Ω_c ≤ 0.95 with ΔΩ = 0.05 and 0° ≤ i ≤ 90° with Δi = 2°. Figure 7 shows the fittings of gravity-darkened models to four Be stars of our sample. We chose to perform the fittings around the He I 4471 Å line profile because it is one of the lines least affected by circumstellar emission. The general results are summarised in Table 1, where T(NRPC) and log g(NRPC) are the effective temperature and surface gravity, respectively, of the non-rotating parent counterpart, which is an object with the same mass as the observed target but without geometrical deformation and homogeneous surface temperatures and gravities (Frémat et al. 2005). T_app and log g_app are the apparent effective temperature and surface gravity – which are the parameters we can infer from the fittings of classical model atmospheres (Hubeny 1988; Hubeny & Lanz 1995) –, i is the inclination angle, and Ω/Ω_c is the rotation rate. These parameters were inferred by fitting the observed spectra to the model grid using χ² goodness-of-fit test.

According to recent studies (Frémat et al. 2005; Zorec et al. 2017a), fast rotation triggers a geometrical flattening, which lowers the stellar luminosity in the inner layers. The faster the star rotates, the more the internal pressure in the stellar layers decreases, mainly in the equatorial regions. This effect decreases the nuclear reaction rate at the stellar core and makes the star remain for longer on the main sequence (Frémat et al. 2005). Figure 8 shows the resulting HR diagrams of the program stars with the fittings of classical model atmospheres or apparent temperatures and luminosities (left panel) and accounting for geometrical flattening and gravity darkening (PNRC parameters, right panel). The stars on the left panel appear more evolved in their trajectories in the HR diagram than their counterparts on the right. Also, stars with their projected rotation velocities estimated through the fitting with the classical model atmosphere models have smaller average υ sin i values than the stars for which the rotation velocity was corrected for gravity darkening (see Fig. 8).

The inference of the probability density function (PDF) of the true rotation velocities from the projected stellar velocities is an ill-posed problem, and we used the algorithm proposed by Lucy (1974). Assuming that the rotation axes of the stars in the sample are randomly distributed, then the discrete probability matrix governing the distribution has its elements associated with the cumulative distribution function given by the integral of the following continuous probability distribution (Lucy 1974): $$P(x\mid \xi )=\frac{x}{\xi }{\left({\xi }^2-x^2\right)}^{-1/2}H(\xi -x),$$(30)

where Η (ξ – x) = 1 if ξ – x ≥ 0 and Η(ξ – x) = 0 if ξ – x < 0. The numerical iterative scheme of Lucy’s algorithm reads $${\psi }_j^{r+1}={\psi }_j^r{\displaystyle \sum _{i=1}^N\frac{{\widehat{\varphi }}_i}{{\varphi }_i}}P(i\mid j),$$(31)

where $\widehat{\varphi }=n\left(x_i\right)/N$ is the initial guess at the histogram of the sampled n(x_i) estimations over the Ν total observed ν sin i values that fall in the closed interval [x_i – 1/2Δx_i, x_i + 1/2Δx_i]. The ψ^r+1 is the r + 1 iteration of the rectified PDF of true velocities and the i, j indices correspond to the x_i =v sin i and true ξ_j = ν values, respectively. The first guess of ${\widehat{\varphi }}_i$ is obtained from the raw data through kernel density estimation (Silverman 1986). Figure 9 shows the histograms of the raw apparent υ sin i (blue histogram) and the histogram of true velocities ν obtained from the fittings with ZPEKTR to the data (red histogram). The dashed blue curve represents the PDF φ of the apparent ν sin i distribution obtained with kernel density estimation. In contrast, the dash-dotted orange curve represents the PDF ψ associated with the rectified distribution. Even though the analysis has been applied to a relatively small stellar sample, and as a consequence their results should be taken carefully, one can see that there is a tendency for apparent velocities to be slightly lower than those obtained with ZPEKTR. Interestingly, the rectified curve ψ seems to follow the pattern of the red histogram of estimates with ZPEKTR. This is a long-known phenomenon (Bernacca 1972), but it serves to illustrate the possibilities of the code itself.

Fig. 7 Line fittings of He I 4471 Å gravity-darkened models with the classical von Zeipel model to the Be stars HD 33328, HD 37795, HD 124771, and HD 158427 from a sample of 31 stars of the BeSOS database. The spectra represent the medians and also show the Mg II 4481 Å line profile, but to determine the PNRC parameters, only the He I 4471 Å profile has been taken into account in the fittings.

8 Discussion

The ZPEKTR code presented in this work was designed following both the general outlines proposed by Collins & Truax (1995) for von Zeipel’s classical picture and the novel formalism by Espinosa Lara & Rieutord (2011). In principle, any grid of synthetic spectra can be adapted to be used as input for the code. In this work, we use a plane-parallel grid of non-LTE spectra as input, which was built with SYNSPEC v.54 departing from model atmospheres generated with TLUSTY v.208 Hubeny (1988); Hubeny & Lanz (1995) for hot stars above 15 000 Κ and a grid of LTE spectra by Kurucz (1979) for cooler stars. Figure 1 shows some illustrative input Eddington flux spectra from this grid in units of erg cm⁻² s⁻¹ Å⁻¹.

The latitudinal distribution of temperatures over the stellar surface -represented in Fig. 2 for T (PNRC) = 15 000 Κ and T(PNRC) = 24000 Κ – is larger for the classical von Zeipel models than for the Espinosa-Lara models. The same occurs for the surface gravity distribution. This is reflected in a depth variation of the helium line profiles shown in Fig. 3. A direct comparison of normalised helium line profiles also shows an increasing dependence of the shape of the helium line profiles on the rotation rate. Helium lines, in general, appear blended with lines of other chemical species, such as oxygen, and the resulting profiles may become very asymmetric with the increasing rotation rate.

Table 1 summarises the stellar fundamental parameters obtained with the fittings of classical and gravity-darkened models. Figure 8 shows the position of the Be star sample in two HR diagrams; the top panel shows the apparent stellar parameter, and the bottom panel shows the PNRC temperatures. The evolutionary tracks and isochrones in both diagrams were calculated by Schaller et al. (1992) for static stars; that is, without the assumption of fast rotation in the solution of stellar structure equations. Generally, for fast-rotating stars, PNRC temperatures are higher than the apparent ones (Frémat et al. 2005). The same occurs for the surface gravities, except when the star rotates slowly; that is, at low rotation rates. These differences become more pronounced as we increase the rotation rate and the inclination angle, ultimately increasing the v sin i values. The same effect was found by Frémat et al. (2005).

In Fig. 9, we see a comparison between apparent and ‘true’ v sin i estimations obtained from the fittings to the BeSOS sample. The apparent v sin i histogram (see Fig. 9, blue bins) was built from the fittings of plane-parallel models (SYN-SPEC/Kurucz) to the sample stars, and is shown together with its PDF distribution (dashed blue line). The ‘true’ v sin i histogram (red bins) shows the velocity values obtained through the fittings of von Zeipel models to the stellar sample. These values are higher than the apparent ones and should correspond to a better v sin i estimation, since the fitting models encompass gravity darkening effects, which affect the fast equatorial zones. Supposing we do not know anything about the red histogram, and only the blue histogram is accessible, we can still infer the ‘true’ PDF through Bayesian statistics. Assuming the stellar rotational axes are randomly distributed, we apply the Lucy rectification scheme (Lucy 1974) using as input the apparent PDF and obtain the red dotted-line curve. This rectified distribution closely follows the histogram of von Zeipel estimates and shows that models without the assumption of gravity darkening tend to underestimate v sin i values, as also pointed out by Frémat et al. (2005).

9 Conclusions

In this paper, we present the code ZPEKTR, aimed to perform the spectral synthesis of fast-rotating stars. The code currently supports two formalisms, the classical formulation by von Zeipel (1924a,b) and the new prescription by Espinosa Lara & Rieutord (2011). Comparing the latitudinal profiles in the two formalisms, we see that in the classical approach, the temperature variation departing from the pole to the equator is larger than that produced with the Espinosa-Lara formulation.

The same effect is seen with the latitudinal distribution of surface gravity. Also, it is worth mentioning that the polar compressibility of the star, concomitantly with equatorial enlargement due to fast rotation, does not preserve the initial stellar volume. Papaloizou & Whelan (1973) estimated the polar radius to decrease by 3% at most, while the equatorial radius increases by around 50%. This large gas expansion in the equatorial zones of critical rotators can reduce the gas density at those regions and the amount of collisional processes in the upper layers of the atmosphere. This would lead to non-LTE radiative processes whose signatures could appear as line core reversals or central emission peaks in higher-level transitions.

The Espinosa-Lara formalism represents a significant improvement over the previous attempts at the theory, including both the requisition of radiative equilibrium among the atmospheric layers and the null divergence of the radiative flux outside the stellar core. By introducing a correction term to the classical formalism, it takes into account the Eddington Sweet flows inside the star. We stress that it is not our goal to compare the formalisms here, as was done by Espinosa-Lara himself in other works (Espinosa Lara & Rieutord 2011 ; Zorec et al. 2017b), but to present some capabilities of the code at this stage. Further explorations of the code will be presented in a subsequent work.

One of the main consequences of fast rotation in stellar evolution is the increase in the time spent by a star on the main sequence. In Fig. 8, we present the HR diagram of the program stars together with the evolutionary tracks and isochrones of non-rotating stars by Schaller et al. (1992). The left panel, representing the apparent parameters, shows the stars situated before the second half of the main sequence. The right panel represents the PNRC parameters and shows the stars displaced in the HR diagram towards regions nearer the ZAMS. Comparing the PNRC temperatures and gravities with the apparent parameters, it can be seen that their differences are small for slow-rotating stars seen at low inclination angles when Ω/Ω₀ < 0.5. Still, these differences become important for higher rotation rates. This behaviour, which has been reported for a larger sample of 130 Be stars (Frémat et al. 2005), occurs due to the decrease in nuclear energy generation rates induced by the reduction in gas pressure caused by centrifugal acceleration. Also, from Fig. 9, one can compare the histogram and PDF of the apparent υ sin i distribution with the shifted histogram of ν sin i values obtained with ZPEKTR, which coincides with the kernel-smoothed rectified PDF obtained from the apparent histogram through iterations with Lucy’s rectification method (Lucy 1974). From this plot, it can be seen that for fast rotation rates, the ν sin i parameter underestimates the stellar equatorial velocity (Frémat et al. 2005).

Another important consequence of gravity-darkening that is not explored in this work but is worth mentioning is the mixing of the outward-moving products of the CNO cycle due to meridional Eddington currents in the upper stellar layers and subsequent enrichment of the photosphere (Porter 1999; Frémat et al. 2005; Abdul-Masih 2023). In the near future, we intend to expand the code to implement a generalised formalism of the Espinosa-Lara model (Zorec et al. 2017b), and the treatment of tidal and radiative interactions in binary stellar systems (Abdul-Masih et al. 2020). We also plan to include differential-rotation effects.

Fig. 8 Top panel: The sample of 31 fast-rotating stars (blue dots) represented in the HR diagram shows the apparent Teff and log L/L⊙. parameters. Bottom panel: HR diagram for the same sample, assuming PNRC parameters obtained with the classical von Zeipel model. Black dashed lines represent isochrones with 107 and 108 yr; coloured continuum lines represent evolutionary tracks from 2 M/M⊙ up to 20 M/M⊙. The evolutionary tracks and isochrones were calculated by Schaller et al. (1992) supposing static (non-rotating) internal structure equations.

Fig. 9 Blue bins: Apparent υ sin i measurements. Red bins: PNRC υ sin i estimations with ZPEKTR. Dashed blue: φ estimate with kernel smooth at the apparent histogram. Dashed-dot: rectified ψ PDF distribution using Lucy’s algorithm.

Acknowledgements

This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00. M.C. and C.A. acknowledge partial support from Centro de Astrofisica de Valparaiso. M.C. also acknowledges support from Centro Interdiciplinario de Estudios Atmosfericos y Astroestadistica, Universidad de Valparaiso. This work used BeSOS Catalogue, operated by the Instituto de Fisica y Astronomia, Universidad de Valparaiso, Chile: https://besos.ifa.uv.cl and funded by Fondecyt iniciacion no. 11130702. M.C., C.A. and I.A. thank the support from FONDECYT project N.1230131. This work has been possible thanks to the use of AWS-U.Chile-NLHPC credits. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02). This project has received funding from the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under the Marie Sklodowska-Curie Grant Agreement No. 823734. E.B.A thanks Universidade Estadual de Feira de Santana for the support received by the Program FINAPESQ (project number 050/2021).

References

All Tables

All Figures

Fig. 1 Model spectra in the temperature range of B stars, evaluated with TLUSTY+SYNSPEC, showing a few illustrative models between Teff = 11 000 K and 33 000 K for log g = 4.0 evaluated with solar abundances, without rotation or microturbulence. By default, output fluxes evaluated with SYNSPEC are expressed as Eddington fluxes. Notable spectroscopic features of hydrogen and other elements are marked.

In the text

	Fig. 2 Colatitudinal variation of surface temperatures for the classical von Zeipel models with β = 0.25 (left panels) and the Espinosa-Lara models (right panels). Top panels show models assuming PNRC Teff = 24000 K, log g = 4.0, and ω = Ω/Ωc = 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, and 0.99. The bottom panels show models with PNRC Teff = 15 000 K, log g = 4.0, and the same rotation rates.
In the text

	Fig. 3 Illustration of some classical von Zeipel (full lines) and Espinosa-Lara models (dashed lines) for He I 4471 + Mg II 4481, He I 4388, He I 4922, and He I 6678 Å. All models are built assuming a PNRC temperature of Τ = 15 000 Κ, M = 4.4 Μ⊙, R = 6.17 R⊙, and i = 35°, and showing the rotation rate ω = Ω/Ωc for three values, ω = 0.687, 0.879, and 0.985.
In the text

	Fig. 4 ZPEKTR grid scheme, assuming Delaunay triangulation on one hemisphere of a spheroid.
In the text

	Fig. 5 Top panel: Model spectra built with ZPEKTR using the triangular mesh (blue) and the rectangular mesh (orange). The two spectra are almost coincident, and minimal differences are detected only at the line cores of strong lines. Bottom panel: The absolute differences between these two spectra are shown in the bottom plot and are expressed in continuum percentages.
In the text

	Fig. 6 Model spectra of a B5 star considering four levels of resolution, from Ν = 300 (first panel) up to Ν = 320 000 tiles (bottom panel). Again, the absolute differences with respect to a reference model with Ν = 500 000 tiles are expressed in continuum percentages. It can be seen that, for models with more than Ν = 180 000 tiles, the residuals are lower than 0.025% within the line cores.
In the text

	Fig. 7 Line fittings of He I 4471 Å gravity-darkened models with the classical von Zeipel model to the Be stars HD 33328, HD 37795, HD 124771, and HD 158427 from a sample of 31 stars of the BeSOS database. The spectra represent the medians and also show the Mg II 4481 Å line profile, but to determine the PNRC parameters, only the He I 4471 Å profile has been taken into account in the fittings.
In the text

Fig. 8 Top panel: The sample of 31 fast-rotating stars (blue dots) represented in the HR diagram shows the apparent Teff and log L/L⊙. parameters. Bottom panel: HR diagram for the same sample, assuming PNRC parameters obtained with the classical von Zeipel model. Black dashed lines represent isochrones with 107 and 108 yr; coloured continuum lines represent evolutionary tracks from 2 M/M⊙ up to 20 M/M⊙. The evolutionary tracks and isochrones were calculated by Schaller et al. (1992) supposing static (non-rotating) internal structure equations.

In the text

	Fig. 9 Blue bins: Apparent υ sin i measurements. Red bins: PNRC υ sin i estimations with ZPEKTR. Dashed blue: φ estimate with kernel smooth at the apparent histogram. Dashed-dot: rectified ψ PDF distribution using Lucy’s algorithm.
In the text