© The Authors 2026

1. Introduction

Stellar pulsations have been studied for more than a century since 1919, when Eddington began to explore the theory of stellar pulsation with particular interest in understanding Cepheid variables (Eddington 1919a,b, 1926). For the first half-century of the study of stellar pulsation, it was assumed that the surface geometry of pulsation modes could be described adequately by spherical harmonics (Cowling 1941), which are the eigenmode solutions to pulsation equations for perfectly spherical stars. The geometry of eigenmodes associated with the spherical harmonics is described by three quantum numbers: n, the radial overtone, which is the number of radial nodes that are spherical shells; ℓ, the degree, which is the number of surface nodes of the mode; and m, the azimuthal order, which gives the number of surface nodes that are lines of longitude. It was further assumed that the axis of pulsation coincided with the rotation axis of all pulsating stars, since that, in most cases, is the principal axis of distortion from spherical symmetry (see, e.g., Kurtz 2022 and Aerts et al. 2010 for reviews).

Traditionally, in stars other than the Sun, the degree of the mode is inferred, rather than observed directly. Surface cancellation of observed sectors of the spherical harmonics has meant that modes of higher degree (ℓ ≥ 3) have been difficult to detect, since these higher-degree modes have very low visibility (Dziembowski 1977). On the other hand, because the Sun has a well-resolved surface, modes of very high degree can be detected to beyond ℓ = 300 (Larson & Schou 2011, Scherrer et al. 2012, Larson & Schou 2015). Observation and knowledge of these high-degree modes place major constraints on the structure of the Sun via helioseismology.

Since we know that the Sun pulsates in very high degree modes, we can reasonably conjecture that other stars do this as well. The question then is how these high-degree modes affect the structure of a pulsating star, its evolution, its internal rotation, and angular momentum transfer. The answers to these important questions can be illuminated by the direct observation of modes of degree ℓ ≥ 3.

Kurtz (1982) discovered the rapidly oscillating Ap (roAp) stars and demonstrated for the first time that some stars have pulsation axes that are inclined to the rotation axis of the star. In the case of the strongly magnetic roAp stars, this pulsation axis is close to the magnetic axis, which is itself inclined to the rotation axis. Kurtz (1982) proposed the oblique pulsator model to explain the pulsation amplitude and phase variations seen to occur in roAp stars with rotation. Briefly, this shows that as an obliquely pulsating star rotates, the geometry of the degree of the pulsation mode becomes visible because the mode is seen from differing aspect. In effect, observers can in thought walk around the star and see the nonradial modes from different aspects, thus exposing the number and positions of the surface nodes.

More recently, data from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) enabled the identification of two δ Scuti stars in binaries with pulsations largely confined to one hemisphere with respect to the tidal axes on the L1 side (Handler et al. 2020; Kurtz et al. 2020). These were dubbed “single-sided pulsators”. HD 74423 had one prominent pulsation mode with a half-dozen components separated by the orbital frequency ν_orb and a systematic change in pulse phase by π radians every orbit. CO Cam had four such modes, each with components also separated by ν_orb. In these systems, there were systematic phase changes when the pulsation amplitude was near zero in each of these modes.

These modes were discovered largely because of their visually striking single-sidedness, but this complication also made it difficult to understand the nature of the underlying modes themselves from the Fourier transforms (FTs). At the time, it was assumed that the axes of these pulsations lie along the tidal axis, that is, that they were tidally tilted modes, wherein the pulsations were tidally trapped on one hemisphere and suppressed on the other (Fuller et al. 2020).

Rappaport et al. (2021) reported a much simpler multiplet mode in a δ Scuti pulsator in a binary system, without a single-sidedness. This mode had two prominent equal-amplitude components separated by 2 ν_orb with two weak peaks on either side of the larger peaks and separated by 1 ν_orb. There were two maxima in this pulsation mode around the orbit (at the maxima of the ellipsoidal light variation; ELV), and phase jumps of π rad at the eclipses. This was interpreted as a tidally tilted Y_11x mode, where the x refers to the pulsation axis lying along the tidal axis (the $\widehat{x}$ direction)¹.

Building on these discoveries, Jayaraman et al. (2022) identified 31 pulsation modes in a pulsating sdB star with a white dwarf companion. These modes typically had between three and five components separated by ν_orb, with one or two π phase jumps per orbit. These were not simple to understand, but were interpreted at the time as tidally tilted Y_1mx or Y_2mx modes.

More recently, TESS data revealed two δ Scuti stars in binary systems with rich pulsation spectra, including 9 and 14 doublets (Zhang et al. 2024; Jayaraman et al. 2024), which we interpreted as simple dipole pulsation modes. Each of these modes consists of just two dominant components separated by 2 ν_orb with either no detectable central component, or with a weak component. Each mode has two amplitude maxima per orbit at either the eclipses or the ELV maxima. Each also has two π phase jumps per orbit at either the ELV maxima or the eclipses. These were interpreted as Y_10x or Y_10y modes, depending on where the modes had their amplitude maxima.

Fuller et al. (2025) modeled these systems as an entirely new type of mode that consisted of linear combinations of Y_1 + 1z and Y_1 − 1z modes that are coupled by tidal, centrifugal, and Coriolis forces in the binary. This modeling revealed that linear combinations of Y_1 + 1z and Y_1 − 1z modes produce, and are mathematically identical to, tidally tilted modes Y_10x and Y_10y. Thus, either interpretation is valid. However, as Fuller et al. (2025) also showed, other higher-order modes with pulsation axes along z (i.e., with ℓ = 2 and ℓ = 3) can also be combined via tidal, centrifugal, and Coriolis forces to produce new eigenmodes of the star. These modes also have multiplet frequencies that are spaced by integer multiples of ν_orb, and can exhibit up to four and six amplitude maxima per orbit, as well as up to four and six π phase jumps, respectively.

