K2 photometry on oscillation mode variability: the new pulsating hot B subdwarf star EPIC 220422705

EPIC 220422705 had been observed by K2 in short-cadence (SC) mode over a period of 78.72 days during Campaign 8. Its assembled light curves were downloaded from Mikulski Archive for Space Telescopes¹¹1https://archive.stsci.edu/k2 (MAST). These archived data were processed through the EVEREST pipeline²²2The EPIC Variability Extraction and Removal for Exoplanet Science Targets as developed by Dr. R. Luger which is an open-source pipeline for removing $\sim 6.5$ hr instrumental systematics in K2 light curves. One can see details through the link: https://archive.stsci.edu/hlsp/everest.. The photometry corrected by EVEREST has comparable precision to the original Kepler mission for targets brighter than $K_{p}\approx 13$ (Luger et al., 2018). EPIC 220422705 is within this brightness range.

The EVEREST flux was firstly shifted to the relative fraction to their mean value. We then used a six-order polynomial fitting to detrend the whole light curve due to residual instrumental drifts. To avoid discontinuities in the light curve across gaps longer than 0.02 days, we separated the light curve piece-wise where such gaps occurred for our fitting. This detrending method will flatten the light curves and dismiss signals with period ($\gtrsim 2$ d) in Fourier transform, which will not have impact on the modulating patterns for our prime aim. We note that those signals are not concerned here due to the fact that the EVEREST pipeline may not recover those signals correctly. Then the light curves were iteratively clipped of a few outliers three times by filtering at 4.5$\sigma$ around the light curve before we produced a Fourier transformation. Figure 1 (top panel) shows the final light curve of EPIC 220422705 which contains 106,444 data points over a duration of 76.43 days with a 1.2 d gap in the middle. The amplitude scatter clearly reveals multiperiodic signals of hours in a close-up view (middle panel). The corresponding Lomb-Scargle periodogram (LSP; Lomb, 1976; Scargle, 1982) up to the Nyquist frequency is shown in the bottom panel where the g-mode frequencies are clearly dominant in a region of [$\sim 100-1000$] $\mu$Hz.

We used the specialized software FELIX³³3Frequency Extraction for Lightcurve exploitation, developed by S. Charpinet, greatly optimizes the algorithm and accelerates the speed of calculation when performing frequency extraction from dedicated consecutive light curves. See details in Charpinet et al. (2010, 2019) and Zong et al. (2016b, a). to perform frequency extraction from the light curves. The frequencies were prewhitened in order of decreasing amplitude until the value of 5.2 times the local noise level, a value that is the median amplitude in the LSP (Zong et al., 2021). This detection threshold is adopted as a compromise between 2-yr Kepler and 27-d TESS photometry (Zong et al., 2016b; Charpinet et al., 2019). The highest peak will be extracted in the case where there are several close frequencies of $<0.4\,\mu$Hz, i.e., about $3\times\Delta f$ ($\Delta f=1/T$, and T $\sim 76.43$ days). We have detected 66 independent frequencies and 13 linear combination frequenciess, with the highest (1172 ppm) frequency at 279.767 $\mu$Hz, which are listed in Table 1.

Pulsating stars with both acoustic p- and gravity g-mode oscillations are excellent candidates to probe their internal profiles since acoustic and gravity waves propagate in different regions of the stellar structure (Aerts et al., 2010; Kurtz, 2022). From spaceborne photometry, some g-mode dominated sdBVs are found with low-amplitude p-mode pulsations (see, e.g., Baran et al., 2017; Zong et al., 2018; Sahoo et al., 2020). A direct and easy way to distinguish the two different types of mode is by their pulsation period. In general, theoretical sdB star calculations suggest that dipole p-modes typically have frequency $>2500$ $\mu$Hz ( $P<400$ s), whereas g-modes $<1000$ $\mu$Hz ($P>1000$ s) (Fontaine et al., 2003; Charpinet et al., 2005, 2011). But p-mode frequencies can decrease below 1700 $\mu$Hz (periods can increase beyond 600 s) as $T_{\mathrm{eff}}$ and $\log g$ decreases (Charpinet, 1999; Charpinet et al., 2001, 2002).

