Mode visibilities in rapidly rotating stars

Mode visibilities in rapidly rotating stars

D. R. Reese Institut d’Astrophysique et Géophysique de l’Université de Liège, Allée du 6 Août 17, 4000 Liège, Belgium
daniel.reese@ulg.ac.be LESIA, Observatoire de Paris, CNRS, UPMC Univ. Paris 06, Univ. Paris-Diderot, 5 place Jules Janssen, 92195 Meudon, France Kavli Institute for Theoretical Physics, Kohn Hall, University of California, Santa Barbara, CA 93106, USA
   V. Prat Universté de Toulouse, UPS-OMP, IRAP, Toulouse, France CNRS, IRAP, 14 avenue Edouard Belin, 31400 Toulouse, France    C. Barban LESIA, Observatoire de Paris, CNRS, UPMC Univ. Paris 06, Univ. Paris-Diderot, 5 place Jules Janssen, 92195 Meudon, France    C. van ’t Veer-Menneret GEPI, Observatoire de Paris-Meudon, CNRS, Université Paris Diderot, 92125, Meudon Cedex, France    K. B. MacGregor High Altitude Observatory, National Center for Atmospheric Research, Boulder, CO 80307, USA
Key Words.:
stars: oscillations (including pulsations) – stars: rotation
Abstract

Context:Mode identification is a crucial step to comparing observed frequencies with theoretical ones. However, it has proven to be particularly difficult in rapidly rotating stars. An important reason for this is the lack of simple frequency patterns such as those present in solar-type pulsators. This problem is further aggravated in Scuti stars by their particularly rich frequency spectra.

Aims:As a first step to obtaining further observational constraints towards mode identification in rapid rotators, we aim to accurately calculate mode visibilities and amplitude ratios while fully taking into account the effects of rotation.

Methods:We derive the relevant equations for calculating mode visibilities in different photometric bands while fully taking into account the geometric distortion from both the centrifugal deformation and the pulsation modes, the variations in effective gravity, and an approximate treatment of the temperature variations, given the adiabatic nature of the pulsation modes. These equations are then applied to 2D oscillation modes, calculated using the TOP code (Two-dimension Oscillation Program), in fully distorted 2D models based on the SCF (Self-Consistent Field) method. The specific intensities come from a grid of Kurucz atmospheres, thereby taking into account limb and gravity darkening.

Results:We obtain mode visibilities and amplitude ratios for models with rotation rates ranging from to of the critical rotation rate. Based on these calculations, we confirm a number of results from earlier studies, such as the increased visibility of numerous chaotic modes at sufficient rotation rates, the simpler frequency spectra with dominant island modes for pole-on configurations, or the dependence of amplitude ratios on inclination and azimuthal order in rotating stars. In addition, we explain how the geometric shape of the star leads to a smaller contrast between pole-on and equator-on visibilities of equatorially-focused island modes. We also show that modes with similar values frequently have similar amplitude ratios, even in the most rapidly rotating models.

Conclusions:

1 Introduction

The space missions CoRoT (Baglin2009; Auvergne2009) and Kepler (Borucki2009) are revealing very rich pulsational spectra in rapidly rotating Scuti stars. For instance, several hundred individual frequencies have been found in HD 50844 and HD 181555, observed by CoRoT (Poretti2009, Michel, private communication), and V2367 Cyg, observed by Kepler (Balona2012). It is becoming increasingly clear that interpreting these spectra will not be a straightforward task and that theory is lagging behind observations. A crucial first step in interpreting this data is correctly identifying the pulsation modes, i.e. finding the correct correspondence between theoretically calculated modes and observed pulsations. Recently, Reese2009b proposed a way to identify acoustic pulsation modes in rapidly rotating stars based on an asymptotic formula which describes the frequencies of low degree modes (see Pasek2012, and references therein). Nonetheless this method runs into trouble if chaotic modes are present in the pulsation spectra, which is expected based on the visibility calculations in Lignieres2009. Furthermore, the pulsation modes in Scuti stars tend to be of low radial order, and may therefore be too far from the asymptotic regime. Lignieres2010 have worked on using the cross-correlation of pulsation spectra. Although it doesn’t yield individual mode identifications, it may provide a way of obtaining the rotation rate and/or the large frequency separation, and explaining recent observations of recurring frequency spacings in rapid rotators (GarciaHernandez2009; Mantegazza2012). Nonetheless, the need remains for methods capable of identifying individual pulsation modes.