These latter Fuller modes are all linear combinations of spherical harmonics with pulsation axes lying along the spin (or orbital) axis. To the lowest order, modes of the same ℓ and differing by Δm = 2 are coupled by the tidal forces (Fuller et al. 2025). We note that for ℓ ≥ 2, the degeneracy is broken, and these are no longer equal to or even mimic simple pulsation modes that have been tidally tilted along either the x- or y-axes. Thus, it is now possible within this paradigm to distinguish between tidally coupled z modes and tidally tilted modes whose pulsation axes lie along $\widehat{x}$ or $\widehat{y}$.

All of the tidally tilted and triaxial modes mentioned above can be reinterpreted as Fuller modes, mostly with ℓ = 1 and ℓ = 2 components that may or may not be substantially suppressed in one of the tidal hemispheres. We report the discovery of a δ Scuti pulsator in a binary system (TIC 287869463) that exhibits a Fuller octupole pulsation mode that cannot be described as any type of tidally tilted pulsation. While ℓ = 3 pulsation modes have been found in red giants (Stello et al. 2016), as well as in some main-sequence solar-like oscillators (e.g., Kjeldsen et al. 2005) and were suggested in some β Cephei stars (e.g., Aerts et al. 1998, Cotton et al. 2022), no such modes have been robustly detected in δ Scuti stars. Moreover, no previously discovered ℓ = 3 sectoral mode is stationary, that is, non-circulating.

In Sect. 2 we present observational evidence for these modes and show how the amplitudes of these modes vary with time, while the phase differences between the upper and lower frequency components remain constant to within the statistical uncertainties. Simulations of these and other Fuller modes are presented in Sect. 3, wherein we highlight the striking similarity between the simulated frequency peaks and what we observe in TESS data. An analysis of the binary system parameters is undertaken in Sect. 4. We summarize and discuss our results in Sect. 5. In the appendix, we discuss O-C diagrams for the pulsations and show that the pulsation frequencies increase steadily over the 5 yr interval of the TESS observations.

2. Discovery of Fuller modes in TIC 287869463

We are currently carrying out an extensive search for Fuller-mode pulsators in TESS data. In particular, we are examining a set of bright likely eclipsing binaries (TESS magnitude T_mag < 13.5) selected by machine learning from 56 sectors of TESS observations with a 10-minute cadence or shorter (Kostov et al. 2025; Powell 2026, private communication; Jayaraman et al., in prep.). In all, there are 51 820 such candidate binaries whose T_eff falls in the range for them to contain possible δ Scuti pulsators (6500 K ≤T_eff ≤ 9000 K). The TESS light curves used for the search are taken from the MIT Quick-Look Pipeline (QLP; Kunimoto et al. 2021, 2022).

For each candidate binary, we produce an automated echelle diagram that plots the pulsation frequency against the echelle phase (the pulsation frequency ν_puls modulo ν_orb, normalized to ν_orb). To do this, we first determined the orbital period P_orb, if any, and then subtracted the first 30 harmonics of ν_orb from the light curve. We then sequentially identified and removed the 75 highest-amplitude pulsation frequencies in the FT from the remaining light curve, or the highest-amplitude pulsations down to a noise floor of five times the rms amplitude in the FT, whichever came first. The echelle diagrams were created from these highest-amplitude peaks. Finally, we examined each echelle diagram by eye for interesting modes that were split by integer multiples of the orbital frequency. Most of the split modes are dipoles (split by 2ν_orb), along with a smaller percentage of quadrupoles (split by 4 ν_orb). However, we identified one target (TIC 287869463, or Gaia DR3 5216608316411701248; Stassun et al. 2019) as having an octupole mode, split by 6 ν_orb.

A short 5 d segment of the TESS light curve for TIC 287869463 is shown in Figure 1 (top panel). The stellar pulsations of the hotter component can be seen superposed on the eclipsing light curve. In all, there are seven sectors of data with a 200 s cadence for a total of about 168 d of data, with 100 d of that being contiguous. A Fourier-reconstructed light curve made from the first 60 orbital harmonics is given in Figure 1 (bottom panel). Table 1 lists all TESS sectors in which TIC 287869463 was observed, alongside the observing cadence. Table 2 enumerates photometric and astrometric information about TIC 287869463 and the system properties (modeled in Section 4).

Fig. 1. Light curves for TIC 287869463. Top: 5-d segment of the raw TESS light curve. The pulsations superposed on the eclipsing light curve are readily apparent. Bottom: Fourier-reconstructed light curve from the first 60 orbital harmonics. We used only the cosine terms in the reconstruction to remove a small time-varying O’Connell effect (O’Connell 1951), presumably arising from star spots on the cooler companion star.

The FT of the 200 s cadence data from sectors 63–66 (100 d of nearly continuous data) is shown in Fig. 2. The data were cleaned of the first 40 orbital harmonics before computing the FT. The FT shows six prominent peaks due to stellar pulsations, all in the frequency range ∼35–40 d⁻¹. The dipole pulsation modes (separated by 2ν_orb) are labeled D1 and D2, and the octupole mode, separated by 6ν_orb, is labeled O1.

Fig. 2. Fourier transform of the TESS data from sectors 63–66 for TIC 28786963. The labels mark the two dipole modes (D1 and D2) and an octupole mode (O1). Each mode has two prominent components, marked a and b. The dipole components are separated by 2 νorb, and the octupole components are separated by 6 νorb.

The echelle diagram for the pulsations of this star is shown in Figure 3. The D1 mode has a small central component, but D2 does not. The octupole mode is split by 6 ν_orb and has only one tiny intermediate frequency peak, which is situated at 1 ν_orb from the missing central frequency. There are also two low-amplitude singlet modes at 35.976 and 41.096 d⁻¹.

