A rapidly evolving high-amplitude $\delta$ Scuti star crossing the Hertzsprung Gap

In this work, we report the discovery of the rapidly evolving high-amplitude $\delta$ Scuti star KIC 6382916 (J19480292+4146558) which is crossing the Hertzsprung gap. According to the analysis of the archival data, we find three independent pulsation modes, whose amplitudes and frequencies vary distinctly in 4 years. The linear period variation rates of the first two modes are about 3-4 times larger than the best seismic model constructed by the standard evolution theory, while that of the third one is about 8 times larger than the first two modes. What is more interesting is that, almost all the combinations of the third mode have frequency peaks $0.0815\ \mathrm{c\ d}^{-1}$ away from them in the frequency domain. A framework is proposed to interpret the markedly large frequency and amplitude variation rates of the third mode, in which we employ a new pulsation mode (resonating integration mode) generated by the resonance between a radial p-mode and a nonradial mixed mode. Moreover, a global analysis of the interactions between the three independent pulsation modes and their harmonics/combinations is performed based on the interaction diagrams of their amplitudes and phases, which would be a useful tool for the future asteroseismology research.


INTRODUCTION
People cannot clearly witness the stellar evolution process of a single star obviously in most cases because of its extremely secular time-scale, except for some special time nodes in it (such as the supernova explosion (Arnett et al. 1989)). However, in some specific evolutionary phases, we have the chance to witness such processes gradually on human timescales. When a star is leaving from the main sequence, the hydrogen nuclei fusion in its core is gradually transferring into the shell. In the Hertzsprung-Russell diagram, its evolutionary phase falls into the Hertzsprung gap, which is one of the most rapidly evolving phases in the life of a star (Kippenhahn et al. 2012).
δ Scuti stars are a class of short-period pulsating variable stars with periods between 15 minutes and 8 hours and spectral classes of A-F, which locate them on the main-sequence (MS) or post-MS evolutionary phases at the bottom of the classical Cepheid instability strip and excited by the κ mechanism (Breger 2000;Handler 2009;Uytterhoeven et al. 2011;Holdsworth et al. 2014). Highamplitude δ Scuti stars (hereafter HADS) are a subclass of δ Scuti stars that usually have larger amplitudes (∆V ≥ 0 m · 3) and slower rotations (Breger 2000). 1 Most of the HADS show single or double radial pulsation modes, and some of them have three radial pulsation modes or even some nonradial pulsation modes (Niu et al. 2013(Niu et al. , 2017Xue et al. 2018;Wils et al. 2008;Xue & Niu 2020). Theoretically, the period variation rates ((1/P )( dP/ dt)) of HADS caused by the evolution should roughly have a value from 10 −10 yr −1 on the MS to 10 −7 yr −1 on the post-MS phase (Breger & Pamyatnykh 1998). According to the seismic models, the observed period variation rates of some HADS (such as AE UMa and VX Hya (Niu et al. 2017;Xue et al. 2018Xue et al. , 2022) can be interpreted by the evolutionary effect self-consistently, and each of these stars should be located in the Hertzsprung gap on Hertzsprung-Russell (H-R) diagram with a helium core and a hydrogen-burning shell. In these works, the linear period variation rates are obtained from the ground-based time-series photometric data accumulated over several decades, but only those of the fundamental pulsation modes have a sufficient precision to be confirmed and tested because of the quality of the data. The high-precision photometric data from the Kepler space telescope lasted for about 4 years and could provide an excellent opportunity to determine the linear period and amplitude variation rates not only of the fundamental pulsation mode, but also of other pulsation modes. These observed quantities would tell us more secrets about the stellar evolution and pulsation in this special rapidly evolving phase.
KIC 6382916 (J19480292+4146558, also known as GSC 03144-595) is an HADS that has been monitored extensively by the Kepler space telescope for its pulsation properties, which shows three independent frequencies (pulsation modes) in its light curves (Ulusoy et al. 2013). These three pulsation modes were at first identified as the nonradial l = 1 pulsation modes in Ulusoy et al. (2013), and then identified as the fundamental, first overtone, and second overtone radial p-modes (Mow et al. 2016). Using the long-cadence (LC) photometric data of KIC 6382916 lasted from BJD 2454953 to 2456424 (Quarter 0-17), we extract the first 23 frequencies with the largest amplitudes (which have a cutoff of 1.3 mmag, see Appendix A and Table A1 for more details) to study their variations and interactions. These frequencies are also identified as the labels of the corresponding pulsation modes in this work, from F0 to F22. Theses 23 frequencies are composed of three independent frequencies, F0, F1, and F5 (also labeled as f 0 , f 1 , and f 2 ) together with their harmonics/combinations.

METHODS AND ANALYSIS
In order to extract the variation information of the amplitudes and frequencies over time, we use the short-time Fourier transformation to deal with the LC photometric data. The prewhitening process in a time window of 120 days is performed when the window is moving from the start to the end time of the LC data, with a step of 20 days. In each step, the frequencies are fixed as that have been extracted in the complete LC data (see Table A1), while the amplitudes and phases are obtained by the nonlinear least-square fitting. Then, we get the amplitude and phase variations in the 23 frequencies. Using the Fourier-phase diagram method (Paparo et al. 1998;Bowman et al. 2016;Xue et al. 2018), the variations in the frequencies can be derived from the variations in the phases.
When the dominating quasi-periodic signals in the variations in amplitudes and phases have been fully handled, we obtain their nonperiodic linear and quadratic variations of them (see Appendix B for more details). 2 The amplitude of f 0 (A 0 ) is relatively stable, that of f 1 (A 1 ) slightly decreases while that of f 2 (A 5 ) distinctly increases by about 28% in 4 years. On the other hand, the phases (subtracted by its average) of f 0 , f 1 , and f 2 vary with the same trend but different levels, which corresponds to the increasing periods of them with different linear period variation rates: (1/P 0 )( dP 0 / dt) = (3.0±1.2)× 10 −7 yr −1 , (1/P 1 )( dP 1 / dt) = (3.2 ± 1.0) × 10 −7 yr −1 , and (1/P 2 )( dP 2 / dt) = (2.4 ± 0.6) × 10 −6 yr −1 . The harmonics/combinations of f 0 , f 1 , and f 2 also show complex amplitude and frequency variations, which are summarized in the appendix, Table A1.
Noting the complex amplitude variations in the F0-F22 pulsation modes, it is quite interesting to study the amplitude interactions between them, which also indicate the energy transfer between them. Here, we introduce the interaction diagram (see Figure 1) to show the amplitude interactions between the pulsation modes. In Figure 1, the color of a small square represents the correlation coefficient 3 between the amplitudes of the labeled pulsation modes whose columns and rows intersect at the square. The dendrograms on the upper and left represent the agglomerative hierarchical clustering (AHC) process (Murtagh & Contreras 2012) in the amplitude interaction space, which is spanned by the vectors (corresponding to each of the pulsation modes) consisting of the correlation coefficients between a specific pulsation mode to all of them. In general, the AHC process reorders the pulsation modes in order to cluster those who have similar behaviors in the interaction space together. More details can be found in Appendix C.
At first glance of Figure 1, the pulsation modes are clustered into two groups (separated between F22 and F7), which correspond to the main features of increase in and periodicity of the amplitudes. In the upper group, all the f 2 related modes are closely clustered together, which have obvious increasing amplitudes. On the other hand, it is obvious that F1, F6, F7, and F10 are f 1 related and anticorrelated with the f 2 related modes, which indicates that they have decreasing amplitudes. What is more interesting is that F12, F21, F19, F8, and F22 are not f 2 related, but clustered together with the f 2 related modes. After a closer look on these pulsation modes in the frequency domain, we find an interesting discovery: almost all these modes have partners ∆ω = 0.0815 c d −1 away from them. 4 Inspired by this phenomenon, we carefully review all the 23 pulsation modes in the frequency domain carefully, and find that all the f 2 related modes have partners ∆ω = 0.0815 c d −1 away from them. 5 It indicates that f 2 could be a nonradial modes with l = 1 6 and all these f 2 related pulsation modes and their partners are locked by the star's rotation (with a rotation frequency of ∆ω = 0.0815 c d −1 ), which are interacting with each other. Additionally, all the f 2 related modes (F5, F11, F13, F14, F15, and F18) have partners ∆ω = 0.0815 c d −1 smaller than them, and the −f 2 related modes (F16) have partners ∆ω = 0.0815 c d −1 larger than them, which implies that f 2 has an azimuthal order of m = 1. If this is true, the linear period variation rates of all f 2 (including −f 2 ) related modes with significant nonzero values should be dominated by the slowed-down rotation of the star, which gives rise to the linear period variation rates that are about an order of magnitude greater than those of the purely f 0 and f 1 related modes without partners. Based on the above inferences, it is easy to understand the significant nonzero (1/P )( dP/ dt) values reaching up to an order of magnitude from 10 −6 to 10 −5 , which are seriously influenced by the slowed-down rotation and the interactions with their partners. In Figure A2, A3, A4, and A5, we show the pulsation modes without significant partners; if their (1/P )( dP/ dt) have significant nonzero values (F0, F1, F2, F4, and F7), these values fall into a range of about 2.4 × 10 −7 yr −1 to 4.0 × 10 −7 yr −1 , which should be the values of global linear period variations that are unaffected (or minimally affected) by the rotation and interactions of the modes. Considering that f 0 and f 1 have almost the same period variation and have a period ratio (P 1 /P 0 = f 0 /f 1 ) of 0.763 (Petersen & Christensen-Dalsgaard 1996;Poretti et al. 2005), they could be considered as the fundamental and first overtone radial pulsation modes, respectively.
Single-star nonrotating evolutionary models are constructed by using different initial masses with three groups of [Fe/H] (0.040 (Mathur et al. 2017), 0.111 (Xiang et al. 2019), and 0.266 (Luo et al. 2019)) from observations (which are labeled by Z = 0.014, Z = 0.016, and Z = 0.022; see more details in Appendix D), from the pre-main sequence to the start of the red giant branch. At each of the steps on the evolutionary tracks, the pulsation frequencies are calculated based on the corresponding stellar structure. Figure 2 shows the best-fit seismic models of the observed independent frequencies (f 0 and f 1 ) together with the corresponding evolutionary tracks for different Z values. The detailed information of the best-fit seismic models and the observations are collected in Table 1.
The results of the seismic models show that KIC 6382916 is located in the post-MS phase with the mass from 1.81-1.87 M . The most incredible thing is that the observed linear period variation rates of f 0 and f 1 are 3-4 times larger than the theoretical predicted ones. If we ascribe these observed linear period variation rates to the stellar evolution, we find that the star is evolving more rapidly than theoretical prediction and the standard stellar evolution theory cannot precisely describe the star's evolution process in this rapidly evolving phase. Another possible scenario is that the large linear period variations might be caused by mass transfer, which could be induced by the stellar wind or an undiscovered binary companion of the star. These possibilities should be studied and confirmed in the future. Table 1. Information of the best-fit seismic models and the observations. Z = 0.014 Z = 0.016 Z = 0.022 Observed f 0 ( c d −1 ) (l = 0, n p = 0) 4.909573 4.909774 4.910073 4.909842 ± 0.000001 f 1 ( c d −1 ) (l = 0, n p = 1) 6.432067 6.431822 6.431710 6.431890 ± 0.000001 (1/P 0 )( dP 0 / dt) ( yr −1 ) 7.78 × 10 −8 8.83 × 10 −8 1.197 × 10 −7 (3.0 ± 1.2) × 10 −7 (1/P 1 )( dP 1 / dt) ( yr −1 ) 7.49 × 10 −8 8.50 × 10 −8 1.150 × 10 −7 (3.2 ± 1.0) × 10 −7 f + 2 ( c d −1 ) (l = 0, n p = 2) 8.  Based on the above best-fit seismic models from the constraint of f 0 and f 1 , we calculate the frequencies of the second overtone mode, which are systematically greater than the observed f 2 and   labeled as f + 2 (see in Table 1). The best-fit results of the nonradial (l = 1, m = 0) modes to the observed f 2 − ∆ω have n p = 2, n g = 23 (for Z = 0.014), n p = 2, n g = 21 (for Z = 0.016), and n p = 2, n g = 19 (for Z = 0.022), which are mixed modes of the p-mode and g-mode types (see in Table 1) and labeled as f − 2 . Although in the Z = 0.016 case the observed f 2 can be roughly fitted by f − 2 + ∆ω and its large linear period variation rate can be ascribed to the slowed-down rotation of the star, the distinctly increasing amplitude of the f 2 mode remains an unsolved problem in this framework.
Here, we propose a framework to interpret all the important characteristics of f 2 . Firstly, a resonance happens between a radial p-mode and a nonradial mixed mode, i.e. f + 2 and f − 2 , and generates a resonating integration mode (RI mode) f ± 2 . The requirements of the resonance can be met by (i) f + 2 and f − 2 having close frequency values and (ii) their radial displacements having close nodes and similar phases along with the radius in the outer layers (see in Figure 3). Secondly, the RI mode f ± 2 should have the properties of radial and nonradial pulsation modes, which splits into two frequencies (f ±1 2 and f ±2 2 , f ±2 2 − f ±1 2 = ∆ω) because of the star's rotation. f ±1 2 resonates and superimposes with the combination of f 0 /f 1 (i.e. −f 0 + 2f 1 ) to form the observed partner of f 2 (i.e. F8), while f ±2 2 corresponds to the observed f 2 .
In the above framework, the frequency of f ± 2 should be modulated between f + 2 and f − 2 (similar to the situation in Goupil et al. (1998)), and then f 2 and its partner (f ±2 2 and f ±1 2 , respectively) should represent the modulations, which are covered by current long-time scale data sets and could be shown out in a short-time scale data set. Because the linear period variation of f 2 is dominantly affected by the slowed-down rotation of the star, it shows a rate about 8 times larger than the global ones (that of f 0 and f 1 ). Taking it a step further, we note that the theoretical linear period variations in f + 2 and f − 2 have opposite directions in the stellar evolution. The distance between f + 2 and f − 2 decreases with the evolution of the star, and the resonance between them will be strengthened and the amplitude of f 2 will increase. 7 About the interactions between the f 2 related modes and their rotationally locked partners, on one hand, all of them have components inherited from the RI mode, which gives rise to increasing amplitudes. On the other hand, each of the observed partners is the resonance result of the f 0 /f 1 related mode and the radial component of the RI mode, which will also cause energy transfer between them and then amplitude modulation. All these factors together produce the interaction relationships in Figure 1. A more detailed quantitative theory about the framework should be constructed in future works.
For KIC 6382916, following the current amplitude variation rates, the amplitude of f 2 will exceed that of f 0 at about BJD 2460098 (2023 June), and exceed that of f 1 at about BJD 2460196 (2023 September). This prediction can be tested in the near future by the following photometric observations, which also could also provide us with an opportunity to witness the stellar evolution process of a single star gradually in this special evolutionary phase.

DISCUSSIONS AND CONCLUSIONS
The period variation rates of HADS have been studied for a long time as their large amplitudes are always dominated by one pulsation mode, and the classical O-C method can be put into practice based on the accumulation of times of maximum light of about several decades (Breger & Pamyatnykh 1998). In recent years, based on the continuous high-precision time-series photometric data from space telescopes like Kepler, the linear period variation rates of HADS can be extracted directly from the light curves, not only for the dominating pulsation modes, but also for other nondominating pulsation modes and the harmonics/combinations of these modes (see, e.g., Breger & Montgomery (2014); Bowman et al. (2016Bowman et al. ( , 2021) . These linear period variation rates ((1/P )( dP/ dt)) generally have the values from ±(10 −9 − 10 −6 ) yr −1 , while some of them are not that reliable and would be updated by additional high-quality photometric data. Theoretically, the period variation of an HADS can be caused by stellar evolution (Niu et al. 2017;Xue et al. 2018Xue et al. , 2022, mass transfer (Xue & Niu 2020), light traveling time effect in a multiple system (Xue & Niu 2020), and pulsational nonlinearity 7 Here, we suppose a mechanism that increases the global linear period variation rates of f 0 and f 1 and also works for f + 2 and f − 2 to some extent, but keeps the signs of them unchanged.  . Radial displacement (represented by rρ 1/2 δr) against the fractional radius (r/R) for the best-fit seismic models with different Z values. r is the radial distance to the center; ρ is the density at r; δr is the radial displacement eigenfunction at r; R is the radius of the star. (Breger & Montgomery 2014;Zong et al. 2016;Bowman et al. 2016). In these possible origins, the first three will affect all pulsation modes globally 8 , while the last one will affect on the pulsation modes involved in the nonlinear interactions individually. In general, the observed period variation of an HADS is the result of a combination of several or all factors above. Especially, as a potential marker of stellar evolution, it is quite interesting and important whether the observed linear period variation of an HADS can be ascribed to the stellar evolution (Niu et al. 2017;Xue et al. 2018Xue et al. , 2022. Harmonics/combinations of independent pulsation modes are common amongst pulsating stars, such as δ Scuti stars (Breger & Montgomery 2014), β Cep stars (Degroote et al. 2009), SPB stars (Pápics et al. 2017), γ Dor stars (Kurtz et al. 2015), and pulsating white dwarfs (Wu 2001). Mathematically, the harmonics/combinations come from the nonsinusoidal properties of light curves, which indicate the nonlinearity of the star's pulsation physically. 9 In general, the harmonics/combinations are not considered as intrinsic stellar pulsation modes (Brickhill 1992;Wu 2001) and should mimic the behaviors of the parent's pulsation modes. As a result, in practice, the harmonics/combinations are always removed in the prewhitening process and not taken as the key information for the asteroseismology. Although the period variation rates of HADS for different pulsation modes have been studied by some pioneering works (see, e.g., Breger & Montgomery (2014); Bowman et al. (2016)), those of harmonics/combinations and the global interactions between harmonics/combinations have received little attention, which could provide us with abundant information about the nonlinear characteristics of both the harmonics/combinations themselves and the interactions between them.
For KIC 6382916, each of the period variations in the 23 pulsation modes are dominated by a quasi-periodic plus a linear signal (for phase variation, quadratic signals). The quasi-periodic signals are not the focus of in this work (which might have instrumental origin, or due to the light-time effect in a binary system), and will be studied in future works. About the linear period variation, unlike in AI Vel (Walraven et al. 1992), the fundamental and first overtone modes (f 0 and f 1 ) have similar period variation rates, which are 3-4 times larger than the theoretical predictions. Considering the similar values on both the f 0 and f 1 modes (and their harmonics/combinations without significant partners), the most possible origins of this result is a faster stellar evolution than that expected in this specific phase for the mass loss of the star, which needs further studies based on observational evidences. f 2 has a (1/P )( dP/ dt) value 8 times larger than those of f 0 and f 1 , which could come from the slowed-down rotation of the star and indicates that f 2 is a nonradial mode. Meanwhile, some of the harmonics/combinations show different period variations from their parent's modes (such as F22), which could be ascribed to the nonlinear interactions to the partners.
Beside the large linear period variation rates, the f 2 mode also has a distinctly increasing amplitude (see in Figure A7). Moreover, in the interaction diagram of amplitudes (see in Figure 1), the f 2 related modes are closely clustered together, which also represents the particularity of f 2 . Inspired by the results of theoretical calculations, we propose a framework to explain the strange behaviors of f 2 qualitatively. In this framework, the observed f 2 is the nonradial component of an RI mode (split by the rotation of the star), while the RI mode comes from a resonance between a radial p-mode and a nonradial mixed mode. In this case, the linear period variation rate of f 2 is dominated by the slowed-down rotation and the increasing amplitude comes from the gradually enhanced resonance along with the stellar evolution. In this sense, we can say that we witness the stellar evolution process of a single star gradually in this special evolutionary phase.
In the above framework, the variation rate of the rotation ( d∆ω/ dt) can be obtained either directly by the observed frequency variations in the f 2 related modes and their partners (the differences between them), or indirectly by the observed frequency variations in the three independent frequencies (f 0 , f 1 , and f 2 ) based on the frequency decomposition relation (see Appendix E for more details). Moreover, if we assume the conservation of angular momentum in this evolutionary phase and a uniform rotation of the star, the variation in the structure along with the stellar evolution will also cause a slowed-down rotation. The variation rates of the rotation frequency ( d∆ω/ dt) obtained from these estimations are shown in Figure 4, more details can be found in Appendix E. In this work, five pairs can be used to calculate d∆ω/ dt: F3-F14, F5-F8, F6-F11, F9-F16, and F15-F22. In Figure 4, all the five pairs indicate a slowed-down rotation, although pairs F3-F14 and F6-F11 show large uncertainties and have no definite answer. Pairs F5-F8 and F9-F16 have consistent results with that obtained from f 0 , f 1 , and f 2 (F0-F1-F5). The different results between the pair F15-F22 and F0-F1-F5 might indicate the nonlinear interaction between F15 and F22 or the influence from a hidden g-mode on pair F15-F22 in this low-frequency domain. On the other hand, it is obvious that the slowed-down rotation frequency (| d∆ω/ dt|) rates from theoretical models are significantly smaller than those from observations, which can be ascribed to a more effective angular momentum transfer based on the differential rotation of the star or an extra angular momentum loss induced by the stellar wind or its undiscovered binary companion. It needs more precise theoretical models and observational data for further study.
Figure 4. Variation rates of the rotation frequency from different estimates. The red dots with error bars represent those obtained from the five pairs. The blue region represents those obtained from the variation of f 0 , f 1 , and f 2 . The three dashed lines represent those from the best-fit seismic models of different Z values.
It is interesting to take a closer look on F22. On the one hand, the frequency variation rate of F22 can be obtained directly from Table B3 (ḟ 22 = (19.12 ± 10.05) × 10 −8 c d −2 ). 10 On the other hand, we have identified f 22 = −2f 0 + 2f 1 , and then we should haveḟ 22 = −2ḟ 0 + 2ḟ 1 . Based onḟ 0 = (−4.0 ± 1.5)×10 −9 c d −2 andḟ 1 = (−5.7±1.8)×10 −9 c d −2 , we can get −2ḟ 0 +2ḟ 1 = (−3.4±4.6)×10 −9 c d −2 . Even taking the large uncertainties into account, it shows thatḟ 22 = −2ḟ 0 + 2ḟ 1 , which indicates that the observed F22 mode is not a simple superposition of f 0 and f 1 . The nonlinear interactions from its partner or a hidden g-mode might be the origin of its significantly large frequency variation. This is not the only case. For δ Scuti stars KIC 9700322 and KIC 5950759, one can also find partners similar to the case in this work (equally spaced frequency multiplets caused by the rotation of the star) beside the dominant radial pulsation modes (Breger et al. 2011;Yang et al. 2018). We believe that this phenomenon is common in slowly rotating HADS, which evolve into a special phase and meet the resonant conditions, and the above framework of the RI mode will be refined in future research.
About KIC 6382916, there are still many interesting questions to be addressed, such as the origin of the quasi-periodic signals in the amplitudes and phases; the origin of the large linear period variation rates of f 0 and f 1 compared with the theoretical prediction; more quantitative details of the large amplitude increases in f 2 (28% in 4 years); whether the amplitude of f 2 will exceed those of f 0 and f 1 in the future; the origin of the observed large slowed-down rotation rate compared with the theoretical predictions; the inconsistency of the frequency variation rates of the harmonics/combinations and their parent's pulsation modes (such as the case in F22). More detailed research and observations are needed.

ACKNOWLEDGMENTS
We would like to thank the anonymous reviewer for his/her professional and detailed suggestions and Jue-Ran Niu for providing us with an efficient working environment. J.S.N. acknowledges support from the National Natural Science Foundation of China (NSFC) (No. 12005124 and No. 12147215 Software: MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018  The long-cadence (LC) photometric data of KIC 6382916 from the Kepler space telescope were used in this work, which covers from BJD 2454953 to 2456424 (Quarter 0-17) (publicly available PDC data Smith et al. 2012)). We downloaded the light curves (in the format of reduced BJD and magnitudes) of KIC 6382916 from the Mikulski Archive for Space Telescope (MAST) 11 , which were then normalized to be zero in the mean for each quarter. An overview of all the normalized LC data in the time domain and frequency domain are shown in Figure A1. All the above normalized LC data were prewhitened to extract the frequencies, amplitudes, and phases of the pulsation modes. This process was cutoff until an amplitude smaller than 1.3 mmg, and we obtained a total of 23 pulsation modes (see in Table A1), which is much higher than the typical noise level in Kepler data of this star (∼ 0.03 mmag). This choice ensured that the uncertainties of the amplitudes and phases were typically smaller than the variations in some pulsation modes Table A1. Frequency Solution of the First 23 Frequencies with the Largest Amplitudes. δ is defined as the difference between the observed frequency of F0-F22 mode and the marked one decomposed by f 0 , f 1 , and f 2 . dA/ dt denotes the linear variation rate of the amplitudes.

ID
Frequency of interest in 4 years, which was determined by the subsequent short-time Fourier transformation results. In each step in the prewhitening process, the Lomb-Scargle algorithm (VanderPlas 2018) was used to help to find the initial value of frequency with largest amplitude, and then a non-linear least square fitting was performed to get the final values of the frequency, amplitude, and phase. The frequencies within 1.0 ≤ f ≤ 24.4 c d −1 were considered in this work, which was determined by the LC Nyquist frequency (Bowman et al. 2016).
In order to show the detailed structures in the 23 pulsation modes in frequency domain, a zoom-in of them is shown in Figure A2, A3, A4, and A5.
The short-time Fourier transformation (Bowman et al. 2016;Zong et al. 2016Zong et al. , 2018 was then performed to the normalized LC data to get the variation in the amplitudes and phases. In this process, a time window of 120 days was moving from the start to the end time of the LC data, with a step of 20 days. In each step, the prewhitening process was performed to extract the amplitudes and phases of the specific 23 pulsation modes, while the frequencies were fixed as the values obtained in the complete LC data in Table A1. At last, the amplitude and phase (subtracted by its average value) for each of the 23 pulsation modes in each of the moving window were collected (with the times that were defined as the midpoints of  Figure A2. Zoom-in of the 23 pulsation modes in the frequency domain, Part I. The vertical dash lines mark the position of the 23 frequencies and that are ∆ω = 0.0815 away from them on both sides. Gray indicates there is no obvious peaks in that position; red indicates the pulsation mode in that position shows a positive value of (1/P )( dP/ dt); blue indicates the pulsation mode in that position shows a negative value of (1/P )( dP/ dt); green indicates the pulsation mode in that position shows an uncertain sign of (1/P )( dP/ dt). The blue circle on the top indicates the amplitude of the pulsation mode is uncertain; the upward-pointing red triangle indicates the amplitude of the pulsation mode is increasing; the downward-pointing blue triangle indicates the amplitude of the pulsation mode is decreasing. the window). The variations in the amplitudes and phases of the 23 pulsation modes are presented in Figure A6, A7, A8, A9, and A10.
In this work, the prewhitening process was performed by the Fourier decomposition, which can be presented by the formula    Figure A4. Zoom-in of the 23 pulsation modes in the frequency domain, Part III.   Figure A6. Variation in the amplitudes and phases of the 23 pulsation modes, Part I. The best-fit results to Equations B4 and B3 are presented by the black solid lines. The 2σ (deep red and blue) and 3σ (light red and blue) bounds are also shown in the panels. In the lower panel of each subfigure, the σ eff is defined as (Q obs − Q cal )/σ, where Q obs and Q cal are the values that come from the observation and model calculation, respectively; σ is the uncertainty of the observed points.    B. VARIATIONS IN THE AMPLITUDES AND PHASES Considering the quasi-periodic signals in most of the amplitudes and phases, which are not the major goals in this work, we fit the variation in the amplitudes by the following formula, The variation in the phases are fit by a similar formula, In such case, the main quasi-periodic signals (from 200 to 800 days) can be fully handled, and the linear and quadratic variations in the amplitudes and phases can be obtained.
The Markov Chain Monte Carlo (MCMC) algorithm (Sharma 2017) was used to determine the coefficients and their uncertainties in the above expressions, and the fitted coefficients are listed in Table B2 and B3.
The variations in amplitudes and phases (subtracted by the average values) are presented in Figure  A6, A7, A8, A9, and A10. In Figure A6, A7, A8, A9, and A10, we find that almost all these amplitudes and frequencies are modulated by quasi-periodic signals, which have the periods from about 300 to 760 days. It might come from the orbital perturbation and other annual systematic modulation from the Kepler space telescope, or the light-time perturbation from a companion in a binary system. Because we focus on the linear variation of the pulsation modes in this work, the detailed discussion about the quasi-periodic signals will be proposed in our upcoming works.
Based on the Fourier-phase diagram method (Xue et al. 2018), the linear period variation rates can be directly derived via the values of c φ (which are equal to the half value of the linear frequency variation rates in c d −2 ).
Because there exist overlaps in the time domain between the data points in Figure A6, A7, A8, A9 and A10 (a time window of 120 days and a step of 20 days) and these data points are not independent, the uncertainties of the parameters in Table B2 and B3 are underestimated. Consequently, we simulated the above process using the data points without overlaps (selecting the next point 120 days away from the last one), and find that the uncertainties of the quadratic term coefficient (c A and c φ ) are about 2-3 times larger than the original ones. For convenience, a uniform factor √ 6 ≈ 2.5 is used to rescale the uncertainties of (1/P )( dP/ dt) and dA/ dt in Table A1. (1/P )( dP/ dt) and dA/ dt are also used as the criteria to mark the amplitudes and frequencies variations in Figure A2, A3, A4, and A5.
One should note that, as the criterion of the increase or decrease in the amplitudes, the dA/ dt in Table A1 is different from the b A in Table B2. The dA/ dt (which is equivalent to the b A0 ) is obtained by the following formula: As a comparison of the results based on short-time Fourier transformation, the frequencies and amplitudes of the 23 pulsation modes in each year (2009)(2010)(2011)(2012)(2013) are summarized in Table B4. Table B2. Fitting results of the amplitudes of the 23 pulsation modes.  Table B3. Fitting results of the phases (subtracted by its average) of the 23 pulsation modes.
The interaction strength of a characteristic quantity (amplitude/frequency/phase) between two pulsation modes is measured by the correlation coefficient of their variation over time. In this work, we use the Spearman correlation coefficient (which is a measure of how well the relationship between two variables can be described by a monotonic function) rather than the Pearson correlation coefficient because of (i) the dataset is nonnormal; (ii) it is more robust for the outliers; (iii) it is more suitable for the measurement of the nonlinear relations, which focus on the relationships of the trends rather than the magnitudes of the variations.
The correlation coefficients of the Fi (0 ≤ i ≤ n − 1, i ∈ Z) mode with all the n pulsation modes (including itself) form a vector (k i = (ρ i 0 , ρ i 1 , ρ i 2 , ..., ρ i n−1 ), (0 ≤ i ≤ n − 1, i ∈ Z)) in the interaction space, which is a n-dimensional Euclid space spanned by n unit vectors e i = (0, ..., 1 (i-th) , ..., 0), (0 ≤ i ≤ n − 1, i ∈ Z). 12 Based on the above structures, one can use the Euclidean distance between the k i s to measure the similarities of each pulsation modes in the interaction space. Moreover, the agglomerative hierarchical clustering (AHC) process (Murtagh & Contreras 2012) can be performed to cluster the pulsation modes who have similar behaviors in the interaction space together (i.e., group a set of modes in such a way that modes in the same group (called a cluster) are more similar (in some sense) to each other than to those in other groups (clusters)).
The color of a small square represents the correlation coefficient of the characteristic quantity between the two labeled pulsation modes whose column and row intersect at the square. The dendrograms on the upper and left represent the AHC process in the interaction space, which reorders the pulsation modes in order to cluster those who have similar behaviors together.
All the relationships and structures shown in the interaction diagram of amplitudes indicate that there exist interactions or energy transfer between the independent pulsation modes and their harmonics/combinations. Moreover, the interaction diagram of amplitudes help us find the particularity of f 2 related modes, discover their rotationally locked modes, and identify f 2 as a nonradial l = 1, m = 1 mode, which indicates its potential in future research. On the other hand, comparing with the complicated interactions shown in Figure 1, the Interactions Diagram of the phases (see Figure C11) shows that almost all the pulsation modes have the same variation trend except for the F16, F22, and F17 modes, which are caused by the decreasing trend of their periods. Excluding the period variation rate of F17 with large uncertainty, that of F16 can be interpreted by its −f 2 component and that of F22 could be interpreted by its interaction with F15. D. THEORETICAL MODEL CALCULATION In order to determine the stellar mass and evolutionary stage based on the single-star evolutionary models, the open-source 1D stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, r15140, Paxton et al. (2011, 2018) was used to construct the structural and evolutionary models. At each step along with the evolutionary tracks, the pulsation frequencies of the specific structure were calculated by the stellar oscillation code GYRE (Townsend & Teitler 2013;Townsend et al. 2018;Goldstein & Townsend 2020).
The initial parameters used to construct pre-main sequence evolutionary models of KIC 6382916 were configured as follows. Different metallicity [Fe/H] with values of 0.040 (Mathur et al. 2017), 0.111 (Xiang et al. 2019), and 0.266 (Luo et al. 2019) dex were considered as the initial metallicity of the evolutionary models. The following formulas were used to calculate the initial heavy element abundance Z and initial hydrogen abundance X: where X = 0.7381 and Z = 0.0134 (Asplund et al. 2009). Eq.(D6) was provided by Mowlavi et al. (1998). Based on the given values of [Fe/H], we got (X = 0.704, Z = 0.014), (X = 0.695, Z = 0.016), and (X = 0.672, Z = 0.022), as the initial inputs of the evolutionary models. At the same time, the values of the option Zbase (which provides the reference metallicity necessary to calculate element variations) were set to be the initial metallicity of the star (0.014, 0.016, and 0.022, respectively), and initial zfracs = 'AGSS09 zfracs' was selected. For the opacity, kap file prefix = 'a09', kap CO prefix = 'a09 co' and kap lowT prefix = 'lowT fa05 a09p' were selected. The initial mass of the models was set in the interval from 1.5 M to 2.5 M with a step of 0.01 M , covering the typical mass range of HADS. The rotation of the star was not considered in the model calculation because of the relative slow rotation of the star (∆ω = 0.0815). The mixing-length parameter was set to be a value of α MLT = 1.89 (Niu et al. 2017). The exponential scheme proposed in Herwig (2000) was adopted to account for the overshooting, and the overshooting parameter was set to the fixed value f ov = 0.015 (Magic et al. 2010).
All the evolutionary tracks were calculated from the pre-MS to the red giant branch. Based on the pulsation frequencies calculated in every step along with the evolutionary tracks, we got the best-fit seismic models (which have the smallest χ 2 with respect to the observed values of f 0 and f 1 ) with the different Z values (see in Table 1 for detailed information). The H-R diagram, log T eff vs log g diagram, and Petersen diagram are shown in Figure 2, Figure D12 and Figure D13 Figure D12. log T eff vs. log g diagram with all the evolutionary tracks and the best-fit seismic models of KIC 6382916 with different Z values. The colored solid evolutionary tracks show from the zero-age MS to the post-MS evolutionary stage, with initial masses from 1.50 M to 2.50 M . The regions surrounded by the black rectangular boxes are selected to zoom in the best-fit seismic models, which are represented by the red stars. The colored dash lines present the evolutionary tracks with initial mass steps of 0.01 M .  Figure D13. Petersen diagram with all the evolutionary tracks and the best-fit seismic models of KIC 6382916 with different Z values. The colored solid evolutionary tracks show the stages from the zero-age MS to the post-MS evolutionary stage, with initial masses from 1.50 M to 2.50 M . The regions surrounded by the black rectangular boxes are selected to zoom in the best-fit seismic models, which are represented by the red stars. The colored dashed lines present the evolutionary tracks with initial mass steps of 0.01 M .
E. VARIATION OF ∆ω Based on the results in Table B3, five pairs (F3-F14, F5-F8, F6-F11, F9-F16, and F15-F22) can be used to calculate the variation rate of the rotation frequency ( d∆ω/ dt), which can be obtained by the differences in the frequency variation rates (2 · c φ ) in each of the pairs. At the same time, all the above pairs follow the frequency relation which can be differentiated on both sides with respect to time: d∆ω dt ≡∆ω =ḟ 0 − 2ḟ 1 +ḟ 2 .
(E9) f 0 ,ḟ 1 , andḟ 2 can also be found in Table B3. About the uncertainty estimation, all the related frequency variation rates were taken as independent variables, and a uniform factor √ 6 ≈ 2.5 was used to rescale the uncertainties as that in Appendix B. The numerical results were collected in Table E5.
For the best-fit seismic models of different values of Z, the variation rates of the rotation frequency were calculated based on an assumption of the conservation of angular momentum in this evolutionary phase and a uniform rotation of the star, whose numerical results were also collected in Table E5.  Model (Z=0.022) −1.4 × 10 −9