Two particularly promising methods for identifying pulsation modes are multi-colour photometric and spectroscopic mode identification. The first approach consists in measuring the amplitudes and phases of a given pulsation mode in different photometric bands, calculating the ratios of the different amplitudes and/or the phase differences, and comparing these to theoretical predictions. The second approach exploits the Doppler shifts caused by the velocity field from the pulsation mode and how it affects observed absorption lines. These methods have been successfully applied to slowly rotating stars (e.g.  DeRidder2004; Zima2006; Briquet2007), but more work is needed before they are applied to rapid rotators. In the present paper, we will focus on mode visibilities in different photometric bands as a first step to multi-colour photometric mode identification, and postpone spectroscopic mode identification to a later paper.

Few studies have dealt with multi-colour photometric mode signatures in rapidly rotating stars, and those that do generally approximate the effects of rotation on the pulsation modes. For instance, Daszynska_Daszkiewicz2002; Daszynska_Daszkiewicz2007 and Townsend2003b used either the perturbative approach or the traditional approximation to calculate their pulsation modes, and in some cases included the effects of avoided crossings. Their calculations of mode visibilities included stellar surface distortion from the pulsation modes and the Lagrangian perturbations to both the effective temperature and gravity. An interesting result from these studies is that contrarily to non-rotating stars, amplitude ratios depend both on , the azimuthal order of the pulsation mode, and , the inclination of the star. More recently, Lignieres2006 and Lignieres2009 calculated geometrical disk-integration factors of acoustic modes in deformed polytropic models by integrating the temperature fluctuations over the visible disk. The effects of rotation were fully taken into account in the pulsation modes, thanks to the 2D numerical approach, but non-adiabatic effects were neglected, thereby making the fluctuations of the effective temperature inaccessible. Furthermore, gravity and limb darkening, the Lagrangian perturbations to the effective gravity, and surface distortion caused by the modes were not taken into account. Nonetheless, important first results were obtained through these articles, namely, that chaotic modes are more visible than their non-rotating counterparts due to irregular latitudinal node placement and may thus be detected, island modes are the most visible modes in a pole-on configuration, leading to a regular frequency pattern (Lignieres2009), and signatures of the large frequency separation and/or rotation rate can show up in the autocorrelation function of the frequency spectrum (Lignieres2010).

In order to obtain more realistic multi-colour photometric mode visibilities in rapidly rotating stars, we derive a new set of equations which take into account the Lagrangian variations to the effective temperature and gravity, as well as the surface distortions induced both by the centrifugal deformation of the equilibrium model and by the pulsation modes. These are applied to adiabatic acoustic modes calculated by the Two-dimension Oscillation Program (TOP, Reese2006; Reese2009a) using rapidly rotating zero-age main-sequence (ZAMS) models based on the Self-Consistent Field (SCF) method (Jackson2005; MacGregor2007). The emergent intensities are calculated from Kurucz atmospheres, taking into account the latitudinal dependence of the equilibrium effective temperature and gravity, thereby including gravity and limb darkening. The main weakness in the present study is the adiabatic approximation, which makes the Lagrangian fluctuations of the effective temperature inaccessible. As was previously done in Lignieres2009, we approximate these by the Lagrangian temperature variations. The following section describes the pulsation calculations, with an emphasis on the improvements and differences with the calculations done in Reese2009a. This section is followed by a derivation of the relevant equations for calculating mode visibilities in rapidly rotating stars. Section 4 then describes various effects of rotation on visibilities in a single band – the CoRoT photometric band. This is followed by Sect. 5 which deals with amplitude ratios in the Geneva photometric system. The paper ends with a short conclusion.

2 Pulsation calculations

In what follows, we review the methods used for obtaining the models and associated pulsations that serve as inputs to the visibility calculations. These closely follow the approach used in Reese2009a but include a number of improvements as described below.

2.1 Equilibrium models

