A diagnostic diagram for γ Doradus variables and slowly pulsating B-type stars

Pulsating variables of γ Doradus type (γ Dor) and slowly pulsating B-type (SPB) stars are found on and near the main sequence with typical periods varying between one and several days, making them rather hard to detect from the ground. It is only with space missions such as CoRoT and Kepler that we became truly capable of determining their oscillation frequencies with enough precision to perform in-depth analyses. Here we present an efficient and easy-to-implement seismic tool, in which the frequency (ν) and the square root of the frequency difference ($ \sqrt{\Delta \nu} $) are plotted against each other as the abscissa and the ordinate, respectively. This allows us to immediately (1) perform mode identification; (2) estimate the average rotation rate and the characteristic period of gravity modes; and (3) recognise certain physical effects, including buoyancy glitches and avoided crossings. This diagnostic tool can only be applied to prograde sectoral g modes. To validate the tool presented here, we used stellar models and also applied it to three γ Dor (KIC 12066947, KIC 5608334 and KIC 4846809) and one SPB star (KIC 3459297), all observed with Kepler. Furthermore, we show that the rotation rates determined using this new tool are consistent with the results of previous studies.


Introduction
Gamma Doradus (γ Dor) stars are a class of pulsating variable stars located on and near the main sequence with spectral types from late-A to early-F. The corresponding mass range is between 1.4 and 2.0 M 1 . They show multi-periodic variability with amplitudes of 0.1 mag and periods of about one day. The average measured projected rotational velocity, v sin i, is ∼70−100 km s −1 , which indicates that the surface rotation rate is of the same order as the pulsation periods. This means that rotation significantly affects the low-frequency oscillation modes in these stars, making them perfect objects for the study of how rotation affects the internal structure and evolution of these stars.
There has been substantial progress in the understanding of γ Dor stars, thanks to space missions such as CoRoT (Baglin et al. 2006) and Kepler (Borucki et al. 2010), which delivered long and continuous observations with unprecedented precision of a vast number of stars. Contrary to ground-based data, space observations do not suffer from daily aliasing issues, which pose a serious problem when studying pulsating stars with periods close to one day. It is only thanks to multiple years of data from Kepler that oscillation periods in γ Dor stars have finally been resolved, leading to a true revolution in the study of these stars (e.g. Van Reeth et al. 2015). While Kepler data have revealed a significant number of oscillation modes in many cases, their exact identification is absolutely crucial for in-depth analyses.
Oscillation frequencies in rapid rotators are fundamentally different from those of slow rotators. For the latter, the rotation can be regarded as a small perturbation and so the average internal rotation rate can be estimated by measuring the rotational splittings in the frequency spectrum (Cowling & Newing 1949;Ledoux 1951). In the case of rapid rotators, the period versus period-spacing (P-∆P) diagram is known to provide useful information about the mode identification and the internal rotation rate (e.g. Bouabid et al. 2013;Ouazzani et al. 2017). The asymptotic formula of the period spacing (Ballot et al. 2012) serves as a useful theoretical basis for extracting information from the frequency spectrum. Nevertheless, the estimation of the rotation rates of these stars is not a simple task. Van Reeth et al. (2016), for example, use a large grid of evolutionary models with an assumed range of rotation rates to find the bestfitting parameters for their sample of stars, including the rotation rates. On the other hand, Christophe et al. (2018) propose a model-independent method based on the discrete Fourier transform. The basic idea is that the pulsation periods will be evenly spaced even when the star is rapidly rotating if the period coordinate is stretched appropriately. Li et al. (2019) developed a different model-independent method in the same theoretical framework as Van Reeth et al. (2016) and Christophe et al. (2018), which is described in detail in Sect. 2.1. Li et al. (2020) then applied the method to as many as 611 γ Dor stars in the Kepler field. The central problem of all the methods mentioned above is the parameter sweep in a multidimensional space, which is generally time-consuming and computationally intense.
In this paper, we show that the internal rotation rate of γ Dor and slowly pulsating B-type (SPB) stars can be estimated in a model-independent way by using a simple linear relation between the oscillation frequency and the square-root of the frequency difference. This method can only be applied to prograde sectoral g modes, which is fortunate as observations (e.g. Saio et al. 2018a;Li et al. 2019) show that γ Dor stars predominantly pulsate in these modes. We confine ourselves to single stars in this paper because tidal effects in binaries (or multiple systems) could complicate the picture significantly. Prograde sectoral g modes consist of the waves that propagate in the same longitudinal direction as the rotation, with large amplitudes confined to the equatorial region. Their main restoring force is buoyancy, while the Coriolis force plays an essential role in the equatorial confinement. These modes, which are found in γ Dor stars, are also excited in SPB stars (Waelkens 1991;Salmon et al. 2014). These stars are known as gravity-mode pulsators on the main-sequence with the mass between 3−9 M . The second most dominant component in the frequency spectrum of γ Dor stars is that of r modes (Papaloizou & Pringle 1978;Van Reeth et al. 2016;Saio et al. 2018b), which appear at lower frequencies than prograde sectoral g modes. Rossby (r) modes are mostly composed of toroidal waves (Rossby waves) in each spherical (or equipotential) layer, which are restored by the Coriolis force. The (phase) propagation of the modes is opposite to the rotation in the longitudinal direction (in the co-rotating frame), while the wave motions in adjacent layers are weakly coupled with each other through the influence of the buoyancy. Observations suggest that these modes are detected not only in γ Dor stars, but also in other various types of stars, including chemically peculiar stars and accreting white dwarfs (Saio 2018(Saio , 2019. The structure of this paper is as follows: the method is described in detail in Sect. 2. It is then applied to (a) synthetic model frequencies and (b) observed frequencies of four stars in Sect. 3. Section 4 is devoted to discussions and conclusion.

Theoretical framework
The present study focuses on eigenmodes of stellar oscillations, specifically, prograde sectoral g modes that have been detected in γ Dor and SPB stars. These modes have significantly longer periods than the dynamical time scale of the star and they are strongly affected by the Coriolis force. This is because their rotation periods are of the same order as the oscillation periods. These properties imply that the modes are composed of gravitoinertial waves, which are waves that are subject to the influence of both buoyancy and Coriolis forces.
The equations that govern the low-frequency eigenmodes of linear adiabatic oscillation in rapidly and uniformly 2 rotating stars can be formulated to a good level of precision in the framework of the traditional approximation of rotation (Eckart 1960). We neglect (1) the perturbation to the gravitational potential (the Cowling approximation after Cowling 1941), (2) the centrifugal deformation (the assumption of spherical symmetry of the equilibrium structure), and (3) the horizontal component of the rotation vector. The derived equations are separable in the spherical coordinates (r, θ, φ), with the direction of θ = 0 aligned with the rotation axis. While the φ parts of the eigenfunctions are described by the sinusoidal function, e imφ , where m is the integral index called the azimuthal order, the θ parts come from the solution of the Laplace tidal equation (the Hough function). On the other hand, the radial parts of the governing equations are of the same form as those of the spheroidal modes of non-rotating stars under the Cowling approximation. The only difference from the non-rotating case is that the term ( + 1), where is the spherical degree, is replaced with the eigenvalue λ , m of the Laplace tidal equation. Assuming further that the radial wavelength of the oscillation is much shorter than the local scale height of the equilibrium structure, we may apply the asymptotic analysis (e.g. Vandakurov 1967) to derive the condition for the eigenfrequency of the low-frequency modes as where the meanings of the symbols are given as follows: k r is the radial wave number of the oscillation; r 1 and r 2 represent the inner and outer boundaries of the region where the constituent waves are propagative; n is the radial (integral) order of modes; α is the phase shift introduced at the inner and outer boundaries of the propagative region. The expression for the radial wave number k r is given by Here N(r) is the Brunt-Väisälä frequency, which is a function of the radius r, and ν co is the cyclic frequency of the oscillation in the co-rotating frame. The eigenvalue λ , m (s) depends on not only and m, but also the spin parameter s, which is defined by The index (spherical degree) of λ , m (s) means that λ , m (s) → ( + 1) as s → 0, while the angular part of the eigenfunctions cannot be described by any single spherical harmonic for s 0. In other words, the of eigenmodes in rotating stars indicates the degree of spherical harmonics that specify the angular dependence only in the limit of no rotation. We note that ν co is related to the oscillation frequency in the inertial frame, ν, by In this paper, we adopt the convention that positive (negative) values of m correspond to prograde (retrograde) modes. Equation (1) implies that the oscillation period in the co-rotating frame is given by in which the parameter P 0 is defined by Equation (5) is in general an implicit expression because s depends on ν co . In the absence of rotation, the series of highorder gravity modes with a given spherical degree are equidistantly spaced by a constant period spacing of P 0 / √ ( + 1) A106, page 2 of 15 provided that the composition profiles have no steep gradient, and that the dependence of r 1 and r 2 on the frequency is negligible. Although the spectrum of the low-frequency oscillations in rapidly-rotating stars does not show constant period spacings, the parameter P 0 still characterises the structure of the spectrum through Eq. (5). In this sense, P 0 can be regarded as a characteristic period of the gravity-mode oscillations 3 .

Prograde sectoral g modes and the
We derive the relation between the frequency and the square root of the frequency difference for prograde sectoral g modes with a given m. Supposing that there is a list of (observed or synthetic) frequencies sorted in the ascending order with the radial orders n k (k = 1, 2, . . .) and the same m, we can take the difference of Eq. (5) multiplied by λ ,m between the (k + 1)th and kth frequencies to obtain in which ∆ k q generally means the difference in q between the (k+1)th and kth frequencies in the list, except for ∆ k n = −(n k+1 − n k ) (≥1). When we derive Eq. (7), we assume that r 1 , r 2 and α are independent of the frequency. For the prograde sectoral g modes, which have = m > 0, there is an approximate relation for large s, (e.g. Bildsten et al. 1996;Townsend 2003). The accuracy of Eq. (8) is depicted in Fig. 1. We may further introduce the linearisation relation, under the assumption of ∆ k ν ν k+ 1 2 , where ∆ k ν = ν k+1 − ν k and ν k+ 1 2 = (ν k + ν k+1 )/2, with ν k being the kth frequency in the list.
With the help of Eqs. (8) and (9), we can reduce Eq. (7) to In Eq. (10), ∆ k ν (the difference between the two adjacent frequencies) and ν k+ 1 2 (the average frequency of the two) are both observables that are computed from the same pair of mode frequencies, while m, ∆ k n, P 0 and ν rot are not known in advance. With estimates of m and ∆ k n, we can construct a diagram, based on only the two observables, in which the abscissa and the ordinate are given by ν k+ 1 2 /m and √ ∆ k ν/(m∆ k n), respectively. We call this the diagram of ν versus √ ∆ν, or simply the ν-√ ∆ν diagram. Although mode identification (i.e., the knowledge of ∆ k n) is not always straightforward, this tool can be used to constrain ∆ k n of these prograde modes as described in Sect. 3. According to Eq. (10), each pair of the two adjacent frequencies should correspond to a point of a common straight line in the diagram. The slope and the abscissa intercept (the intersection point with the horizontal line of √ ∆ k ν/(m∆ k n) = 0) of the linear fit are given by √ P 0 and ν rot , respectively. Using this result inversely, we can estimate ν rot and P 0 by plotting the modes in the diagram and fitting them with a straight line. The simplest procedure is to perform linear least-squares fitting under the assumption that the errors of each point are uncorrelated with each other. This assumption is acceptable for rough estimates. However, as shown in the Appendix A, ∆ k ν and ∆ k+1 ν are correlated with each other since they both depend on ν k+1 .

Revised method of estimating ν rot and P 0
The formula given by Eq. (10) can be improved by discarding the two approximations that are given by Eqs. (8) and (9). The revised formula can directly be obtained from Eq. (7) as in which The correction factor f k , which is a function of ν rot , is practically close to 1 because Eqs. (8) and (9) are generally good approximations. This property can be utilised to construct the following iterative procedure to obtain better estimates of ν rot and P 0 : 1. set f k to f (1) k = 1 and perform the least-squares fitting based on Eq. (11) to estimate ν rot = ν (1) rot and P 0 = P (1) 0 (cf. Appendix A); this step is called iteration 1.

Error estimates
The errors in ν rot and P 0 obtained by the procedure in Sect. 2.3 can originate from (1) the observational errors in the frequency measurement and (2) the systematic errors in Eq. (11), which arise from differential rotation, steep chemical composition gradient (discussed in Sect. 3.1), and deviation from the traditional approximation of rotation in the asymptotic limit. Because the relative errors in the frequencies of γ Dor and SPB stars in the Kepler data are typically of the order of 10 −5 or below, we regard that the dominant uncertainties come from (2). In order to estimate the systematic errors accurately, we need a detailed comparison between the frequencies computed under the traditional approximation of rotation in the asymptotic limit and those determined by two-dimensional computations that take full account of the effects of the rapid rotation. We leave this expensive comparison to future work, however, for a first approximation, we determine uncertainties in ν rot and P 0 from the residual sum of squares of the least-squares fitting. These estimates are based on the assumptions that (a) the uncertainties in the frequencies are not correlated and that (b) they follow a Gaussian distribution with the same standard deviation. For the sake of simplicity, we assume that there is no error correlation, but we do point out that the systematic errors in Eq. (11) are likely to be highly correlated. We set the default confidence level to 99% for the uncertainties in the estimates of ν rot and P 0 . More details are given in Appendix A.

Results
As described in Sect. 1, the ν-√ ∆ν diagram method was developed based on Kepler observations. It indicates that most of the observed frequencies in γ Dor stars can be identified as prograde sectoral g modes and r modes. Therefore, we first address these two types of modes in this section, while the other types are discussed in Sect. 3.3.2.
Since analysing γ Dor and SPB stars requires very similar techniques, in this section we choose to concentrate on the former, while the latter is discussed in Appendix B.

Tests based on models
We construct the ν-√ ∆ν diagrams and estimate ν rot and P 0 for the following three models of γ Dor stars (all having the same mass of 1.86 M ), which are described in detail in Christophe et al. (2018): (1) Model A (with smooth profiles of chemical composition); (2) Model B (with steep gradients of the chemical composition profiles just above the convective core); (3) Model A with differential rotation (DR) (ν rot = 15 µHz and 7 µHz in the convective core and the envelope, respectively). For each model, we used the same set of eigenfrequencies as Christophe et al. (2018), that is, those of prograde sectoral g modes with m = 1 and the g mode radial orders between −50 and −20. All mode frequencies were computed with the acor oscillation code (Ouazzani et al. 2012(Ouazzani et al. , 2015. For the sake of simplicity, we took only the spherically-symmetric part of the centrifugal deformation into account. In Figs. 2-4 we show the ν-√ ∆ν diagrams for these models. The corresponding estimates of ν rot and P 0 are summarised in Table 1. The ν-√ ∆ν diagrams in Figs. 2-4 clearly show that all the open circle symbols are well aligned and fit a straight line. This confirms that Eq. (10) is a good approximation in all of the models. In models (1) and (2), the spin parameter of the eigenmodes ranges between 1 and 3, for which Eq. (8) is accurate to only about 10-20% (see Fig. 1). The inaccuracy of Eq. (8) results in the deviation of f k in Eq. (12) from 1. This is the main reason for the differences between the dashed and solid lines in Figs Open circles describe the modes with an azimuthal order m = 1 and ∆ k n = 1, which implies consecutive radial orders. The blue dashed and black solid lines describe the output of iterations 1 and 5, respectively. The open square symbol shows the abscissa intercept with the solid line, which gives the converged estimate of the rotation frequency.   A106, page 4 of 15 Table 1. Estimates of ν rot and P 0 for the three theoretical models from the synthetic eigenfrequencies of the prograde sectoral g modes with the azimuthal order of m = 1. ( Notes. The process to estimate ν rot and P 0 converges after 5, 5, and 4 iterations for models (1), (2), and (3), respectively, with the requirement that the relative differences in ν rot and P 0 are both less than 10 −4 between two successive iterations. The errors that are given for the converged values correspond to the 99% confidence intervals. The results of Christophe et al. (2018) based on the discrete Fourier transform are also given and described in the table as C18. The true value of ν rot for Model A with differential rotation (DR) indicates the range of the rotation rate in the model. and 3 or, equivalently, the differences between the outputs of iteration 1 and the converged values for these two models (see Table 1).
The converged values for the rotation rate based on the iterative linear least-squares fitting described in Sect. 2.3 are generally consistent with the results based on the discrete Fourier transform (DFT) map developed by Christophe et al. (2018), as shown in Table 1. Although the differences between the converged estimates and the true values are generally small, the true values of ν rot and P 0 are found outside the given confidence intervals for all cases, which indicates that the unaccounted systematic errors of Eq. (11) can be significantly larger than the uncertainties derived here.
As an example of the unaccounted error sources, we may list the steep gradients in the chemical composition profiles just above the convective core in model (2). This cannot be treated by the standard asymptotic analysis because it assumes that the wavelength of the oscillation is much shorter than the scale height of the equilibrium structure. The corresponding variation in the Brunt-Väisälä frequency (buoyancy glitch) is known to cause a wavy structure in the period spacing (e.g. Miglio et al. 2008;Bouabid et al. 2013;Kurtz et al. 2014). Such an oscillatory behaviour in ∆P caused by a glitch can also be found in √ ∆ν as shown in Fig. 3. The case of model (3) (see Fig. 4) can be used to check what happens if the method is applied to the case of differential rotation, although the uniform rotation is assumed in the (standard) traditional approximation of rotation. Broadly speaking, we find no critical difference from Fig. 2 because most of the points are aligned. This is, at first glance, consistent with the argument in Appendix E, which claims that the effect of weak differential rotation (in the radial direction) can be taken into account only by replacing the uniform rotation frequency by its average weighted by N/r over the propagation region. However, the weak differential rotation cannot be assumed for model (3) since the variation between 7 and 15 µHz is large. This means that for a strongly differentially rotating star, we cannot simply interpret the derived ν rot , even though the value of 11.80 ± 0.04 µHz for model (3) is clearly within the range of values in the model.
It is remarkable that our method allows us to highlight the two avoided crossings in Fig. 4 at a frequency of ∼15.5 µHz and ∼17.5 µHz, respectively. It should be noted that methods using the traditional approximation framework cannot account for this phenomenon. This is shown by the filled circles deviating from the solid line (base line). Since these points cannot be interpreted in the framework developed in Sect. 2, we simply discard them in the fitting to estimate ν rot and P 0 . Lee & Saio (1989) demonstrated that an avoided crossing can be caused by the horizontal component of the rotation vector, which is neglected in the traditional approximation of rotation. However, there may be other reasons too. In principle, it can occur due to any neglected terms that make the governing equations variable-inseparable into the radial and angular parts. In the case of differential rotation, those that include the gradient of the rotation frequency can also induce the avoided crossing.
Another important difference in model (3) is that some modes are missing, meaning that not all modes are of consecutive radial order. More specifically, the open diamonds aligned with the blue dotted line labelled "× √ 2" in Fig. 4 correspond to the jumps in the radial order n with ∆ k n = 2. This demonstrates that this method allows us to not only determine ν rot and P 0 but also ∆ k n. When a point is found on the line that has the same abscissa intercept as the base line but a slope steeper than the base line by factor of √ j (with j = 2, 3, . . .), we may interpret this as a missing j−1 radial order(s) with ∆ k n = j. With the value of ∆ k n, we can calculate the left-hand side of Eq. (10), which can be used in the fitting. This correction factor of 1/ √ ∆ k n is indicated by each of the vertical cyan lines in Fig. 4. We only use the open circles in Fig. 4 to estimate ν rot and P 0 .
We apply the same method and construct the ν-√ ∆ν diagram for r modes of model (1) (see Fig. 5). In sharp contrast to Fig. 2, these modes are not distributed along a straight line with a positive slope. This very clear difference can be used to distinguish r modes from prograde sectoral g modes. Thus, the diagram is also useful for mode identification in this respect.

KIC 12066947
The first γ Dor star analysed here is KIC 12066947. The observed frequencies and the Lomb-Scargle periodogram are provided in Table A Here we used the frequencies between 27 µHz and 35 µHz which were already identified to be prograde sectoral g modes with m = 1.
University ECHO (Extraction of CoHerent Oscillations) pipeline described in Antoci et al. (2019). Analyses show that the extracted frequencies are separated in two clusters, Cluster A between 20 and 24 µHz and Cluster B between 27 and 35 µHz.
The ν- Fig. 6 is constructed using 14 frequencies in Cluster B, which could be identified as m = 1 modes. It can be seen that except for one peak, all modes follow the three lines that are displayed. The solid line is the base line, while the two dotted lines labelled "× √ 3" and "× √ 6" indicate slopes that are larger than the base line by factors √ 3 and √ 6, respectively. This is consistent with the interpretation that these frequencies are those of the prograde sectoral g modes. The outlier at ν k+ 1 2 , √ ∆ k ν = (28.4, 1.08) is discarded when we perform the least-squares fitting, which converges after three iterations. Table 2 provides the converged estimates of ν rot and P 0 with their 99% and 1σ confidence intervals. These results are consistent with those of Van Reeth et al. (2016) and Christophe et al. (2018) within 1σ, while those of Li et al. (2019) are within 99% confidence intervals. We note that the uncertainties of Li et al. (2019) are smaller than the others because they include not only prograde sectoral g modes but also r modes in their analysis. The ν-√ ∆ν diagram constructed using the frequencies in Cluster A is shown in Fig. 7. It is clear that we fail to observe the modes showing the expected pattern for prograde sectoral g modes. We therefore conclude that these peaks are those of other modes, which is consistent with Van Reeth et al. (2016) who identify these as r modes (cf. Fig. 5).

KIC 5608334
Saio et al. (2018a) study KIC 5608334, which shows remarkable frequency groups (fg) in the amplitude spectrum (cf. Fig. 3 in the paper). We adopted their mode identification of the four frequency groups, fg1-fg4, as the prograde sectoral g modes with m = 1−4, respectively. From 192 frequencies that are determined with the ECHO pipeline of Antoci et al. (2019) and listed in Table A.1 of Saio et al. (2018a), we first excluded those with amplitudes smaller than 5 ppm (all of which belong to fg3 and fg4). We chose to do so because the small peaks show significant outliers in the ν-√ ∆ν diagram, which may imply that many of them do not correspond to prograde sectoral g modes or additional effects such as glitches are present. The physical origin of them should be examined in a separate study. In addition, if multiple frequencies are found in a very narrow frequency range for all of fg1-fg4 (cf. Fig. 4  There are a few possible physical reasons for the discarded points, which are indicated by filled circles in Figs. 8 and 9: they might originate from an avoided crossing (as in Fig. 4) or from a contamination of other types of frequencies (e.g. those of tesseral modes or combination frequencies) than those of the prograde sectoral g modes. These outliers are caused by the frequencies with small amplitudes: they disappear if we exclude those with amplitudes smaller than 100 and 20 ppm in fg1 and fg2, respectively. We observe a small deviation of some of the open circles from the base lines in Figs. 8-11. It is possible that they also originate from an avoided crossing or the influence of the nonlinear mode interaction, as discussed in Saio et al. (2018a). Although the estimates of ν rot and P 0 for fg1-fg4 are completely consistent with each other within the errors, we cannot conclude uniform rotation only from these results because the theoretical model of Saio et al. (2018a) indicates that there is little difference in the propagation cavity of the modes for the different    Table 3 as non-significant. Saio et al. (2018a) also note an isolated peak at 2.2397 d −1 (25.922 µHz) in the frequency spectrum, which could be due to rotational modulation.
If this is indeed the case, the derived surface rotation rate agrees well with the average internal rotation rate in Table 3 (cf. Appendix E). Li et al. (2020) also analyse this star, but based only on peaks belonging to fg1. They obtain ν rot = 26.0±0.1 µHz (2.25 ± 0.01 d −1 ) and P 0 = (4.4 ± 0.2) × 10 3 s. These are completely consistent with our results in Table 3.

Slow rotation
While the fundamental approximation of the present analysis, √ λ ≈ m (cf. Eq. (8)), can be justified only when the spin parameter satisfies s 1, we show in Appendix C that in principle the ν-√ ∆ν diagram method may also work for s 1. However, we should stress that it is essential not to use the approximation of √ λ ≈ m in the case of slow rotation, but to utilise the realistic values shown in Fig. 1.

Mode identification
As described above, the ν-√ ∆ν diagram allows us to distinguish prograde sectoral g modes from r modes (cf. Sect. 3.1), however other types of modes may complicate the mode identification process. In this section, we discuss possible issues that may influence mode identification.
Retrograde modes with = 1 and m = −1 behave differently from the prograde modes when plotted in both, the ν-√ ∆ν as well as the P-∆P diagrams (cf. Bouabid et al. 2013). More specifically, the sign of the slope is different: negative for retrograde and positive for prograde modes in the ν-√ ∆ν diagram. One relevant question may refer to how we discriminate prograde sectoral g modes with different azimuthal order m. We believe that this is not a real problem. As in the case of KIC 5608334 (Saio et al. 2018a), these modes with different m usually appear in different groups in the frequency spectrum, with smaller amplitudes for the group with higher frequencies.
Since the effect of the geometrical cancellation is more significant for the modes with larger m, we may identify those groups as the modes with m = 1, 2, . . . in the order of frequency. In particular, if there is only one group observed, and the ν-√ ∆ν diagram confirms that it consists of the prograde sectoral g modes, it is quite likely that the group is that of m = 1 modes.
Some possible confusion can be caused by the zonal modes (m = 0), which are claimed to be observed in rare cases (cf. Van Reeth et al. 2016;Li et al. 2020). In Fig. 12 we show the = 1 zonal modes of model (1) in the ν-√ ∆ν diagram. We clearly observe very similar alignment of the points on a straight line to that in Fig. 2. If we erroneously identify these as prograde sectoral g modes with m = 1, we obtain the rotation rate of ν rot = 3.6 µHz (0.31 d −1 ), which is quite different from the true value of 7 µHz (0.60 d −1 ). This example shows that it is dangerous to apply the method blindly. We therefore suggest for the mode identification to use the information about the amplitude in addition to the ν-  Fig. 15 of the paper), in which the shorter-period (higher-frequency) series have significantly larger amplitudes. Both series cannot be prograde sectoral g modes, because frequencies in the higher-frequency group (with higher azimuthal order, m) are more likely to have smaller amplitudes due to the more severe effect of the geometrical cancellation. We therefore argue that the higher-frequency (shorter-  Other types of modes, such as = 2 and m = ±1, can also introduce confusion similar to that caused by the zonal modes. They appear as a series with a negative slope in the P-∆P diagram, which translates into one with a positive slope in the ν-√ ∆ν diagram. However, they are unlikely to have much larger amplitude than prograde sectoral g modes when the stellar rotation is fairly high and are, apparently, identified in only few stars (Li et al. 2020). Such exceptional cases should carefully be investigated in a separate study.
Having considered all of the above, we may propose a general procedure for mode identification as follows: 1. The first step is to identify frequency groups, each composing a cluster in the frequency spectrum of a given star. If this is not possible, it may be the case that the star rotates slowly and that frequency splittings may be identified.
2. If one or more groups are present, mode identification is required. From a practical point of view, it seems to be a good initial guess that the group with the largest amplitude corresponds to the prograde sectoral g modes of m = 1. For these modes we expect to observe larger amplitudes because the geometrical cancellation effect is significantly smaller in comparison to other types of modes. We may find an exception for r modes but only if the star rotates rapidly, or more precisely if the spin parameter satisfies s 1 (Saio et al. 2018a).
3. If the modes of the first group do not align in the ν-√ ∆ν diagram and show a positive slope consistent with prograde sectoral g modes, then the frequency group with the second-largest amplitudes should be examined.
4. This process should be repeated until prograde sectoral g modes with m = 1 are found (or all of the frequency groups have been checked). If there are two (or more) frequency groups with similar amplitudes at each step, we may assume that the one with the highest frequency originates from the prograde sectoral g modes with m = 1.
5. Once the identification of the first group is successful (i.e. m = 1), one may look for m = 2, 3, . . . modes that should be located around the frequencies of the first group multiplied by m.
We regard this procedure as a general guideline that is applicable to the majority of the stars. The further details of the mode identification should be examined carefully based on a large sample of stars.

Discussion and conclusion
We propose the ν-√ ∆ν diagram as an asteroseismic tool to estimate the average rotation rate (ν rot ) and the characteristic period of the gravity-mode oscillations (P 0 ) of γ Dor and SPB stars from the frequencies of the prograde sectoral g modes. Our method is endorsed by the fact that the prograde sectoral g modes are most commonly observed in (rapidly rotating) γ Dor stars and SPB stars (e.g. Van Reeth et al. 2015;Pápics et al. 2017). Although the method is based on the same asymptotic formula of the low-frequency oscillation under the traditional approximation of rotation as in previous works (Van Reeth et al. 2016;Christophe et al. 2018;Li et al. 2019Li et al. , 2020, the biggest advantage here is in its efficiency and simplicity. It is easy to get rough estimates of ν rot and P 0 by plotting the diagram and applying a linear fit, which the majority of the modes follow ideally. The estimates of ν rot and P 0 can be derived from the abscissa intercept and (the square of) the slope of the line, respectively. If no clear alignment is visible, it is unlikely that the observed frequencies correspond to prograde sectoral g modes. The diagram thus also allows us to perform mode identification. In addition, this tool can also be used to discuss the physics of oscillation modes including the buoyancy glitch (Fig. 3) and the avoided crossing (Fig. 4). The first rough estimates of ν rot and P 0 can be improved by just a linear least-squares fitting with a few iterations. It should be stressed that the method is free from any time-intensive parameter search in two-(or more) dimensional space. There is no need to constrain the ranges of any parameters in advance or give any initial guesses. These properties are important, especially when analysing a large sample of stars.
The method proposed in this paper depends on the eigenvalue of the Laplace tidal equation of the prograde sectoral g modes being approximately equal to the square of the azimuthal order, m 2 , if the spin parameter s is much larger than 1 (cf. Fig. 1). This implies that the mode frequencies are nearly independent of ν rot in the co-rotating frame (cf. Eq. (5)). This is because the dispersion relation of the constituent waves (internal equatorial Kelvin waves) is independent of ν rot when s 1 and also because they propagate only along the equator in the same limit. More details about the physical properties of the modes, which may be called Kelvin g modes, are discussed in Appendix D. While the mode frequencies are insensitive to ν rot in the co-rotating frame, we may infer ν rot from those in the inertial frame through the advection effect induced by rotation (Doppler shift), which is represented by the second term on the right-hand side of Eq. (4). The ν-√ ∆ν diagram allows us to represent this effect graphically.
While it is standard to assume that the rotation rate is constant in the traditional approximation of rotation, this is not necessary in the asymptotic limit. In the case where the rotation rate weakly depends on the radius, the obtained estimate of ν rot should be interpreted as the average of the internal rotation rate weighted by N/r (the Brunt-Väisälä frequency divided by the radius), as is proposed by Ouazzani et al. (2019) and justified in Appendix E.
We confirm that the estimate of ν rot is highly correlated with that of P 0 (cf. Van Reeth et al. 2016;Christophe et al. 2018): the correlation coefficient a, which is defined in Eq. (A.29), satisfies a ≥ 0.98 in most cases presented in this paper. Such high values can be understood based on the ν-√ ∆ν diagram. Supposing we have obtained the best-fit line on the diagram, if we slightly increase (or decrease) the abscissa intercept (ν rot ) from the bestfit value and perform a linear fit under this constraint, the slope of the fitted line ( √ P 0 ) is necessarily steeper (or shallower) than the best-fit line because all the points in the diagram are located in the upper right side of the abscissa intercept (the domain of ν/m > ν rot and √ ∆ν > 0). The reason for the high correlation coefficient is that we take only prograde sectoral g modes into account. We can, in principle, reduce the correlation by including other types of modes, if available (cf. Li et al. 2019). In this context, a natural extension of the present method would be to include the other types of modes by deriving similar formulae to Eq. (11).
We adopt conservative estimates of the formal uncertainties in the present analysis, corresponding to 99% confidence intervals calculated based on the goodness of the fit. Nevertheless, Table 1 shows evidence that these values can significantly underestimate the true uncertainties. This is probably a common problem in any method based on the asymptotic frequency formula under the traditional approximation of rotation and it should be examined in detail in a future work. In order to evaluate the total uncertainty, it is necessary to improve Eq. (11), which is the fundamental relation of the present work. This can be done by taking into account additional physical effects, including centrifugal deformation, correction to the traditional approximation of rotation, higher-order terms in the asymptotic frequency formula, and the possible steep variation of the chemical composition profiles just outside the convective core.

A.1. Simplifying assumptions
Here we describe the method of estimating ν rot and P 0 with their uncertainties based on Eq. (11) in detail. We paid special attention to the treatment of the error correlation. The following three assumptions are made for the sake of simplicity: (1) the errors { k } of the mode frequencies {ν k } follow a Gaussian distribution with the same (but unknown) standard deviation, σ, and those of different modes are not correlated with each other; supposing that the number of the mode frequencies is K + 1, and that indicates the (K + 1)-element vector whose kth element is equal to k , this assumption implies in which u T , E (a) and I k generally mean the transpose of u, the expected value of a and the identity matrix of size k, respectively; (2) because the relation ν k+ 1 2 |∆ν k | holds for any data sets we use in this paper, the frequency errors affect only the vertical position of the modes (symbols) in the ν-√ ∆ν diagram; (3) since f k , which is defined by Eq. (12), is generally close to 1, the frequency errors have negligible impact on f k .

A.2. Formulation in the vector form
We assume for the moment that the K + 1 mode frequencies have consecutive radial orders and the same azimuthal order m, so that there are K relations in the form of Eq. (11). This assumption is not essential (cf. the second paragraph of Appendix A.3). Neglecting the second-and higher-order terms in { k }, we may rewrite those relations as y = Xθ + e, (A.2) in which e represents the error term caused by { k }. The definitions of the symbols in Eq. (A.2) are given as follows: the kth element of the K-element vector y is given by the 2-element vector θ is provided by the kth element of the K-element vector e is given by in which We note that e can be expressed by respectively.

A.3. Error correlation
The variance-covariance matrix of the error vector e can be calculated from Eqs. (A.1) and (A.8) as in which C 0 = HQQ T H. We separate C 0 from C because C 0 can be calculated from the mode frequencies, whereas σ is unknown.
The fact that C is a tridiagonal matrix means that any two successive elements of e are correlated with each other. The inverse of C can be expressed by in which the (i, j) element of C −1 0 is given by It has so far been assumed that the K relations in the form of Eq. (11) are all included in the analysis. We may have neglected some of them because they produced outliers in the ν-√ ∆ν diagram. The possible reasons include avoided crossings, which cannot be described by the traditional approximation of rotation, and contamination of other types of modes than the prograde sectoral g modes. In this case, the expression of C (and its inverse) should be modified accordingly. For example, if we discard the kth element of Eq. (A.2) for 1 < k < K, the remaining relations are separated into two groups, one with the 1, . . ., (k − 1)th elements and the other with the (k + 1), . . ., Kth elements. Since e k−1 is not correlated with e k+1 (cf. Eq. (A.6)), there is no error correlation between the two groups. The modified expression of C −1 0 is composed of two diagonal blocks both in the form of Eq. (A.13) with K replaced by the block size, which is equal to k − 1 and K − k for the first and second blocks, respectively. This treatment can easily be generalised as follows: supposing that the relations are separated into J different groups with no error correlation, C −1 0 should be composed of J diagonal blocks in the form of Eq. (A.13). Based on this, it is possible to include the mode frequencies with different azimuthal order m in the analysis because the errors of any two relations with different m in the form of Eq. (11) are not correlated with each other. The different values of m correspond to different diagonal blocks of C −1 0 , and each can generally consist of multiple diagonal blocks in the form of Eq. (A.13).

A.4. Linear least-squares fitting
The problem of the linear least-squares fitting in the presence of error correlation is discussed in the literature (e.g. Eadie et al. 1971;Gough & Sekii 2002). The key idea is to multiply Eq. (A.2) by an appropriate K × K regular matrix from the left side so that any two elements of the resultant error vector are not correlated with each other. Once the error correlation disappears, we can adopt the standard procedure of the least-squares fitting. The procedure can be rewritten in terms of the original variables. In fact, we should find such θ that minimises (y − Xθ) T C −1 0 (y − Xθ) (A.14) for given y, X and C −1 0 . The solution is given bŷ where A is a 2 × 2 matrix defined by Vectorθ is the best linear unbiased estimator of θ, whose variance-covariance matrix is provided by On the other hand, the residual sum of squares divided by K − 2 yields the minimum variance unbiased estimator of σ 2 , The corresponding estimator of ν rot is naturally found to bê Neglecting the second-and higher-order terms in { k }, we can show thatν rot is unbiased. Its variance is accordingly calculated as The uncertainty inν rot can be estimated based on the fact that follows Student's t-distribution with K − 2 degrees of freedom.
Here, s (ν rot ) stands for the standard error ofν rot defined by The q% confidence interval of ν rot is provided bŷ in which t 0 is such a value that the probability of |t| ≤ t 0 is equal to q/100. Alternatively, t 0 can be set (for K > 4) to which corresponds to one standard deviation.
The way of estimating P 0 is analogous. All we need to do is replace the expressions for the estimator, its variance and the standard error witĥ respectively. Finally, the covariance and the correlation coefficient betweenν rot andP 0 can be computed as and a = cov ν rot ,P 0 var (ν rot ) var P 0 , (A.29) respectively.

Appendix B: Application to SPB stars
In this section, we apply the ν-√ ∆ν diagram method to SPB stars. First we validate the method based on a typical evolutionary model, and then apply it to KIC 3459297.

B.1. Validation based on a model
As in the case of Model A for γ Dor stars (cf. Sect. 3.1), we use the stellar evolution code cles (Scuflaire et al. 2008) to construct a 4 M model of an SPB star. The model parameters are summarised in Table B.1. The details on how this model is calculated are the same as in Model A (cf. Christophe et al. 2018), however we note that turbulent diffusion is included to approximate the effect of rotationally-induced mixing, which smoothes the chemical composition profiles outside the convective core. The eigenfrequencies of the prograde sectoral g modes with m = 1 are computed using the acor oscillation code (Ouazzani et al. 2012(Ouazzani et al. , 2015 with a uniform rotation frequency of 7 µHz (0.60 d −1 ). The range of radial orders is between −50 and −18. The results are shown in Fig. B.1.
We confirm that also in the case of an SPB star all modes follow a nearly perfect linear trend, which is similar to model (1) shown in Fig. 2. Our estimates for the rotational frequency ν rot and the characteristic period P 0 are provided in Table B.2. Strictly speaking, these estimates are slightly smaller than the true values by 0.9% and 4% for ν rot and P 0 , respectively, which is similar to the case of Model A listed in Table 1. We find that this is an intrinsic limitation of the method, which relies on the asymptotic formula under the traditional approximation of rotation. Nevertheless, since the two parameters are recovered reasonably well, we consider that the method can be applied to not only γ Dor stars, but SPB stars as well.  Notes. The columns describe the following parameters: T eff is the effective temperature; X c is the hydrogen mass fraction at the centre; Y ini and Z ini are the initial mass fractions of helium and heavy elements, respectively; α MLT is the mixing length parameter of convection; D t is the coefficient of turbulent diffusion. We use the frequencies of KIC 3459297 that are listed in Table A.2 of Christophe et al. (2018). While those frequencies distribute in three separate groups (cf. Fig. 8 of Christophe et al. 2018), we ignore the lowest-frequency group below 1.5 µHz because it probably originates from instrumental effects or the artefacts introduced during data analysis. The group between 9 and 14 µHz and that between 23 and 31 µHz are referred to as fg1 and fg2, with 24 and 8 frequencies, respectively. We identify the frequencies in fg1 and fg2 as those of prograde sectoral g modes with m = 1 and m = 2, respectively. The ν-√ ∆ν diagram is shown in Fig. B.2. As in the case of KIC 5608334, which is discussed in Sect. 3.2.2, we discard four frequencies (10.219, 10.354, 11.233 and 24.903 µHz) from the list because there is another frequency very close to each of them with higher amplitude. We also exclude the two red points in Fig. B.2 in the ν rot and P 0 fitting because they do not follow the same linear trend as the other modes. These peaks might indicate the presence of a buoyancy glitch (cf. Fig. 3) or an avoided crossing (cf. Fig. 4). The red and blue circles correspond to the frequencies in the first group (with m = 1) and the second group (with m = 2), respectively. The dashed and solid lines are the fitting results of all the modes in iterations 1 and 5, respectively. The two red filled circles are excluded from the fit.  frequencies of prograde sectoral g modes with m = 1 and radial orders between −31 and −17. The corresponding range of the spin parameter is between 0.70 and 0.37. Figure C.1 shows the ν-√ ∆ν diagram for this slowlyrotating model (Model As). Although Eq. (8) is not very accurate for these modes (cf. Fig. 1), we still observe that all modes follow the linear trend in the diagram. We thus cannot find any qualitative difference from the fast-rotating case in Fig. 2. Table C.1 provides the estimates for ν rot and P 0 , which are larger than the true values by only 3% and 0.5%, respectively. Comparing these results with the case of fast rotation in Table 1, we find the following: (1) larger relative uncertainties for the case of slow rotation by about an order of magnitude for ν rot and multiple factors for P 0 ; (2) larger differences between the outputs of the initial and last iterations in spite of the similar number of iterations necessary for the convergence (5 and 6 for model (1) and Model As, respectively). The reason for point (1) is the difference in the number of frequencies that are considered in the analyses. If we extend, to include all the modes between n = −17 and n = −52, the uncertainties of ν rot and P 0 are reduced to values similar to model (1). Point (2) is simply because the deviation from Eq. (8) is more significant. We may still regard that this approximation remains reasonably good even in the limit of slow rotation. In fact, the limiting value of √ λ = 4 √ m (m + 1) for the non-rotating case is different from the approximation of √ λ = √ m by only about 20% for m = 1. Therefore, the factor f k given by Eq. (12) changes during iterations by maximum this amount. This property helps the iterative procedure to be robust. However, in the case of slow rotation it is crucial to use the accurate numerical values of λ, which are shown in Fig. 1, rather than its approximate form of m 2 for s 1. A failure to do so will impact the accuracy of ν rot and P 0 .
The method also works for a slow rotation because the asymptotic formula for the oscillation frequencies under the traditional approximation of rotation (cf. Eq. (5)) can be justified even in the limit of |s| 1. In fact, it can generally be shown that, for |s| 1, λ is given by Higgins 1968), from which we can derive 2) is the standard formula of the rotational splitting for the high-order g modes in the case of very slow rotation (cf. Cowling & Newing 1949;Ledoux 1951). Of course, it is more practical in that case to use the rotational splittings rather than the ν-√ ∆ν method to determine ν rot , while P 0 can be measured directly from the period spacing.
In conclusion, we find no serious problem in applying the method of the ν- Prograde sectoral g modes are composed of internal equatorial Kelvin waves when the spin parameter s is much larger than 1 (cf. Townsend 2003). This can be understood from the common physical properties between the Kelvin waves and the prograde sectoral g modes, which are described in meteorology and oceanography literature (e.g. Gill 1982). These waves are trapped in the equatorial region, with their displacement predominantly parallel to the equator. Their horizontal propagation is in the same longitudinal direction as the rotation, while they also propagate in the radial direction because of the buoyancy force. Unlike the other gravito-inertial waves, they are restored only by the buoyancy force, whereas the Coriolis force does not influence the wave propagation except for controlling the range of the equatorial confinement by contributing to the force balance (geostrophic balance) in the latitudinal direction. This explains why the dispersion relation does not depend on the rotation rate. Having understood the preceding properties, we may call the prograde sectoral g modes: "Kelvin g modes". We note that Kelvin modes are sometimes referred to as f modes, oscillation modes of non-rotating stars that consist of surface gravity waves (e.g. Chandrasekhar & Lebovitz 1964;Cox 1980). There should be no confusion between the two types of Kelvin modes of totally different physical characters. The Kelvin waves in the rotating case originate from Thomson (1880), while Thomson (1863) studies the oscillation modes of a liquid sphere (a homogeneous incompressible spherical body), which correspond to the surface gravity modes.

D.2. Variation in λ
When there is no rotation (with the spin parameter of s = 0), the Laplace tidal equation is reduced to the standard associated Legendre differential equation, so that we obtain the eigenvalue of λ m, m = m (m + 1) for the prograde sectoral g modes. As demonstrated in the previous studies (e.g. Bildsten et al. 1996;Lee & Saio 1997), λ m, m monotonically decreases as s becomes large, and is asymptotically equal to m 2 as s → ∞.
We may interpret the transition from λ m, m = m (m + 1) for s = 0 to λ m, m = m 2 for s → ∞ as being caused by the change in the geometry. When s = 0, the system is totally spherically symmetric, so that the horizontal propagation of the waves can reach all over the (two-dimensional) sphere. From a mathematical point of view, the eigenvalue, m (m + 1), is that of the angular part of the Laplace operator, whereas it is physically determined by the quantisation condition that the constituent waves are in phase before and after travelling around the sphere (e.g. Gough 1993). On the other hand, the propagation region of the waves degenerates towards the equator (one-dimensional sphere) in the limit of s → ∞. In addition, the geostrophic balance forces the latitudinal (θ) component of the displacement to be zero, leaving only the longitudinal (φ) component. Therefore, only the longitudinal waves along the equator are allowed. The eigenvalue of λ m, m = m 2 can again be obtained by the quantisation condition, which requires the same phase of the waves before and after going around the equator. We can formally show that the Laplace tidal equation for p (the Eulerian perturbation to the pressure) is reduced to −m 2 p = −λp , which means ∂ 2 p ∂φ 2 = −λp , (D.1) at the equator (θ = π/2) under the condition of the geostrophic balance. Equation (D.1) is equivalent to the spatial part of the wave equation at the equator.

D.3. Contribution of the spheroidal and toroidal components
We may discuss the decrease in λ as s increases from a different point of view than the one in Appendix D.2. We first note that λ can be regarded as a measure of the divergence of the horizontal components of the displacement, ∇ h · ξ h = 1 sin θ ∂ ∂θ sin θ ξ θ + 1 sin 2 θ ∂ξ φ ∂φ , (D.2) in which ξ θ and ξ φ are the θ and φ components of the horizontal displacement ξ h , respectively. Under the traditional approximation of rotation, the angular dependence of the eigenfunctions is described by the Hough functions, Θ (cos θ), which are the solutions of the Laplace tidal equation. Based on this, the following can be shown: where Ξ h indicates the (common) radial dependence of ξ θ and ξ φ . The denominator on the right-hand side of Eq. (D.3) is just a normalisation factor. Next, we may suppose to express ξ h by a sum of two parts, the spheroidal and toroidal components, which are linear combinations of the terms in the form of respectively. Here, Y , m (θ, φ) is a spherical harmonic of degree and azimuthal order m, while e r , e θ and e φ indicate the unit vectors in the r, θ and φ directions, respectively. Given that ξ h is expanded into the series, we may regard that λ is fixed by an average of the contribution of each term in the both components, specified by Eqs. (D.4) and (D.5). In the absence of rotation, the eigenfunction, ξ h , is contributed by only one term of the spheroidal component, which satisfies (D.6) We can thus understand λ = ( + 1) in this case. Although the toroidal component also contributes to ξ h when s 0, ∇ h · ξ h has no sensitivity to the component (divergence free). From these considerations, we may interpret that the gradual decrease in λ as s increases is caused by the decaying contribution of the spheroidal components to ξ h , which implies the growing contribution of the toroidal components.