Fig. 3. Echelle diagram for TIC 287869463 constructed from the TESS data of sectors 63–66. This shows the frequencies of the pulsations vs. echelle phase, i.e., the frequency modulo νorb and normalized to νorb. The labeling of the components follows Fig. 2.

The six most significant frequencies seen in the periodogram and echelle diagram, related to D1, D2, and O1, are enumerated in Table 3. These frequencies were produced from phase tracking the pulsations (see, e.g., Rowman et al. 2016, and also Appendix A) over all 12 sectors of 200 s and 600 s cadence data (where the pulsation frequencies were below the Nyquist frequency). All frequencies are referenced to a common epoch, BJD = 2 460 000. In the fifth column of Table 3, we list the mean frequency split between the two components of the particular mode. In the last column, we indicate the (orbital) period inferred from the splitting, and these values all agree to within 10⁻⁵ d of the calculated orbital period².

The amplitudes and phases of the six pulsations as a function of TESS sector are given in the top and bottom panels of Fig. 4, respectively. The amplitudes of the two components of each of the dipoles are clearly correlated, as are the two amplitudes of the octupole mode. However, the latter are less variable than for the D1 and D2 modes. The two amplitudes for each pulsation mode are always very similar and do not vary independently.

Fig. 4. Top: Variations in the dipole (red and blue) and octupole (green) component amplitudes with sector number across ∼5 yr. For each mode, the two component peaks have very similar amplitudes, and their variations are clearly tightly correlated. Bottom: Difference in phase between the two components of each of the three modes at the times of primary eclipse. To fit all three curves on the same scale, we subtracted 180 degrees from the phase difference of modes D2 and O1. Note that the phase differences remain approximately constant at either 0 or π rad.

In the phase plot (Fig. 4, bottom panel), we show the phase difference between the two components of each mode as a function of TESS sector, all referenced to a time of primary eclipse (BJD = 2 460 014.4757). For the D1 mode, the two components are mostly in phase at the primary eclipse. The two components of the D2 and O1 modes are essentially always out of phase by 180° at the primary eclipses. (This in turn means that these modes are at an amplitude minimum at primary eclipse.) The standard errors in the phase differences are only 1.1°, 0.9°, and 3.7° for the D1, D2, and O1 modes, respectively.

We next force-fit a three-element triplet in frequency (separated by ν_orb) to each of the dipole modes and a seven-element septuplet to the octupole mode. We used the amplitudes and phases of the extracted mode elements to reconstruct the amplitudes of the three pulsation modes as a function of orbital phase. These results are shown in the top panel of Fig. 5. An amplitude-scaled orbital light curve is shown superposed for phasing reference. Each of the dipoles has two amplitude maxima per orbit, while the octupole has six such maxima. We find that D1 has a maximum amplitude at the eclipses, while D2 has a minimum there. The octupole mode also has a pulsation minimum at the eclipses. These all agree with the discussion above of the phase differences of the two elements of each mode (see also Fig. 4) at the time of primary eclipse. The bottom panel of Fig. 5 shows the run of reconstructed pulse phase versus orbital phase. The dipole modes have two π phase jumps each orbit, while the octupole mode has six π phase jumps.

Fig. 5. Reconstructions of pulsation amplitude and phase vs. the orbital phase. The color-coding, top to bottom, is red = D1, blue = D2, and green = O1. A scaled version of the orbital light curve is superposed to guide the eye on the orbital phasing. We note the six maxima in pulsation amplitude and six π rad phase jumps per orbit for the octupole (O1).

Several arguments indicate that the octupole mode in this star (the first such stationary ℓ = 3 sectoral mode discovered) is a single-pulsation mode split by 6 ν_orb, as opposed to two independent modes that are coincidentally separated by ∼6 ν_orb. First, the separation of the octupole mode elements corresponds to an integer multiple of the measured orbital period to one part in 10⁵. The odds of this occurring at random for the two significant non-dipole frequencies over the pulsation range of 30–40 d⁻¹ is only 1.7 × 10⁻⁶ (verified by Monte Carlo simulations). Second, the very similar amplitudes of the two mode components and their correlated amplitude and frequency variations over a 5-year interval (see Appendix A) also strongly suggest that these are not from independent modes. Third, a minimum or maximum at the time of the primary eclipse (minimum in this case) in this mode would occur at random only about 25% of the time³. We also note that a well-developed model (Fuller et al. 2025) predicts such modes, and this adds confidence to the interpretation that this is an octupole mode.

3. Simulations of the modes and interpretation

Our analysis of the TESS data for this star showed that we have found an octupole mode of the form

$$\begin{array}{cc}\hfill Y_{33+}& =(Y_{3+3z}+Y_{3-3z})e^{i\omega t}\hfill \\ \hfill & =(y^3-3x^{2}y)sin(\omega t),\hfill \end{array}$$(1)

where the Ys on the right-hand side of the top equation are ordinary spherical harmonics with ℓ = 3 and m = ±3. The additional z designation indicates that the axis of these modes lies along the rotation vector of the pulsating star, which in turn we assume is aligned with the orbital angular momentum vector (along $\widehat{z}$). A traveling wave propagates around the stellar equator for each of these spherical harmonics. Fuller et al. (2025) showed how such modes can be coupled by perturbations on the pulsating star in a binary from tidal, Coriolis, and centrifugal forces into a stationary mode (standing wave) in the rotating stellar reference frame.