The equilibrium models are calculated via the SCF method (Jackson2005; MacGregor2007). This method is an iterative procedure which alternates between solving Poisson’s equation and the equations of mass, momentum and energy conservation before converging onto a 2D centrifugally deformed stellar model. These models are chemically homogeneous ZAMS models with a cylindrical rotation profile, although throughout the rest of the article, we will work with uniformly rotating SCF models, even if the formulas in the visibility calculations are established for general (non-cylindrical) rotation profiles. Given the rotation profile, the structure of the model is barotropic, i.e. all thermodynamic quantities remain constant on isopotentials, which are calculated from the sum of the gravitational and centrifugal potentials. Finally, we wish to make the distinction between the critical rotation rate, , and the Keplerian break-up rotation rate, :

(1)

Although very similar, the former uses the true gravity (excluding the centrifugal force) at the equator, , based on the actual distribution of matter, to calculate the break-up rotation rate, whereas the latter uses its Keplerian approximation, , which amounts to assuming spherical symmetry for the distribution of matter. As pointed out in Roxburgh2004, the Keplerian approximation slightly underestimates the true gravity, so that . Table 1 gives the relative differences between these two quantities for different rotation rates, calculated two different ways. The first method is based on global quantities provided with the models, whereas the second involves recalculating the gravitational potential from the density distribution and using this to calculate the equatorial gravity. A comparison of columns two and three, and five and six gives an idea of the uncertainty on these values. We also note that the theoretical value of at is , since the star is spherically symmetric and the Keplerian approximation is exact.

Method 1 Method 2 Method 1 Method 2
Table 1: Relative differences, , for selected rotation rates, calculated with two different methods (see text for details).

Before being used in the pulsation calculations, the models need to be interpolated onto a new grid, and a number of supplementary equilibrium quantities have to be derived, including a variety of geometric terms as well as gradients of different equilibrium quantities. Since Reese2009a, a number of improvements have been incorporated into these procedures. For instance, the stellar models are now interpolated onto a non-uniform radial grid which becomes dense near the stellar surface. This allows the pulsation code to correctly resolve the rapid spatial variations of acoustic modes near the surface, resulting from the decrease in sound velocity. Instead of interpolating the density and pressure directly, their logarithm is interpolated. This leads to more accurate and consistent values near the surface, where these quantities are several orders of magnitude smaller than in the centre, and ensures they remain positive. The effective gravity is calculated via Poisson’s equation. This avoids taking the ratio of the pressure gradient divided by the density, both of which are small quantities subject to relatively large uncertainties. Furthermore, the equipotentials, and hence the geometric structure of the star, are recalculated using the solution from Poisson’s equation thus removing some numerical inaccuracies in the original models. The derivative of equilibrium quantities is now correctly calculated. In Reese2009a, the derivative was mistakenly calculated with respect to rather than . Although this changed the quantitative results, the qualitative conclusions from Reese2009a remain unaltered. Finally, the profile is not derived from Eq. (16) of Jackson2005, but rather from the equation of state, which is based on the formula of Eggleton1973. Furthermore, the quantity , which intervenes in the visibility calculations described below, is also calculated via the equation of state.

2.2 Pulsation equations

A new set of variables is used in the pulsation equations:

(2)

where is the Lagrangian displacement, the Lagrangian pressure perturbation divided by the equilibrium pressure profile, the Lagrangian density perturbation divided by the equilibrium density profile, and the Eulerian perturbation to the gravitational potential. Throughout this article, the subscript “” denotes equilibrium quantities. We assume that the time and dependence of these variables takes on the form , where is the azimuthal order. As such, we use what could be called the “retrograde convention”, i.e. modes with positive azimuthal orders, , are retrograde.

Based on these variables, the continuity equation becomes

(3)

and Poisson’s equation is

(4)

where is the gravitational constant. Euler’s equation takes some more manipulations (see, for example, Eq. A.3 of Reese2009a):

(5)

where is the rotation profile, the distance to the rotation axis, and the effective gravity. The term in curly brackets cancels out because is parallel to in a barotropic stellar structure. Finally, the adiabatic relation takes on the following very simple form:

(6)

This last equation is then used to eliminate in favour of throughout the differential system and thus to reduce the size of the problem compared to what is obtained in Reese2009a. Explicit expressions for Eqs. (3), (4) and (5), using spheroidal coordinates (see Sect. 2.4), are given in App. A.1.