Figure 2 shows preliminary classification for the p- and g-mode regions based only on the period. We detect 53 independent frequencies in the range [$\sim 80-1000$] $\mu$Hz which are clearly g-mode ($P>1000$ s) pulsations (two top panels). Another 7 independent frequencies are found in the high frequency p-mode ($P<400$ s) region, [$\sim 2600-3800$] $\mu$Hz (bottom panel). There are six independent frequenciesin the region of [$\sim 1100-1400$] $\mu$Hz or [$\sim 715-910$] s, which might be low-order high-degree ( $\ell>3$) g-modes or mixed modes that need further classification (see, e.g., Charpinet et al., 2011, 2019). Those frequencies, hardly directly classified to be p- or g-mode by merely of their frequency value, can be used to penetrate a much larger portion of stellar interior or to detect the differential rotation in radial or longitude. However, determining the exact modes requires an exploration of seismic models. We note that a few frequencies were detected in this intermediate region in sdB stars observed with Kepler photometry, for instance, KIC 3527751 and KIC 10001893 (Foster et al., 2015; Uzundag et al., 2017). In addition, we have resolved 13 linear combinations with frequencies $<1400~{}\mu$Hz which could be intrinsic resonant modes (Zong et al., 2016a) or non-linear effects from the linear eigenfrequencies (Brassard et al., 1995).

From linear perturbation theory, an eigenmode of oscillation can be characterized by spherical harmonics that are described by three quantum numbers: the radial order $n$, the degree $\ell$, and the azimuthal order $m$. When a star rotates, the degenerated $m$ components will split into $2\ell+1$ multiplets. Referring to Ledoux (1951), their frequencies are related by,

where $\nu_{n,l,0}$ is the frequency of the central $m=0$ component, $\Omega$ is the solid rotational frequency, and $C_{n,\ell}$ is the Ledoux constant. For acoustic p-mode, $C_{n,\ell}$ is very near to zero and can be ignored, whereas it is estimated as $C_{n,\ell}\sim 1/\ell(\ell+1)$ for high-radial order gravity g-modes.

To resolve any rotational split multiplet from spaceborne photometry, a minimum criterion is that the observations should cover at least twice the rotation periods. Charpinet et al. (2018) and Silvotti et al. (2021) present the distribution of rotation periods for sdB stars determined from Kepler photometry. Frequency multiplets found rotation periods from a few days to near one year with most having periods a bit longer than one month. This indicates that K2 photometry can likely resolve frequency multiplets in sdB stars.

Figure 3 shows the frequency spacings of 12 groups of frequencies. We first consider six g-mode frequency groups that are detected with close frequency spacings around 0.17 $\mu$ Hz, which we consider to be dipole modes, i.e., $f_{5}\sim$ 123.09 $\mu$Hz, $f_{14}\sim$ 180.39 $\mu$Hz, $f_{41}\sim$ 198.89 $\mu$Hz, $f_{23}\sim$ 288.56 $\mu$Hz, $f_{19}\sim$ 324.16 $\mu$Hz and $f_{30}\sim$ 501.25 $\mu$Hz. The weighted (by $1/\sigma f_{i}$) average value is $0.168\pm 0.016~{}\mu$Hz. Those dipole modes give a rotational frequency of 0.33 $\mu$Hz which would mean quintuplet splitting of $0.28~{}\mu$Hz using $C_{n,1}=1/2$ and $C_{n,2}=1/6$. Three g-modes, $f_{18}\sim$ 89.48 $\mu$Hz, $f_{54}\sim$ 622.79 $\mu$Hz and $f_{24}\sim$ 697.63 $\mu$Hz, have close frequency spacings of $\sim 0.23$ $\mu$Hz, which could be rotational quintuplets, considering frequency uncertainties. We also resolve three $\ell=1$ p-mode multiplets, $f_{52}\sim$ 3706.12 $\mu$Hz, $f_{39}\sim$3720.77 $\mu$Hz and $f_{62}\sim$3741.96 $\mu$Hz, with low-amplitude peaks at a frequency distance of $\sim 0.4~{}\mu$Hz.