The Y₃₃₊ mode is a standing-wave sectoral octupole mode in a frame rotating with the orbital motion. Here, the plus indicates the sum, rather than the difference, between the terms on the right side of Equation (1). This mode differs from standard sectoral Y_3 + 3z or Y_3 − 3z modes by its stationary nature. Throughout the orbital cycle, we can therefore observe this mode from varying angles. This gives the same advantage of providing information on the geometry of the mode (and its viewing aspect) as is gained from oblique pulsation in roAp stars. Such geometric information cannot be extracted from rotationally perturbed normal sectoral octupole modes pulsating about the rotation axis.

Previously, we identified numerous stationary dipole modes formed from Y_11± modes (Y_1 + 1z ± Y_1 − 1z) (Zhang et al. 2024; Jayaraman et al. 2024). We are also working on a smaller set of other pulsators with related quadrupole modes formed as Y_22± = Y_2 + 2z ± Y_2 − 2z (see, e.g., Handler et al. 2025). As described in Fuller et al. (2025), there are three other perturbed stationary ℓ = 2 modes and five other ℓ = 3 modes in addition to the two described above. However, the components of these other ℓ = 3 Fuller modes are split by either two or four times the orbital frequency and can be difficult to distinguish from Y_11± or Y_22± modes.

In order to make some of the properties of the new stationary perturbed modes concrete, we developed a code to simulate these pulsations in a binary star system. We used two stationary spherical stars separated by the semimajor axis of the binary a. We uniformly distributed 7000 unit vectors on each star and then viewed the system from the perspective of an observer orbiting around the binary at a frequency ν_orb at an inclination angle i. Any pulsation mode can then be assigned to the fixed pulsating star; this mode can be either circulating or stationary⁴. The pulsation amplitudes were arbitrary since we have no way of calculating these a priori. At a given time, we evaluated the dot product between the viewing direction unit vector and the unit vectors over the surface of both stars. When the dot product was negative, then the observer sees that part of the star; all the visible flux (as calculated in this manner) was then summed. Eclipsed regions were also excluded from the sum. This same dot product was also used to compute the limb darkening and the projected area of the surface element. We also added a simple cosine term of 2ν_orb to represent an ellipsoidal light variation for aesthetic purposes.

This set of calculations was made every 120 s (to match the TESS 2-min cadence) for 50 000 steps (to roughly match three sectors of TESS data). The orbital modulations were removed from the simulated time series, as is done for the real dataset. We then took the FT of the simulated data. The entire process was then repeated for 18 orbital inclination angles between 90° and 5°. An illustrative set of simulated light curves for a Y₃₃₊ mode is shown in Fig. 6. Each light curve has six pulsational maxima per orbit with minima (in this particular case) at both eclipses. In general, the relative pulsation amplitude decreases with decreasing inclination angle and reaches zero at i = 0, as expected. We highlight the curve at i = 75° because it is close to the inclination determined for TIC 287869463 (see Sect. 4).

Fig. 6. Simulated light curves of TIC 287869463 with a Y33+ pulsation mode. The orbital inclination angles in degrees are written next to the right y-axis. The orange curve for i = 75° represents the simulated light curve for the approximate inclination angle of TIC 287869463. The curves for different inclinations are shifted vertically by arbitrary amounts for clarity.

FTs of simulated data for four different ℓ = 3 modes are shown in Fig. 7. In the upper left panel, we display the FTs for a Y₃₃₊ mode of the type detected for TIC 287869463 as a function of orbital inclination angle. There are two main peaks that are separated by 6 ν_orb. The smaller peaks for i ≳ 80° arise because the eclipses periodically modulate the intensity of the pulsations, a phenomenon known as spatial filtration (Gamarova et al. 2003; Reed et al. 2005; Biró & Nuspl 2011; Johnston et al. 2023; Van Reeth et al. 2023). For an inclination of ∼75° (as observed for TIC 287869463), these frequency peaks would not have been detected. In the upper right panel, we show the corresponding FTs for a Y₃₃₋ mode, and the results are similar to those for a Y₃₃₊ mode. The main difference is that the small intermediate peaks are somewhat more pronounced because the pulsation amplitudes are maximum at the eclipses, rather than minimum, and this enhances the effects of the eclipses on the visibility of the pulsations. Even at i ≃ 75°, however, these smaller peaks would have been only barely detectable, however.

Fig. 7. Simulated FTs for different examples of ℓ = 3 pulsation modes. The modes represented in the four panels are Y33+, Y33−, Y33x, and Y31x, clockwise starting from the upper left panel. The first two modes are defined in Sect. 3 and Equation (1). The pulsation axis of the latter two modes has been tilted into the orbital plane and lies along the tidal (or x) axis. The FT peaks are arbitrarily shifted to the left for decreasing inclination angles to avoid overlapping peaks.

A surface plot of the Y₃₃₊ mode and the Y₃₃₋ mode, as they might appear in the context of TIC 287869463, is shown in Figure 8. Note the pulsation minimum (maximum) at the time of superior conjunction in the Y₃₃₊ (Y₃₃₋) mode.

Fig. 8. Diagrams showing the flux perturbation across the surface of a star pulsating in the Y33+ mode (top) and Y33− mode (bottom). The diagram corresponds to orbital phase zero for an observer viewing at at inclination i = 55° relative to the orbital axis. The observed phase of the amplitude modulation of the O1 mode is consistent with that of a Y33+ mode.

In the lower two panels of Fig. 7, we show the FTs for simulated data of two different ℓ = 3 modes, namely a Y_33x and Y_31x mode. Here, the x signifies that the pulsation axis is tilted into the orbital plane and lies along the tidal axis. As mentioned above, these are tidally tilted pulsations, and traveling waves propagate around the x-axis. As shown by the simulations, there are numerous intermediate peaks out to ±3 ν_orb. We recall that the central peak of each multiplet is the only pulsation mode frequency; all of the other components of the multiplet are Fourier descriptions of the observed pulsation amplitude and phase modulating with changing orbital aspect. These multiplets in the lower two panels are very different from what is observed in TIC 287869463. We simulated all combinations of traveling and stationary modes around all three axes x, y, and z with ℓ = 3, and the only mode that fits the observational data (FT and the phasing of amplitudes and pulse phases with respect to the eclipses) for TIC 287869463 is a Y₃₃₊ octupole Fuller mode.