Besides reducing the computational cost, using the variables in Eq. (2) allows us to obtain a much cleaner derivation of , and hence , near the surface. Indeed, if one were to calculate this quantity from the Eulerian pressure perturbation, they would apply the following relation:

Near the surface, this involves the sum of two nearly opposite terms, divided by a small quantity, thus leading to poor numerical results.

2.3 Non-dimensionalisation

Contrarily to Reese2009a, we non-dimensionalise the SCF models in a more classical way. The following reference quantities are used as units of length, density and pressure:

(7)

where is the equatorial radius and the mass. These lead to the following time scale:

(8)

where is the Keplerian break-up rotation rate. Hence, the non-dimensional frequencies are directly . With this non-dimensionalisation, the preceding pulsation equations remain unchanged except for Poisson’s equation, which becomes

(9)

2.4 Spheroidal geometry

As was done in Lignieres2006 and Reese2006, a surface-fitting coordinate system, , based on Bonazzola1998, is introduced. This system is related to the usual spherical coordinates, , via the relation

(10)

where corresponds to the surface and is a measure of the oblateness, being the polar radius, and is comprised between and . As can be seen, corresponds to the stellar surface. A second domain, with a spherical outer boundary, is added around the star so as to simplify the boundary condition on the gravitational potential:

(11)

where . For , Eq. (11) coincides with the stellar surface, whereas for , it yields a sphere of radius (or in dimensional form). For conciseness, we will use the subscripts “” and “” to denote derivatives of with respect to these variables. For example, and .

2.5 Boundary conditions

The pulsation equations are supplemented by a number of boundary equations. Regularity of the solutions is imposed in the centre. The perturbation to the gravitational potential is made to go to zero at an infinite distance from the star. This condition is imposed by extending into the second domain and matching it to a vacuum potential on the second domain’s outer spherical boundary as described in Reese2006. When extending into the second domain, both it and its gradient need to be kept continuous across the perturbed stellar surface. This can be achieved by imposing the continuity of the Lagrangian perturbation to the gravitational potential and its gradient at :

(12)

where the superscripts “int” and “ext” correspond to “just below” and “just above” the stellar surface, respectively. Given that the gradient of the equilibrium gravitational potential is continuous, the first condition simplifies to

(13)

The second condition can be simplified, using Poisson’s equation, applied to the equilibrium gravitational potential:

(14)

Here, is the radial component of the Lagrangian displacement, when decomposed over the alternate basis introduced in Reese2006. The first equation is different from what was applied in Reese2009a, since it takes into account the contribution from a non-zero surface density. However, given its low value, the resultant difference is quite negligible. The latter two equations are implied by Eq. (13), making it unnecessary to impose them.

Finally, the usual mechanical boundary condition on the Lagrangian pressure perturbation, , has been replaced by a slightly different condition:

(15)

where is a vector perpendicular to the surface (see Eq. (51)). The above, modified mechanical condition corresponds to setting the vertical gradient of , rather than its value, to zero. We note that Pesnell1990 and Dupret2002 applied a similar condition for spherically symmetric stars. In order to avoid having a boundary condition with a radial derivative in it, we calculate the dot product between Euler’s equation and and cancel out the term corresponding to the vertical gradient of . This leads to a complicated expression which is given in spheroidal coordinates in App. A.2.

The main purpose in using Eq. (15) is to obtain a non-zero value for , useful for visibility calculations as described below. Indeed, when combined with a non-zero surface pressure and the adiabatic relation, the simpler condition, , leads to , as illustrated in Fig. 1 (left panels). One may then wonder if Eq. (15) has an important effect on the frequencies and on the displacement at the surface. Numerically, it turns out the frequencies vary little when using either boundary condition, at least in the present study. The middle panels of Fig. 1 also show that the displacement is hardly affected.

Figure 1: Lagrangian temperature variations and Lagrangian displacement at the stellar surface for two modes, using either or Eq. (15) as an external boundary condition. The figures to the right display the meridional cross-section of the Eulerian pressure perturbation of the two modes, divided by the square-root of the equilibrium density.

2.6 Numerical method

This set of equations and boundary conditions is projected onto the spherical harmonic basis and discretised in the radial direction before being solved using the code TOP (Reese2006; Reese2009a). Besides the improved way of treating the equilibrium model, the present calculations also benefit from a new form of finite differences. This form achieves 4th order accuracy for 1st order derivatives in spite of using windows with 4 rather 5 grid points. More importantly, this approach is robust to problems like mesh drift and spurious solutions.