In order to determine the rotational period in a quantitative way, we propose a new approach that defines the rotation period associated with errors by their probability of occurrence. As shown in Figure 3, we adopt a Gaussian distributions, $N\sim(\delta f_{i},\sigma f_{i}^{2})$, to represent the probability of a group of resolved frequencies at their shifted values. Here $\delta f_{i}$ is the relative frequency to the central component. The probability of 68.27% (i.e., 1$\sigma$ in $N\sim(0,1)$) was calculated to define the values and the uncertainties of frequency spacings. We obtained the values of $0.170\pm 0.049$ $\mu$Hz, $0.234\pm 0.035$ $\mu$Hz and $0.399\pm 0.132$ $\mu$Hz for $\ell=1$, and $2$ g-modes and p-modes, respectively. The corresponding rotation periods are $34.04_{-7.12}^{+13.78}$ d, $41.21_{-5.35}^{+7.22}$ d and $28.86_{-7.21}^{+14.34}$ d, respectively.

Our result suggests that a slightly differential rotation occurs in EPIC 220422705 as g- and p-modes probe stellar interior under different depth (see, e.g., Kurtz, 2022). We note that our results are completely based on only marginally-resolved frequencies with low-amplitudes near the detection limit. We do not detect multiplets in the four highest-amplitude frequencies or 12 of the 13 highest-amplitude frequencies. Nevertheless, the rotation period we determine is consistent with that of a typical sdB star (see details in Charpinet et al., 2018; Silvotti et al., 2021). Foster et al. (2015) claims to detect a differential rotation in KIC 3527751 whose core rotates slower than the envelope, which, however, was challenged by an independent analysis of the same photometry by Zong et al. (2018), citing that the claimed rotational p-mode multiplets had missing components under significant confidence. In reverse, Kawaler & Hostler (2005) suggest that the core might rotate faster than the envelope of sdB stars from evolutionary models. In the case of EPIC 220422705, a binary system but with a poorly measured orbital period (Barlow et al., 2012), a slightly faster rotating envelope could be interpreted that the orbital companion has accelerated it via the tidal force. Theory predicts that angular momentum transportation leads to radiative envelope first synchronized then gradually proceeds to the inner part (Goldreich & Nicholson, 1989). In combination with the poorly-determined orbital and rotation period (because of low-amplitude multiplets), it would be unwise to speculate too much on this star. There are other sdBV stars which would be better for such work. To be cautious, EPIC 220422705 can still be a rigid object if the uncertainties of rotational periods are fully considered.

For g-mode pulsations in the asymptotic regime, consecutive high-radial orders ($n\gg\ell$) follow a pattern of equal period spacing (see, e.g., Aerts et al., 2010), which depends on the structure. Seismology theory provides the following relationship,

with $\Pi_{0}$ defined as,

where $N$ is the Brunt-Väisälä frequency and $r$ is the radial coordinate. For the period spacing of $\ell$ = 2 sequence, it is related to the $\ell$ = 1 sequence as,

Previous analysis of sdBV stars from Kepler and TESS reveals that the period spacing is about 250 s and 150 s for dipole and quadruple modes, respectively (see, e.g., Reed et al., 2011; Sahoo et al., 2020). To find the spacing periods in EPIC 220422705, we performed the popular Kolmogorov-Smirnov (KS) test on the independent g-mode frequencies. The KS test returns spacing correlations as highly-negative values for the most common spacings in a dataset (Kawaler, 1988). We first apply the KS test to the six rotational (incomplete) triplets, which has a deepest trough at 63 s = 252/4 (Figure 4 a). In addition, we perform a linear fitting to those 6 periods with a result of $\sim$265.5 s. Then we applied another KS test for the 53 independent frequencies lower than 1000 $\mu$Hz as probable g-modes, which gives a value around 276.8 s for dipole modes (Figure 4 a). All values are consistent with that of dipole modes in sdB stars (Reed et al., 2011; Sahoo et al., 2020).

Based on the preliminary period spacings, we have identified nine modes as dipole, 13 as quadruple, and additional 8 frequencies which fitted both period sequences. We note that the four identified modes include one of the above 6 dipole modes, $f_{5}\sim 123.1~{}\mu$Hz. We obtained the average period spacing of 268.5 $\pm$ 2.8 s and 159.4 $\pm$ 0.6 s for $\ell=1$ and $\ell=2$ modes, respectively, via linear fitting to 17 (4 + 5 + 8) dipole modes and 21 (13 + 8) quadruple modes (Figure 4 b). We list $\ell$ and relative $n$ values in Table 1. We note that the real radial order can only be obtained through seismic modeling. Our results suggest that EPIC 220422705 has a somewhat large period spacing among the known sdB variables, in a range of [220, 270] s for the dipole mode (see, e.g., Reed et al., 2011; Sahoo et al., 2020; Uzundag et al., 2021). As stellar models presented in Uzundag et al. (2021), the evolution tracks suggest that the lower value of period spacing, the lower value of $\log g$. This agrees well to a low $\log g$ derived for EPIC 220422705. However, atmospheric parameters, with a much higher precision, are encouraged for EPIC 220422705 to test the results of Uzundag et al. (2021) in future.