4. Analysis of the binary system

To evaluate the properties of the two stars in TIC 287869463, we simultaneously fit the spectral energy distribution (SED) curve and the TESS light curve. We employed a Markov chain Monte Carlo (MCMC) algorithm to fit for the stellar masses M₁ and M₂ (where M₁ is the more massive pulsating star), the system age, the orbital inclination angle i, interstellar extinction, A_V, and the distance. The latter two parameters were fixed at the values from Gaia. In total, there were between four and six free parameters. We assumed that the binary had evolved in a coeval fashion, that is, no prior mass transfer or mass loss from the system. We were thus able to use stellar evolution models (Choi et al. 2016; Dotter 2016) to determine the radii and T_eff of both stars, which are uniquely determined given their masses and age.

We collected the available archival spectral flux measurements (SED points) from 0.15 to 11.6 μm. These data were accessed from VizieR⁵ (Ochsenbein et al. 2000), which in turn uses the systematic sky coverage of such surveys as Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016), Sloan Digital Sky Survey (SDSS; Gunn et al. 1998), Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), Wide-field Infrared Survey Explorer (WISE; Cutri et al. 2013), Galaxy Evolution Explorer (GALEX; Bianchi et al. 2017), and Global Astrometric Interferometer for Astrophysics (Gaia; Gaia Collaboration 2021). More details about this type of SED fitting for binary star systems are given in Jayaraman et al. (2024), Handler et al. (2025), Yakut et al. (2025a,b).

The input information we used included (i) 23 available SED points spanning 0.15 μm to 11.6 μm, (ii) an approximate value for K₁ of the radial velocity (RV) curve⁶, (iii) the eclipsing light-curve profile, (iv) the Gaia distance of 1086 ± 3 pc (Bailer-Jones et al. 2021), and a Gaia value for A_G of 0.709⁷. Additionally, we used the MESA Isochrones and Stellar Tracks (MIST) evolution tracks (Paxton et al. 2011, 2015; Choi et al. 2016; Dotter 2016; Paxton 2019) to relate the stellar mass and age to its radius and T_eff. In order to compute the model stellar spectra from the radii and T_eff of the stars, we used the Castelli & Kurucz model stellar atmospheres (Castelli & Kurucz 2003), assuming solar metallicity. Finally, we corrected for interstellar extinction using the prescription presented in Cardelli et al. (1989).

When a composite SED of two or more stars is fit, it is important to have additional constraints among the properties of one star and the other(s). In particular, the eclipsing light curve in this system provides important information about parameters such as T_eff, 2/T_eff, 1, R₂/R₁, R₁/a, and the orbital inclination angle, where a is the semimajor axis of the system. Since well-developed light-curve emulators such as PHysics Of Eclipsing BinariEs (PHOEBE; Prša & Zwitter 2005) are not designed to work jointly with SED fitting codes, we incorporated our own simple light-curve emulator to help the SED fitter separate the light contribution from the two stars more properly (Jayaraman et al. 2024). This light-curve emulator uses two spherical stars with limb darkening to generate eclipses. The out-of-eclipse behavior is described by three sinusoids of the form cos(ω_orbt), cos(2ω_orbt), and sin(ω_orbt) to represent effects such as ellipsoidal light variations, illumination effects, Doppler boosting, and corotating spots on the cooler star (see, e.g., Kopal 1959; Carter et al. 2011).

The best-fitting SED and light-curve models are shown in the top and bottom panels of Fig. 9. The best-fitting system parameters derived from these fits are summarized in Table 4. The luminosity of the primary star, with a mass of 2.2 M_⊙, radius 2.6 R_⊙, and T_eff, 1 = 8734 K, is typical of a δ Scuti star, with T_eff near the blue border of the δ Sct instability strip. The pulsator parameters place it between the Balona et al. (2015)δ Scuti-category boxes 4 and 7 in their Fig. 2. The pulsations of pulsators in these categories can reach ∼50 d⁻¹, and there are many such systems with pulsations in a range similar to what we see in TIC 287869463.

Fig. 9. Top panel: Spectral energy distribution fit for TIC 287869463. The blue curve shows the model spectrum for the primary star (the pulsator), and the red curve represents the model spectrum for the cooler smaller secondary star. The black curve that nearly coincides with the red curve is the sum of the model fluxes. The orange points are measured SED points from the literature (for the SED fitting techniques and the origin of the data points, see Sect. 4). Bottom panel: Light-curve fit using our custom light-curve emulator. The black curve shows the Fourier-reconstructed light curve (i.e., the data), and the red curve shows the model fit (see text for details).

Finally, we point out that even though the secondary star contributes less than 1% of the system luminosity, we can be confident that we measured its parameters because there are numerous constraints between the secondary and the primary other than from the SED, for instance, K₁, R₂/R₁, and T₂/T₁. Therefore, the primary star largely accounts for matching the SED, while the properties of the secondary star (including its contribution to the SED) come from the relative relations listed above, and from a common age.

5. Discussion

While carrying out a large-scale search for Fuller-mode pulsators in TESS data, we discovered an unprecedented octupole mode in a δ Scuti star. The mode is characterized by two peaks in the Fourier transform separated by 6 ν_orb. This splitting corresponds to six times the orbital frequency to within an accuracy of one part in 10⁵. We argue that this cannot be just two random and unrelated modes of the pulsating star, but must comprise a single mode.