3 Mode visibilities

At this point, we will switch to using the spherical vector basis in which the polar or z-axis is lined up with the star’s rotation axis. Furthermore, in order to make the equations more compact, we will prefer the notation “” to “” when designating the stellar surface, although it should be understood as in what follows. Similar implicit arguments also apply to other geometric terms such as and . We introduce the unit vector which points from the star to the observer. Furthermore, we will assume that the vector lies in the meridional plane . Let be the inclination angle, i.e. the angle between and , where is lined up with the rotation axis111As opposed to a non-rotating star, the inclination cannot be arbitrarily set to to simplify the calculations.. An explicit expression for in terms of the usual spherical basis is:

(16)

The radiated energy, received by an observing instrument from a non-pulsating stars, is:

(17)

where is the distance to the star, the cosine of the angle between the outward normal to the surface and , and the effective gravity and temperature, and the specific radiation intensity, multiplied by the instrument’s and/or filter’s transmission curve, and integrated over the wavelength spectrum,. As can be seen, the integral is carried out over the visible surface.

In a pulsating star, this quantity is perturbed as follows:

(18)

where denotes the Lagrangian perturbation, the real part, and the amount by which the visible surface is modified due the modifications of the surface normal induced by the oscillatory motions. Furthermore, we assume a complex form for the eigenfunctions, hence the reason for taking the real part of the above expression.

The first term in Eq. (18) is proportional to the square of the displacement and therefore neglected (e.g.  Dziembowski1977).

The Lagrangian perturbation to the specific intensity which intervenes in the second term may be developed as follows:

(19)

The partial derivatives of are obtained from model atmospheres and will be dealt with in Sect. 3.4. We note that the above expression is not exact as it neglects the slight Doppler shifts caused by rotation and the oscillations. Such shifts modify the position of the emerging flux with respect to the instrument’s and/or filter’s transmission curve thereby modifying , but are expected to play a negligible role compared to other effects. We are therefore left with a number of geometrical terms to calculate as well as the Lagrangian perturbation to the effective temperature and gravity.

3.1 Geometrical terms

In what follows, we will use the following expression for the displacement:

(20)

where is the usual vector basis in spherical coordinates. Note: we have used subscripts rather than superscripts, for the letters , , and , to distinguish these components from those given in Eq. (53). In what follows, we will be using the following relations:

(21)

We start by calculating a surface element on a rotating, non-pulsating star. This is given by the following expression:

(22)

We then calculate the Lagrangian perturbation to a surface element:

(23)

The approach used to obtain Eqs. (22) and (23) is essentially the same as that of Buta1979 and Townsend1997, but the effects of horizontal Lagrangian displacements are also included in the latter equation. We then calculate :

(24)

where . The Lagrangian perturbation to is calculated as follows:

(25)

where we have used the relation . In the spherical limit, the above expressions become

(26)
(27)
(28)
(29)

in full agreement with the expressions previously obtained by Heynderickx1994 for a non-rotating star with , provided one uses the relations:

(30)

3.2 Lagrangian perturbation to the effective gravity

Although the pulsation equations are given for a cylindrical rotation profile, we will relax this assumption in this section and allow for a general 2D rotation profile when establishing an expression for the Lagrangian perturbation to the effective gravity. Before giving the Lagrangian perturbation to the effective gravity, it is useful to recall various expressions for the unperturbed surface effective gravity:

(31)

The minus sign in the last equation comes from the fact that gravity is pointed inward, being the norm of . Furthermore, contrarily to what happens in the stellar interior, the derivative of vanishes at the surface since the surface is in pressure equilibrium. Any of the above expressions may be used to evaluate the surface effective gravity, although provides the most accurate numerical results. It is also worth noting that the above expression neglects any contributions from meridional circulation and viscous forces to the equilibrium model (e.g.  Rieutord2009).

The relative Lagrangian perturbation to the effective gravity is then given by:

(32)

where we have used the simplification . The quantity represents the vectorial Lagrangian perturbation to the effective gravity. It includes the Lagrangian perturbation to the gradient of the gravitational potential, and the acceleration of a particle tied to the surface, resulting from the oscillatory motions:

(33)

Rather than working with the above expression, it is more useful to introduce the equilibrium effective gravity by adding and subtracting :

(34)

After some simplifications based on Poisson’s equation and a lengthy derivation described in App. B, one obtains the following explicit expression for a general rotation profile:

(35)

It is worth noting that in the above expression, all of the terms involving a derivative of one of the perturbed quantities are divided by . The remaining terms involve no derivatives whatsoever. This characteristic is what one would expect for a quantity which is independent of the mapping.

If a cylindrical rotation profile is used instead, the terms in square brackets simplify to the following expressions, respectively:

(36)

In the non-rotating limit, Eq. (35) reduces to:

(37)

in full agreement with Dupret2002.

Given the complexity of Eq. (35), it is interesting to see if it can be approximated by a simpler expression. We consider the following approximation:

(38)

Figure 2 compares and for three different modes – two gravito-inertial modes and one p-mode. As can be seen in this figure, is a very good approximation for p- and g-modes with a sufficiently high frequency. It is only for low frequency g-modes that the difference between the two expressions becomes non-negligible.

Figure 2: Comparison between and for three modes. The upper row shows the meridional cross-section of the three modes, whereas the lower row compares the two calculations for the perturbation to the effective gravity at the stellar surface, between one of the poles and the equator . For the middle and right mode, the two curves nearly overlap, making it hard to distinguish them.

3.3 Lagrangian perturbation to the temperature

We then turn our attention to the relative Lagrangian perturbations of the effective temperature, . Given that the pulsation calculations are done in the adiabatic approximation, this quantity is not readily available. We will therefore content ourselves with the approximation . We note that this approximation is different from the approximation given that is latitude-dependant whereas is not in barotropic stellar models. As was pointed out in Sect. 2.5, a modified version of the mechanical boundary condition was needed in order to obtain a non-zero value for at the surface. One may then wonder if other approximations may be more suitable, such as applying the boundary condition high up in the stellar atmosphere and extracting at a deeper optical depth, or working with the Eulerian temperature fluctuation instead. The main difficulty with the first option is that the stellar models used in the present study do not include an atmosphere222The Kurucz atmospheres described in the following section are only used to calculate the emergent intensities as a function of and , and are not “joined” to the present models.. Hence, it is not clear at what depth (which is likely to be latitude-dependant), one should extract the Lagrangian temperature variations. Using the Eulerian temperature variations also has problems of its own. Indeed, at the surface, the Eulerian temperature variations are dominated by the the advection term, thereby making them an order of magnitude larger than the Lagrangian variations, and with the opposite sign for the most part. Such variations do not seem the most appropriate from a physical point of view, since the intensity variations are more likely to result from local temperature variations physically present in the fluid, i.e. Lagrangian variations. Hence, given the limitations of our models and pulsation calculations, using the Lagrangian temperature variations along with the modified mechanical boundary condition seems like the best choice. Nonetheless, according to Dupret2002 and Dupret2003, adiabatic calculations of are not reliable in the superficial layers and consequently lead to a poor approximation of . In order to obtain accurately, one would need to do a full non-adiabatic calculations including a detailed treatment of the atmosphere such as what is done in Dupret2002, using 2D rotating models in thermal equilibrium, such as what is being developed in the ESTER project (Rieutord2009; Espinosa2010).

The quantity can be deduced very simply from via the adiabatic relation:

(39)

where is the second adiabatic exponent, given by the expression .

3.4 Specific intensities

Contrarily to the non-rotating case, the effective temperature and gravity of the equilibrium model depend on the latitude. Consequently, it is necessary to find the appropriate set of specific intensities and their derivatives for each latitude. We therefore calculated a grid of Kurucz atmospheres with solar composition which spans the relevant temperature and gravity ranges. We used the ATLAS 9 code333See http://kurucz.harvard.edu (Kurucz1993) in a modified form so as to include the convective prescription of Canuto1996, known as CGM (for more details see Heiter2002 and Barban2003). Due to convergence problems, some of the grid points were missing and had to be interpolated from neighbouring points. The final grid is illustrated in Fig. 3 (upper left panel) as is as a function of for a set of models with rotation rates ranging from to of the critical break-up rotation rate.