The échelle diagrams for two sequences are presented in Figure 4 (d) where the $\ell=2$ sequence is more consistent than the $\ell=1$ sequence. There are two $\ell=1$ frequencies with larger period deviations, $f_{19}\sim 324.16\,\mu$Hz, $f_{30}\sim 501.25\,\mu$Hz whereas only one occurs in the $\ell=2$ sequence, $f_{8}\sim 218.36\,\mu$Hz. Asymptotic theory indicates that period spacings are determined by the size of the pulsation resonant cavity (see e.g. Tassoul, 1980) and could be affected by the extent of the convective core (Smeyers & Moya, 2007). The ideal pattern for a period sequence in the échelle diagram is a vertical ridge for central components of a star with internally-homogeneous composition. Some deviations from the mean period spacing are to be expected in g-mode pulsating sdB stars – and indeed have already been unambiguously detected in some cases (Østensen et al., 2014) – due to the phenomenon of mode trapping by steep composition gradients. Such deviations reflect properties of the stellar interior, and in particular could reveal the character of the mixing processes at work near the core boundary (see, e.g., Charpinet et al., 2011; Ghasemi et al., 2017). Figure 4 (c) shows the deviation of period spacing as a function of the reduced pulsation period. We only observe a large deviation that might be associated to a trapped mode at the fifth order. Indeed, seismic models suggest that strong trapping are more likely found for the lower order (higher frequency or shorter period) modes than the higher order g-modes (Charpinet et al., 2014).

A quick look of all AM/FM patterns occurring in these 22 frequencies suggest that most of them exhibit simple or quasi-regular variations (Figures 5, 6, 7 and 8). In order to quantitatively characterize the modulation patterns, we apply three simple types of fittings: linear, parabolic, and sinusoidal waves. We first calculate standard deviations of our sub-set data, $\sigma_{0}$. Then we adopt a simple AM/RM fitting, typically linear first and then a second type on the residuals. This fitting method may be iterated several times until the residuals present no clear structure and look like random noise. For each fitting, a standard deviation of residuals will be calculated as:

where $G_{k}(t)$ defines as,

Here $k=1,2,3,4,5$ and $N$ denotes number of data points, $m_{i}$ are the measured values of amplitude or frequency.

In principle, the freedom degree of $G_{k}(t)$ increases as $k$ increases, which in turn results in a lower $\sigma_{k}$. However, to avoid overfitting of the modulation patterns, for instance, by including a linear or parabolic fitting, we follow the statistical test by Pringle (1975) who defines $\lambda$ as,

Here $D_{i}$ is the freedom degree of the fitting function, e.g., 2 and 3 for linear and parabola fitting, respectively. A significant improvement for the higher freedom degree fitting should have the parameter that meets the $F-$distribution, $\lambda>F_{P=99.75\%}(D_{2}-D_{1},N-D_{2})$. For each AM/FM, we prefer to keep the fitting function with a lower freedom degree as indicated by the parameter $\lambda$. The results of our fitting for each AM/FM is listed in Table 2 where most frequencies have a sinusoidal component, $G_{3}(t)$, or with an additional linear fitting, $G_{4}(t)$.

After the fitting function was set, we adopted the posterior distributions based on the Bayesian frame to estimate the best fitting parameters and their uncertainties. The posterior distributions of parameters are sampled by the Markov Chain Monte Carlo (MCMC), which is implemented by the EMCEE code (Foreman-Mackey et al., 2013). The MCMC sampling is performed with more than 22 chains, according to the number of free parameters. It proceeds until the chains converge to values that are inferred by the auto-correlation time module of the EMCEE. Parameters of the best fittings are given by the medians of marginalized posterior distributions, associated with the corresponding errors that are calculated by the half-widths between the 16^th and 84^th percentiles of the distributions.