This octupole mode is actually a combination of ordinary Y_3 + 3 and Y_3 − 3 modes that are aligned with the stellar rotation axis. These are combined into a new eigenmode of the star by binary perturbations, and we call it a Fuller Y₃₃₊ mode (Fuller et al. 2025). We presented a number of compelling reasons why the octupole mode cannot plausibly be two independent modes whose frequencies are accidentally aligned in the echelle diagram. This is the first securely identified octupole mode in any δ Scuti star, and the first Fuller-type ℓ = 3 stationary sectoral mode found in any star, including the Sun.

As we show in the appendix, the pulsation frequencies of the dipole and octupole modes were found to steadily increase with time, all at different rates. However, the splits in the mode frequencies remain constant at an integer multiple of the orbital frequency.

Because of geometric cancellation of higher-order ℓ modes over the visible hemisphere of pulsating stars, the visibility of a pulsation mode decreases sharply with increasing ℓ. To investigate the visibility of higher-degree modes, we followed the approach described in Daszyńska-Daszkiewicz et al. (2002). As a first step, we approximated the mode visibility by the disk-averaging factor

$$\begin{array}{c}\hfill b_{\ell }={\displaystyle {\int }_0^1}h(\mathrm{\mu })P_{\ell }(\mathrm{\mu })\mathrm{\mu }d\mathrm{\mu },\end{array}$$(2)

where μ is the cosine of the angle between the observer line of sight and the local surface normal, and h is the limb-darkening function. This should be reasonably accurate when the effects of rapid rotation and tidal interaction are neglected. In a more general approach, as in the case of the star we studied, the mode visibility may also depend on the azimuthal order m.

We show a plot of the visibility of pulsation modes in a star like TIC 287869463 in Fig. 10 as a function of ℓ for four different filters, including a red filter, which is a reasonable approximation to the TESS bandpass (the results are also given numerically in Table 5). Dipole and quadrupole modes (ℓ = 1 and ℓ = 2, respectively) have relatively high visibilities of 0.7 and 0.3, respectively. However, for octupole modes (ℓ = 3), the visibility drops abruptly to 0.037 (for the red filter). In TIC 287869463, the dipole mode amplitudes are 1.8 mmag and 0.8 mmag. However, the amplitude of the octupole mode at 0.4 mmag is not too far below this. When we renormalized this according to the relative visibilities, the intrinsic amplitude of the octupole mode was found to be ∼10 mmag. Pulsation amplitudes of this size are fairly common among δ Scuti stars.

Fig. 10. Mode visibility diagram as a function of ℓ. These visibilities were calculated as in Daszyńska-Daszkiewicz et al. (2002) for the BVRI passbands; the red passband matches that of TESS best.

This exercise also showed that even hexadecapole modes with ℓ = 4 are expected to be as straightforward to detect as the octupole mode found here. However, higher modes with ℓ ≥ 5 will be very difficult to detect.

It is also worth noting that according to the simulations of Daszyńska-Daszkiewicz et al. (2006), the octupole mode should be relatively easy to identify from radial-velocity variations. Therefore, time-series spectroscopy might be used to validate our identification of the octupole mode for TIC 287869463.

Finally, we show model evolution tracks in the HR diagram for several sets of model parameters in Figure 11, along with a comparison to the stellar values we derived for the TIC 287869463 pulsator. The evolutionary models were computed with the Warsaw-New Jersey code (e.g., Pamyatnykh 1999) adopting the OPAL opacity tables (Iglesias & Rogers 1996) and the chemical mixture of Asplund et al. (2009). With the accurately determined stellar parameters, including the mass and radius, we expect strong constraints from seismic modeling of this star, for example, on the metallicity Z, the efficiency of outer-layer convection α_MLT, and the overshooting distance α_ov. This will be the subject of our follow-up studies of this object.

Fig. 11. Six different evolution tracks in the HR diagram compared to the measured luminosity and Teff of TIC 287869463. The parameters that are varied in the different tracks involve metallicity Z, αMLT, mass, rotation velocity, and convective overshooting, αov (see text for details).

Acknowledgments

RJ is currently supported by a Klarman Fellowship from the College of Arts & Sciences at Cornell University. This project has received partial funding from the HUN-REN Hungarian Research Network. T.B. acknowledges the financial support of the Hungarian National Research, Development and Innovation Office – NKFIH Grant OTKA K-147131. VBK is grateful for financial support from NASA grant 80NSSC22K0747. GH thanks the Polish National Center for Science (NCN) for support through grant 2021/43/B/ST9/02972. JDD thanks the Polish National Center for Science (NCN) for support through grant 2023/50/A/ST9/00144. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Science Mission Directorate. The QLP data used in this work were obtained from MAST (Huang 2020), hosted by the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work also presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The Gaia mission website is https://www.cosmos.esa.int/web/gaia. The Gaia archive website is https://gea.esac.esa.int/archive/. Some of the SED fluxes and magnitudes were obtained with the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Additionally, some of the SED fluxes and magnitudes were obtained with the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We used the Simbad service operated by the Centre des Donnés Stellaires (Strasbourg, France). This research has also made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. (2000).

References

Appendix A: Changes in the pulsation frequencies and the orbital period

The pulsation frequencies and the orbital period appear to be nearly constant across the 200-s cadence data from the four contiguous TESS sectors 63-66, or even to some extent when we include the final three sectors of 200-s cadence—separated from the former sectors by about 700 d. However, the apparently constant pulsation frequencies are seen to be clearly varying once we add five earlier sectors with 600-s cadence data (sectors between S30 and S39). Overall, after analyzing the available TESS data between sectors S30 and S93, the pulsation frequencies are all found to be significantly increasing with time. The orbital frequency can be tracked to even earlier times by using six sectors between S3 and S13 with 30-min cadence, and ν_orb is also found to be increasing with time. We thus sought to characterize the frequency changes in the pulsations as well as in the orbit.