For each grid point, we then approximated the dependence through a Claret type law (Claret2000) using a least-squares fit. This yielded a set of five coefficients from which it is possible to calculate and for any . Using this approach rather than dealing directly with the original values for yields better numerical results, both for interpolating and especially for taking its derivative.

We then successively applied cubic spline interpolations, first as a function of , then as a function of . This yielded b-spline coefficients from which it is possible to deduce the Claret coefficients and their derivatives for any value of and within the relevant range. From this representation of , we then calculated the Claret coefficients for , and as a function of latitude for each stellar model. The remaining three panels of Fig. 3 display these quantities as a function of and colatitude, , for a model at of the critical rotation rate. As is illustrated in these panels, both limb and gravity darkening are taken into account. The latitude-dependant Claret coefficients were then projected onto the spherical harmonic basis, truncated at 40 terms. This last step allows us to easily interpolate these coefficients onto a denser latitudinal grid during the visibility calculations.

Figure 3: (Upper left panel) Effective temperatures and gravities of the Kurucz model atmospheres (symbols) and the models (continuous lines; the lower longer lines correspond to more rapid rotation). The pluses represent initial Kurucz atmospheres whereas the diamonds represent missing points which were interpolated from the neighbouring atmosphere models. (Remaining panels) Specific intensities and its partial derivatives for a , model, as a function of colatitude and .

3.5 Visibility integrals

Once the above quantities are calculated over the stellar surface, the visibilities may be evaluated numerically. However, rather than using the time dependant form given in Eq. (18), it is more useful to extract the amplitude and phase of . To do so, we cast Eq. (18) into the following schematic form:

(40)

where and are real. This leads to:

(41)

where is the amplitude and the phase. Based on symmetry considerations with respect to the variable , the term cancels out in the adiabatic case, thereby leading to and or . If, however, one tries to simulate non-adiabatic effects by introducing a phase-lag in the temperature variations, then the term no longer cancels out.

The integrals and are evaluated numerically. The surface is discretised in the direction, using a Gauss-Legendre collocation grid with typically points, and in the direction using a uniform grid with a typical resolution of . The integration weights in the direction are deduced from a Gauss quadrature, whereas those in the direction are uniform. Figure 4, shows the differences which result from increasing either (left panel) or (middle panel) in a stellar model, rotating at of the critical rotation rate. As can be seen, these differences are several orders of magnitude smaller than the visibilities themselves.

At each point on the surface, the following condition is evaluated to determine whether or not the surface element is facing the observer and should be included in the integral:

(42)

where is given in Eq. (24). We note that this condition is only valid if no obstacles are present between the surface element and the observer, such as what could happen, for instance, in the more distorted SCF models with concavities at the poles (due to highly differential rotation) and an orientation which is not too far from equator-on. For such configurations, a more general condition based on a z-buffer or ray-tracing approach would be required.

One possible concern is that using a Gauss quadrature implicitly assumes a spectral decomposition of the function which is being integrated, which could lead to a Gibbs phenomena near the cutoff between the visible and hidden side of the star. However, given that the Gibbs phenomena is oscillatory, one can expect its integrated contribution do be small. Nonetheless, we carried out a test on the same model, using both Gauss quadrature and a trapezoidal integration. As can be seen in the right panel of Fig. 4, the two approaches yielded very similar results, showing that the Gibbs phenomena does not play an important role.

Figure 4: Visibility calculations and differences for a set of modes, obtained using two different resolutions, and (left panel), two different resolutions, and (middle panel), and two different integration methods in the theta direction (right panel). These modes were calculated in a model at , with an inclination angle . If is used instead, the differences are reduced.

3.6 Comparison with 1D case for non-rotating stars

It is also important to check that in the non-rotating case, visibility calculations based on the above 2D integrations agree with the results obtained with a simpler 1D integration, as is commonly used in other works (e.g.  Heynderickx1994). In order to obtain the 1D integrals, one starts with a given mode, which in this case will be proportional to a particular spherical harmonic, , re-expresses this spherical harmonic in terms of a new coordinate system with the -axis aligned with the line of sight via the following decomposition (Edmonds1960):

(43)

and calculates the integrated contribution from each of the . In this new coordinate system, the visible surface corresponds to a hemisphere. Consequently, all of the non-axisymmetric terms cancel out and only the term remains. This leads to the following expression:

(44)

where