Figure 5 is an example of the MCMC application for AM and FM occurring in the frequency $f_{3}\sim 328.75~{}\mu$Hz. We clearly see that both AM and FM exhibit periodic variations but with different periods and phases. In general, we find that almost all measurements are consistent with the fitting curves, accounting for uncertainties. The best fittings return $A=8.61\pm 0.54$ ppm and $6.1\pm 0.3$ nHz, $T=14.4\pm 0.1$ d and $20.4\pm 0.2$ d, $\phi=5.2\pm 0.1$ and $1.4\pm 0.1$, $b=-0.91\pm 0.05$ and $-0.08\pm 0.02$, and $c=514\pm 1$ and $2.4\pm 0.4$ for AM and FM, respectively.

In this section, we describe the modulation patterns in AM/FM of all 22 significant frequencies except the frequency of $f_{3}$ which has already been described. Figure 6 is a gallery of 16 frequencies with various modulation patterns, which are identified by different fitting function $G_{k}(t)$ and summarized in Table 2. In general, the fitting periods of AMs and FMs are on timescale of months.

We now specifically describe these AMs and FMs. A few modes have been observed with simple modulation patterns of only linearly decreasing or increasing in amplitude or frequency. For instance, the amplitude of $f_{9}\sim 143.4~{}\mu$Hz decreases from 200 ppm to 100 ppm over the time interval of about 50 days. In other observed minor cases, a few modes can be characterized by a simple parabolic fitting, such as the AM/FM of $f_{6}\sim 306.2~{}\mu$Hz where both the amplitude and frequency reach to their vertex values near the time BJD = 2457420 d. The majority of modulations we observed have quasi-sinusoidal patterns, most with extra linear and a few with parabolic fittings. For example, the amplitude and frequency of $f_{8}\sim 218.4~{}\mu$Hz exhibit completely sinusoidal pattern and evolve in phase. Whereas the frequency of $f_{13}\sim 299.0~{}\mu$Hz demonstrates a sinusoidal pattern with an additional linear fitting. We only find $f_{16}\sim 139.9~{}\mu$Hz whose frequency and $f_{10}\sim 187.5~{}\mu$Hz whose amplitude show a sinusoidal plus a parabolic pattern. We note that the amplitude of $f_{14}\sim 180.434~{}\mu$Hz decreases below the detection threshold from $\mathrm{BJD}-24554000=200$ d to 225 d, which holds even if our subsets span 70 d. For most frequencies, a variation of frequency scale is around 20 nHz but there are a few exceptions with large values around 100 nHz. For instance, $f_{20}\sim 2781.672~{}\mu$Hz, $f_{19}\sim 324.164~{}\mu$Hz and $f_{23}\sim 288.547~{}\mu$Hz. The variations of amplitude can span up to a few hundred ppm ($f_{1}\sim 279.767~{}\mu$Hz) or down to a few ppm ($f_{4}\sim 261.269~{}\mu$Hz). We note that the amplitudes of AM and FM are not strictly proportional to each other. For instance, the amplitude of $f_{1}$ increases from about 50 ppm to 250 ppm, but the frequency only varies around 1 nHz.

A very interesting feature of the observed AM/FM is that several frequencies are found with (anti-) correlations between AM and FM. We thus calculate the values of correlation for all frequencies as listed in Table 2. The frequencies $f_{23}\sim 288.561\mu\mathrm{Hz}$ and $f_{6}\sim 306.223~{}\mu$Hz all exhibit strong correlation with a coefficient $|\rho_{\mathrm{a,f}}|>0.8$ between their amplitude and frequency variation. For instance, $f_{23}$ and $f_{6}$ show clear anti-correlation between AM and FM, with derived coefficients of -0.867 and -0.824, respectively. For $f_{23}$, the amplitude exhibits a slightly decreasing trend whereas the frequency is opposite. For $f_{6}$, its AM and FM are both represented by parabolic fits but with anti-phase evolutions. We also note several frequencies with correlated AM and FM patterns showing sinusoidal variations, such as $f_{15}\sim 173.4~{}\mu$Hz and $f_{8}\sim 218.3~{}\mu$Hz.