To track the pulsation phases ϕ (see also Rowman et al. 2016), we fit a sine curve of the form Acos[ω(t−t_ecl)−ϕ] to represent a pulsation over an interval of, say, half of a full TESS sector (∼12 d). We needed to consider sufficiently long intervals because some of the pulsation peaks are within a fraction of 1 d⁻¹ of other peaks. The reference time, t_ecl is the time of a primary eclipse, e.g., BJD 2460014.476. We then fit for A and ϕ and stepped the fitting window forward in time by a fraction of a TESS sector.

The results of tracking the phase for each of the 6 pulsation frequencies are given in Figure A.1. Each panel displays the equivalent Eclipse Timing Variation (ETV) curve for both components of a given pulsation mode: D1a and D1b, D2a and D2b (top two panels) and O1a and O1b in the lower left panel. Note that for mode D1, the a and b components track each other very well, and we saw in Fig. 4 that indeed these two components are always in phase at the primary eclipses. By contrast, we can see from the ETVs of D2 and O1 that the ETV curves for the two components are separated by close to half a pulsation period, consistent with the fact that the two components of both modes remain out of phase at the primary eclipses over all the sectors with 10-min or better cadence.

Fig. A.1. ETV curves for each component of the three pulsation modes, D1, D2, and O1, as well as from the binary orbital period. All four curves show decreasing periods with time, but not at the same rates. All phases are referenced to the times of a primary eclipse, i.e., at BJD - 2,456,014.4757. Note that the vertical scale for the D1 ETVs is exactly 1/3 that for the D2 and O1 ETVs. In turn, the vertical scale for the orbital ETV curve is 4/5 that of the D1 ETVs.

The final panel in Fig. A.1 shows the ETV curve for the orbital period. This shows that P_orb is also apparently decreasing with time. If we interpret this as an orbital decay, then the decay timescale of just ∼2 × 10⁵ yr seems rather implausible. Alternatively, we might interpret the curve as a portion of a Doppler delay orbit caused by an orbiting third body. Unfortunately, the duration of available TESS data is insufficient to determine an outer orbit for a putative third body. However, for one of the shorter allowed orbits of 15 yr, the semimajor axis is 950 light-seconds (for an assumed circular orbit), and the corresponding mass function is f(M)≃0.03 M_⊙. For a binary mass of 3 M_⊙ (see Table 4), this yields a value of M₃sin³i of about 0.75 M_⊙. A third body of such a mass would not have contributed sufficient light to the SED to have been unequivocally detected, especially for outer inclination angles of ≳75°. We note that the known secondary in the binary contributes only ≲0.01 of the system luminosity, and a comparable third star would barely perturb the SED.

Returning to the systematic nonlinear behavior in the ETV curves for the pulsations, we examine what the consequences would be for different assumptions about the nature of the nonlinear ETV curve observed for the orbit. Suppose that a Fuller mode pulsation frequency is ν_p, and that it is physically changing with time as ${\dot{\nu }}_{\mathrm{p}}$. The splitting with the orbit will be ±m ν_orb (see, e.g., Shibahashi & Kurtz 2012). Now, we also allow that the orbital frequency is physically changing as ${\dot{\nu }}_{\mathrm{orb}}$. This means that the two components of the mode (upper and lower) should follow different functions of time:

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{p},\mathrm{u}}& ={\nu }_{\mathrm{p}}+m{\nu }_{\mathrm{orb}}+{\dot{\nu }}_{\mathrm{p}}t+m{\dot{\nu }}_{\mathrm{orb}}t\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{p},\mathrm{l}}& ={\nu }_{\mathrm{p}}-m{\nu }_{\mathrm{orb}}+{\dot{\nu }}_{\mathrm{p}}t-m{\dot{\nu }}_{\mathrm{orb}}t\hfill \end{array}$$(A.2)

In turn, if we were to measure the ${\dot{\nu }}_{\mathrm{p}}$ of the upper and lower components, we should find

$$\begin{array}{cc}\hfill {\dot{\nu }}_{\mathrm{p},\mathrm{u}}& ={\dot{\nu }}_{\mathrm{p}}+m{\dot{\nu }}_{\mathrm{orb}}\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill {\dot{\nu }}_{\mathrm{p},\mathrm{l}}& ={\dot{\nu }}_{\mathrm{p}}-m{\dot{\nu }}_{\mathrm{orb}}\hfill \end{array}$$(A.4)

with a difference in ${\dot{\nu }}_{\mathrm{p}}$ of $2m{\dot{\nu }}_{\mathrm{orb}}$. From Table A.1 we can see that the differences between the $\dot{\nu }$ terms of the upper and lower frequency components are not statistically significant, and the uncertainties are ∼10⁻⁷ d⁻². For the octupole mode, we might expect a change in the difference of the $\dot{\nu }$s by about 6 × 9 × 10⁻⁹ = 5.4 × 10⁻⁸ d⁻² due to the changing orbit. This is smaller than the uncertainty in the difference in $\dot{\nu }$ listed in Table A.1 for the octupole mode. Thus, we cannot quite use this effect to test for a decaying orbit.

For Doppler shifts of the binary (which include the pulsating star) the quantity $\dot{\nu }/\nu$ would be constant. So, we should find that all the pulsation components have ${\dot{\nu }}_{\mathrm{p}}/{\nu }_{\mathrm{p}}={\dot{\nu }}_{\mathrm{orb}}/{\nu }_{\mathrm{orb}}$, or ${\dot{\nu }}_{\mathrm{p}}\simeq 50{\dot{\nu }}_{\mathrm{orb}}\simeq 0.5\times {10}^{-6}$ d⁻². In all cases, this is less than the values of $\dot{\nu }$ for the pulsations listed in Table A.1, but would be readily detectable if we could separate the Doppler from other natural changes in pulse frequency.