Figure 7 shows another 4 frequencies, $f_{5}\sim 123.087~{}\mu\mathrm{Hz}$, $f_{7}\sim 80.595~{}\mu\mathrm{Hz}$, $f_{11}\sim 312.0~{}\mu\mathrm{Hz}$ and $f_{39}\sim 3720.7~{}\mu\mathrm{Hz}$, whose FM patterns still present modulation structures after fitting by a $G_{k}(t)$ function with $k=4,5$. To remove those FM residuals, they needed additional sinusoidal function, i.e., three fitting functions to present FM patterns. For a strong correlation frequency, $f_{5}$, with $\rho_{\mathrm{a,f}}=0.862$, exhibits regular AM and FM with period of $101^{+14}_{-14}$ and $40^{+1}_{-1}$ days, respectively, derived from EMCEE. The other three frequencies, $f_{7}$, $f_{11}$ and $f_{39}$ are determined by EMCEE with either different fitting types of functions or different periods of sinusoidal patterns between AM and FM. For instance, the periods of AM and FM in $f_{7}$ are calculated as $23^{+1}_{-1}$ and $49^{+1}_{-1}$ days, respectively, almost with a ratio of $1:2$. The EMCEE returns periods with values of $\sim 10.3,19.0,17.4$ and $11.3$ days for the additional sinusoidal fittings of the FM residuals of $f_{7}$, $f_{5}$, $f_{11}$ and $f_{39}$ after extracting $G_{k}(t)$, respectively. We note that FM periods of these residuals are shorter than the periods fitted by $G_{k}(t)$.

Figure 8 presents the AM and FM of the frequency $f_{24}\sim 697.62~{}\mu$Hz that shows a quasi-regular behavior but cannot be described by simple fittings of $G_{k}$ functions, which is also suggested by the fine profile of the LSP and the sliding LSP (sLSP). Thus we do not perform fitting and EMCEE on its AM and FM. We observe that both the AM and FM began with a decreasing trend: the amplitude went down from $\sim 80$ ppm to a local minimum value of $\sim 60$ ppm and the frequency varied from +120 nHz down to -100 nHz relative to its averaged frequency. Then the frequency and amplitude had experienced an increasing trend. The amplitude reaches to its maximum with a time interval of about 30 days but passing across one stationary point, whereas the frequency generally went up to around the average value with a back and forth trend.

As presented in Section 2.3, we detected several potential rotation multiplets with frequency spacings of about 0.2 to 0.4 $\mu$Hz, which is comparable to the frequency resolution determined by the observation duration, $\Delta f\approx 1/T$. This finding suggests that close frequencies, not only the multiplets we found, may be present in EPIC 220422705 with very high probability. In literature, many sdBV stars are resolved with such nearby frequencies from the original 4-yr Kepler observations (see .e.g, Zong et al., 2018; Foster et al., 2015; Baran et al., 2012). Those close frequencies will induce amplitude variations not frequency variations, called beating, if they are unresolved, as demonstrated in Zong et al. (2018). Here we only take two close frequencies as an example,

If we consider two comparable amplitudes, $A_{1}\sim A_{2}=A$, and a very small frequency separation, $\Delta\omega\ll\omega_{0}$, the equation (8) will be reduced to,

The above equation indicates that two close frequencies will generate amplitude modulations if they were not resolved. We thus have to evaluate this effect on the observed AMs in EPIC 220422705.

A series of simulations had been performed to quantitatively compare the observed and calculated AMs. The simulation process are similar to that of Zong et al. (2018), including light curve construction and frequency injection, but adopted to K2 observations. In practice, we set three parameters $A_{1}$, $A_{2}$ and $\Delta\omega$ to be constant in construction of each light curve based on equation (8) and phase to be zero for simplicity. We then make all parameters as variants in a series of 1331 light curves, with $A_{1},A_{2}\in[0.5,0.6,...,1.5]\times A_{0}$ and $\Delta\omega\in[0.1,0.12,0.14,...0.3]~{}\mu$Hz. Then we transform the simulated light curves into sLSP and directly compare them with the observed sLSP via direct subtraction. Both types of sLSPs have to be transformed into the same time step and window length and scaled to the same maximum amplitude. The flatness of the residual sLSP defines the goodness of similarity between the two sLSPs. We finally select an optimal sLSP for the simulation and return the parameters.