Regarding the observed (or apparent) changes in the orbital ETV curve, we point out that many such close binaries containing an early-type star exhibit ETVs that have nonlinear variations on long timescales (decades or even centuries) that are almost certainly neither actual orbital period changes (i.e., orbital decay) nor third-body induced light travel-time effects. Some examples of very irregular and unexplained ETVs are the cases of XZ And (Yuan & Qian 2019), X Tri (Lee et al. 2023), RZ Cas (Lehmann et al. 2020), and Y Cam (Çelik & Kahraman Aliçavuş 2024). While we are confident in the nonlinear behavior of the TIC 287869463 orbital ETV curve, we cannot be certain what causes this behavior. Therefore, at least tentatively, we take all the $\dot{\nu }$s for the pulsations to be independent of the orbit, and for the orbital ${\dot{\nu }}_{\mathrm{orb}}$ to be due to neither an orbital decay nor an orbital Doppler shift.

All Tables

All Figures

Fig. 1. Light curves for TIC 287869463. Top: 5-d segment of the raw TESS light curve. The pulsations superposed on the eclipsing light curve are readily apparent. Bottom: Fourier-reconstructed light curve from the first 60 orbital harmonics. We used only the cosine terms in the reconstruction to remove a small time-varying O’Connell effect (O’Connell 1951), presumably arising from star spots on the cooler companion star.

In the text

	Fig. 2. Fourier transform of the TESS data from sectors 63–66 for TIC 28786963. The labels mark the two dipole modes (D1 and D2) and an octupole mode (O1). Each mode has two prominent components, marked a and b. The dipole components are separated by 2 νorb, and the octupole components are separated by 6 νorb.
In the text

	Fig. 3. Echelle diagram for TIC 287869463 constructed from the TESS data of sectors 63–66. This shows the frequencies of the pulsations vs. echelle phase, i.e., the frequency modulo νorb and normalized to νorb. The labeling of the components follows Fig. 2.
In the text

Fig. 4. Top: Variations in the dipole (red and blue) and octupole (green) component amplitudes with sector number across ∼5 yr. For each mode, the two component peaks have very similar amplitudes, and their variations are clearly tightly correlated. Bottom: Difference in phase between the two components of each of the three modes at the times of primary eclipse. To fit all three curves on the same scale, we subtracted 180 degrees from the phase difference of modes D2 and O1. Note that the phase differences remain approximately constant at either 0 or π rad.

In the text

	Fig. 5. Reconstructions of pulsation amplitude and phase vs. the orbital phase. The color-coding, top to bottom, is red = D1, blue = D2, and green = O1. A scaled version of the orbital light curve is superposed to guide the eye on the orbital phasing. We note the six maxima in pulsation amplitude and six π rad phase jumps per orbit for the octupole (O1).
In the text

	Fig. 6. Simulated light curves of TIC 287869463 with a Y33+ pulsation mode. The orbital inclination angles in degrees are written next to the right y-axis. The orange curve for i = 75° represents the simulated light curve for the approximate inclination angle of TIC 287869463. The curves for different inclinations are shifted vertically by arbitrary amounts for clarity.
In the text

Fig. 7. Simulated FTs for different examples of ℓ = 3 pulsation modes. The modes represented in the four panels are Y33+, Y33−, Y33x, and Y31x, clockwise starting from the upper left panel. The first two modes are defined in Sect. 3 and Equation (1). The pulsation axis of the latter two modes has been tilted into the orbital plane and lies along the tidal (or x) axis. The FT peaks are arbitrarily shifted to the left for decreasing inclination angles to avoid overlapping peaks.

In the text

	Fig. 8. Diagrams showing the flux perturbation across the surface of a star pulsating in the Y33+ mode (top) and Y33− mode (bottom). The diagram corresponds to orbital phase zero for an observer viewing at at inclination i = 55° relative to the orbital axis. The observed phase of the amplitude modulation of the O1 mode is consistent with that of a Y33+ mode.
In the text

Fig. 9. Top panel: Spectral energy distribution fit for TIC 287869463. The blue curve shows the model spectrum for the primary star (the pulsator), and the red curve represents the model spectrum for the cooler smaller secondary star. The black curve that nearly coincides with the red curve is the sum of the model fluxes. The orange points are measured SED points from the literature (for the SED fitting techniques and the origin of the data points, see Sect. 4). Bottom panel: Light-curve fit using our custom light-curve emulator. The black curve shows the Fourier-reconstructed light curve (i.e., the data), and the red curve shows the model fit (see text for details).

In the text

	Fig. 10. Mode visibility diagram as a function of ℓ. These visibilities were calculated as in Daszyńska-Daszkiewicz et al. (2002) for the BVRI passbands; the red passband matches that of TESS best.
In the text

	Fig. 11. Six different evolution tracks in the HR diagram compared to the measured luminosity and Teff of TIC 287869463. The parameters that are varied in the different tracks involve metallicity Z, αMLT, mass, rotation velocity, and convective overshooting, αov (see text for details).
In the text

Fig. A.1. ETV curves for each component of the three pulsation modes, D1, D2, and O1, as well as from the binary orbital period. All four curves show decreasing periods with time, but not at the same rates. All phases are referenced to the times of a primary eclipse, i.e., at BJD - 2,456,014.4757. Note that the vertical scale for the D1 ETVs is exactly 1/3 that for the D2 and O1 ETVs. In turn, the vertical scale for the orbital ETV curve is 4/5 that of the D1 ETVs.

In the text