© The Authors 2025

1. Introduction

The development of high-precision photometric space missions such as Convection Rotation and Planetary Transits (CoRoT) (Baglin et al. 2006), Kepler (Koch et al. 2010), Transiting Exoplanet Survey Satellite (TESS) (Ricker et al. 2015), CHaracterising ExOPlanets Satellite (CHEOPS) (Benz et al. 2021), and PLAnetary Transits and Oscillations of stars (PLATO) (Rauer et al. 2014) was primarily motivated by the quest for extrasolar planets by detecting transits and for the characterization of these planets. In addition to the orbital inclination and information on the limb darkening, however, the transit method just yields the ratio of the planetary to stellar radius. To determine absolute dimensions, knowledge about the size of the planetary host star is required.

This is where asteroseismology comes into play. Asteroseismology determines the interior structure of stars by using stellar oscillations as seismic waves (e.g., Aerts 2021; Kurtz 2022, and references therein) and is at the same time capable of constraining the absolute parameters of its target stars (e.g., Huber et al. 2013). As asteroseismic investigations themselves also greatly benefit from high-accuracy photometric data, almost all of the planet-search space missions have an asteroseismology program as well (e.g., Gilliland et al. 2010).

The outcome from those missions in terms of asteroseismic research has been rich. To mention just a few results from TESS alone, more than 150 000 oscillating red giant (Hon et al. 2021) and over 100 000 pulsating A–F type stars (Gootkin et al. 2024) have been discovered, period-spacing patterns were detected that permitted asteroseismology of dozens of δ Scuti stars for the first time (Bedding et al. 2020), low-frequency oscillations in massive supergiant stars appear to be ubiquitous (Bowman et al. 2019), and the interiors of individual massive stars were studied (Burssens et al. 2023).

The availability of data of unprecedented quality clearly also facilitated the discovery of unexpected phenomena, such as a star that stopped to pulsate (Kurtz et al. 2025), or stars in close binary systems whose pulsation axes are tilted with respect to their rotation axes through the gravitational pull of their companion star (Handler et al. 2020). Several more of these systems were reported soon afterwards (Kurtz et al. 2020; Rappaport et al. 2021; Kahraman Aliçavuş et al. 2022; Jennings et al. 2024). Whereas the pulsating star in the latter binaries was of the δ Scuti type, a tidally tilted subdwarf B pulsator was also discovered (Jayaraman et al. 2022a). From the point of view of theory, the pulsations of these objects were explained via tidal mode coupling that both tilts the pulsation axes and confines the pulsations to certain regions of the star (Fuller et al. 2020).

More recently, two of the tidally tilted pulsators have been found to be triaxial pulsators with dipole pulsations about three orthogonal axes, including the tidal axis (Zhang et al. 2024; Jayaraman et al. 2024). Based on these discoveries, Fuller et al. (2025) developed a theoretical description of how these modes of pulsation may appear in triaxially deformed stars. Henceforth, these are called tidally tilted standing modes (TTS; sometimes Fuller modes for short). These are modes that arise when tidal, centrifugal, and Coriolis perturbations couple modes of equal spherical degrees (Fuller et al. 2025). These modes were found to be triaxial for dipole configurations, while for higher multipole modes, their geometry is more complex. These modes typically receive nearly equal contributions from spherical harmonics of spherical degrees +m and −m, which means that the modes are primarily standing modes that do not propagate around the rotation axis. Instead, their pattern is nearly stationary with respect to the companion star.

Furthermore, based on this theory, Fuller et al. (2025) predicted how quadrupole triaxial pulsation modes would manifest themselves observationally and provided the theoretical expressions for quadrupole and octupole TTS pulsations. We report the discovery of a textbook example of a pulsator with two dipole modes around different axes in the orbital plane plus a quadrupole Fuller mode that occurs in the primary component of the eclipsing binary system EL CMi.

2. Observations and analysis

2.1. Photometry

The source EL CMi (TIC 321254859) has been observed with TESS in full-frame images taken during Sectors 7 (January 2019), 34 (January 2021), 61 (February 2023), and 88 (January 2025) with cadences of 1800, 600, 200, and 200 s, respectively. TESS-SPOC light curves are available for the first two sectors, and we made our own lightkurve reductions for the Sector 61 and 88 data. All light curves show some drifts in the mean brightness that were filtered out, and we removed segments with poor data quality. EL CMi was initially identified as the high-scoring output of a neural network that was designed to find eclipses in TESS light curves. This is described further in Powell et al. (2021) and Kostov et al. (2025). Upon manual review, tidally tilted or triaxial pulsations were suspected in the Sector 61 data. Consequently, the system was subjected to closer examination regarding the periodic content in its light curves.

To perform frequency analyses, we used the software package Period04 (Lenz & Breger 2005). This package applies Fourier analysis and simultaneous multifrequency sine-wave fitting. Advanced options such as the calculation of optimal light-curve fits for multiperiodic signals including harmonic, combination, and equally spaced frequencies are available. Fits composed of multiple sine waves were subtracted from the data, and the residuals were examined for further periodicities.

TESS data come in different flavors, in particular, simple aperture photometry (SAP; Twicken et al. 2010; Morris et al. 2020) and pre-search data conditioning SAP (PDCSAP; Stumpe et al. 2012; Smith et al. 2012) fluxes. The latter were detrended using co-trending basis vectors and are often cleaner than SAP fluxes. On the other hand, the detrending procedure can distort light curves with intrinsic sharp features, such as eclipses, and the choice of the optimal data version therefore sometimes depends on the goals and demands of the scientific analysis to be carried out.

The eclipses in the Sector 7 and 34 SAP data for EL CMi appear to be shallower than in the PDCSAP and in the Sector 61 and 88 data. This can be explained by blending of nearby stars into the photometric apertures, as implied by the corresponding CROWDSAP parameters (0.8784 and 0.8959). Therefore, the analysis of the binary light curve was based on the Sector 34 PDCSAP data, which are expected to no longer be contaminated by the companion. On the other hand, the Fourier noise level in the range of the pulsation frequencies is considerably lower (by about 30%) in the SAP data than in the PDCSAP light curves. Consequently, the SAP data were used to detect the pulsational frequencies, with the caveat that their amplitudes may be systematically lower than in reality. This is not relevant for the results of the analyses performed on these data, however.

In a first step, the ephemeris for the times of primary minimum from the TESS data was determined by computing a multiharmonic fit to the light curves merged into 1800 s bins to give them all equal weight. This ephemeris was

$$\begin{array}{c}\hfill T_I=2459229.2182(2)+1.05384979(7)E,\end{array}$$

where the zeropoint in time relates to the center of the first primary eclipse in the Sector 34 data. The corresponding orbital phase curve is shown in Fig. 1.

Fig. 1. Phase-folded TESS Sector 34 PDCSAP light curve of EL CMi without the pulsation.

Only the Sector 34, 61 and 88 data were used to analyze the pulsations because the pulsation frequencies are above the Nyquist limit of the Sector 7 data. An essential question in the quest for the detection of periodic signals in a dataset is where to stop the analysis. To this end, several criteria have been proposed. We adopted the method that was originally put forward by Breger et al. (1993). In this framework, a Fourier peak must exceed the mean amplitude in the amplitude spectrum by a factor of 4 in the local frequency domain to be considered significant. For space-based data, this can lead to an overinterpretation of the periodic content, however (Balona 2014). For TESS data, Baran & Koen (2021) have proposed stricter detection criteria depending on the sampling of the dataset. Interpolating their values, we adopted thresholds with a signal-to-noise ratio (S/N) of 4.8 for the Sector 34 data and of 5.0 for the Sector 61 and 88 data. They correspond to false-alarm probabilities of 0.1%. The mean amplitude of the noise was evaluated in a 5 d⁻¹ wide window centered at the frequency under consideration.

We detected the pulsation frequencies based on the residuals of a multiharmonic fit accounting for the eclipses. We then added the pulsation frequencies to a fit that simultaneously included the eclipses and the pulsations. When no significant signals were left, this process was stopped. At this point, an échelle diagram was constructed with respect to the orbital frequency (see Jayaraman et al. 2022a, for more details) to evaluate which of the detected signals form multiplets. The frequencies of confirmed multiplet members were then fixed to their predicted values using the orbital frequency, and the fit to the complete light curve was optimized and re-evaluated to verify that no additional signals had been left behind. The resulting lists of frequencies are given in Table 1, and the resulting échelle diagram is shown in the upper panel of Fig. 2. The pre-whitening process is illustrated in Fig. 3.

Fig. 2. Échelle diagram for EL CMi from Sectors 34 (black), 61 (red), and 88 (blue). Top: All data. Bottom: Primary eclipses removed from the light curve. The sizes of the plot symbols relate to the logarithm of the pulsation amplitude; the cross marks the expected location of the centroid of mode ν2.

Fig. 3. Fourier transforms of the Sector 34 TESS data for EL CMi with subsequent pre-whitening steps (first, only the orbital harmonics, then the orbital harmonics plus 7, 18, and 19 pulsation frequencies indicated by the blue arrows, respectively). The inset in the top panel shows the spectral window function of this dataset. The level of the S/N curve is applicable for all of these amplitude spectra, but is plotted only once for clarity.

The échelle diagram derived in this way is fairly simple, but shows a faint extended ridge around mode ν₃ (Fig. 2, upper panel). This is a sign of spatial filtering (eclipse mapping), which is a known phenomenon among pulsating stars in eclipsing binaries. In brief, when part of the surface of a nonradially pulsating star is covered during an eclipse, only the light variations in the visible part of the star can be measured. This modifies the observed pulsation amplitudes and phases depending on the types of the pulsation modes that are excited (e.g., Van Reeth et al. 2023; Bókon et al. 2025, and references therein).

The effect of spatial filtering occurs when the pulsating star is attenuated by its companion; when the pulsator blocks some light of the constant star, only the background light changes in strength. For this reason, the frequency analysis was repeated using light curves without the primary, and then without the secondary eclipses. As there was no change in the derived set of pulsation frequencies in the latter case, the hotter component of the binary was identified as the pulsator. Figure 4 shows the frequency analysis after the primary eclipses were removed from the Sector 34 data. The list of pulsation-related frequencies derived from all sectors (with sufficient cadence) without the primary eclipses is given in Table 2, and the corresponding échelle diagram can be found in Fig. 2 (lower panel). The results from each individual sector are highly consistent when we acknowledge that the pulsation amplitudes are somewhat variable over time. This precludes the detection of some of the weaker signals in Sectors 61 and 88.

Fig. 4. Same as Fig. 3, but for the data without the primary eclipses. The inset in the top panel shows the spectral window function of this dataset.

With this information in hand, the runs of pulsation amplitude and phase for the three dominant modes were reconstructed using the method described by Jayaraman et al. (2022a) (their Sect. 3.2). Figure 5 shows the results.

Fig. 5. Run of the pulsation amplitude and phase for the three main modes of EL CMi. Black: Sector 34 data. Red: Sector 61 data. Blue: Sector 88 data. Upper panel: Orbital light curve. Second and third panels: All data. Lower two panels: Reconstructions from the data without the primary eclipses. The frequencies are for the multiplet centroids.

Whereas the amplitudes of the three pulsation modes varied from one sector of data to the next, the runs of the pulsational phases over the orbit were essentially identical. The higher amplitude of ν₃ and lower noise level in Sector 34 make the effects of spatial filtering during primary eclipse detectable, whereas there are only indications of this effect in the Sector 88 data, and they are lost in the observational noise in Sector 61. To identify the pulsation modes by comparison with the pulsation amplitude and phase curves theoretically predicted by Fuller et al. (2025), it is more convenient to use the reconstructions from the data without the primary eclipses.

The amplitude and phase modulation of mode ν₁ over the orbit is consistent with that of a Y_10x dipole mode, whereas the run of mode ν₂ is that of a Y_10y dipole mode (see Fuller et al. (2025) for visualizations of these modes). These modes were reported previously (Zhang et al. 2024; Jayaraman et al. 2024). The behavior of the pulsation amplitude and phase depending on orbital phase is clearly different for mode ν₃; it is that of a quadrupole Y₂₂₋ mode, which is thus observationally detected for the first time.

As discussed by Fuller et al. (2025), the Y₂₂₋ mode is a superposition of Y_2 + 2z and Y_2 − 2z modes. This behaves like an l = m = 2 pattern that is wrapped around the z-axis, with maxima and minima offset from the tidal axis by an angle of π/4, and it is a standing mode that does not propagate around the z-axis.

2.2. Spectroscopy

Spectroscopic observations of EL CMi were carried out with the Alhambra Faint Object Spectrograph and Camera (ALFOSC) instrument in spectroscopy mode at the 2.56 m Nordic Optical Telescope (NOT) on the Observatorio del Roque de los Muchachos, La Palma, Spain, during the nights from 2025 January 23/24 to 27/28. Grism 18 was used and resulted in a wavelength coverage of 3450–5350 Å with a spectral resolution of 2000 using a 0.5″ slit. Integration times of 300 s were used and resulted in a S/N of 65–145 around 4500 Å, depending on weather conditions. Each observation consisted of two or three consecutive integrations.

The data were reduced with the package PypeIt (Prochaska et al. 2020a,b), using standard day-time calibrations as well as arc frames taken immediately after the on-target integrations in order to avoid spurious variability due to flexure, and applying barycentric corrections.

Radial velocities (RVs) were extracted from these data using the Image Reduction and Analysis Facility (IRAF) task fxcor. To facilitate the cross-correlation, we computed a synthetic spectrum with T_eff = 7750 K, log g = 4.0, and v sin i = 90 km s⁻¹ with the code SPECTRUM (Gray & Corbally 1994). The radial velocities so obtained are listed in Table 3. Using the ephemeris in Sect. 2.1 and assuming a circular orbit yields a RV amplitude of K₁ = 102.3 ± 3.0 km s⁻¹ for the primary star and a systemic velocity γ = +28.1 ± 3.0 km s⁻¹. The stated uncertainties include the formal error of the fit and the wavelength calibration error of the spectrograph. Figure 6 shows the RV measurements in comparison to this and a model fit to be derived in the next section.

Fig. 6. Orbital NOT radial velocities (black points) with the simple Keplerian (red line) and the PHOEBE (blue line) fits that also include the predicted Rossiter-McLaughlin signature (Rossiter 1924; McLaughlin 1924). The uncertainties on the measurements are about the size of the plot symbols.

Finally, the individual spectra were corrected for RV shifts, were coadded and compared to model spectra of different temperatures and surface gravities. The best fit was found for T_eff = 7750 K, log g = 4.0, v sin i = 90 km s⁻¹ (hence the previous choice of template). The primary star appears to be somewhat underabundant in heavy elements compared to the Sun ([M/H]= − 0.3 ± 0.2). Our spectra are insufficient to determine whether this is a true underabundance or if the primary star has some mild chemical peculiarity of the λ Bootis type, however.

3. System parameters

3.1. Data availability

The Gaia DR3 catalog (Gaia Collaboration 2016, 2023) gives a G magnitude of 11.82 and (B_P − R_P) = 0.34, which correspond to the mean values out of eclipse. The Gaia DR3 parallax of the system is 0.977 ± 0.019 mas. The reddening maps of Green et al. (2019) imply E(g − r) = 0.025 ± 0.015, consistent with Schlegel et al. (1998) (A_g = 0.139), which was taken as an upper limit to the interstellar reddening of the system given its Galactic latitude of +15 deg.

With the reddening relations by Zhang & Yuan (2023) and the transformations from the Gaia photometric system given in Table 5.9 of the Gaia DR3 release notes, we derived E(B_P − R_P) = 0.030 ± 0.018 and V = 11.86. The standard relations by Pecaut & Mamajek (2013) then suggest T_eff = 7670 K, which agrees very well with the spectroscopic estimate above. With a bolometric correction of BC = +0.015, A_v = 0.08 and the Gaia parallax M_bol = +1.74 ± 0.07 were derived (the error budget is dominated by the uncertainties in interstellar extinction and the V magnitude). With the bolometric magnitude of the Sun (Torres 2010), M_bol, ⊙ = 4.74, a mean out-of-eclipse system luminosity of L = 15.8 ± 1.1 L_⊙ results.

Archival spectral flux measurements from 0.15 to 11.6 microns were collected. These data were taken from VizieR SED¹ (Ochsenbein et al. 2000), which uses a systematic sky coverage of surveys such 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) (Wright et al. 2010), and Galaxy Evolution Explorer (GALEX) (Bianchi et al. 2017). For EL CMi, we found 26 spectral points that we were able to use.

3.2. General approach

In what follows, the EL CMi system parameters are derived using two different approaches. In the first approach, the PHOEBE2 software is used to analyze the eclipsing light curve obtained with TESS plus the RV data from the primary star. In the second approach, our own code for fitting the spectral energy distribution (SED) coupled with a relatively simple eclipsing light-curve emulator is used to jointly analyze the TESS light curve, the available SED data, and the radial velocities for this source. Each approach is detailed in the following subsections.

3.3. PHOEBE analysis

The new-generation Wilson-Devinney code PHOEBE2 (Prša et al. 2016) was used to analyze the TESS eclipsing light curve of EL CMi and to simultaneously take into account the measured RV data points of the primary star. A Monte Carlo Markov chain (MCMC; Foreman-Mackey et al. 2013) wrapper was used to evaluate the uncertainties in the fitted parameters by sampling over a broad range of uninformed priors on the ratio of effective temperatures, primary radius, primary mass, mass ratio, secondary radius, orbital inclination, and systemic velocity.

The fitted light curve is shown in Fig. 7. The best-fitting parameters that result from this fit are shown in the first results column of Table 4. We simply fixed the T_eff of the primary at 7750 K, as deduced from the spectroscopy (see Sect. 2.2). Both stars were modeled with Castelli et al. (2003) model atmospheres, where these atmospheres were used to derive the emergent intensity of each stellar surface element as well as the limb darkening via the PHOEBE2 interpolated limb-darkening scheme, and with a default gravity brightening (β = 0.32) and bolometric albedos (α = 0.6).

Fig. 7. Upper panel: Fit to the phase-binned TESS PDC light curve from Sector 34 (black points) compared to the fit derived with PHOEBE2 (blue line) and the combined SED + light-curve fitting approach (red line, data and fit shifted by 10% in intensity for clarity). Lower panel: Residuals from these fits. Blue: PHOEBE2, and red: SED + light-curve fitting, shifted by 4% in intensity for clarity.

The mass and T_eff of the pulsating δ Scuti star are typical for these pulsators. On the other hand, the secondary, with a mass of only 0.93 M_⊙ and radius of nearly 2 R_⊙, strongly indicates a star that has lost its envelope during a prior episode of mass transfer. This is discussed further in the next section.

3.4. Light curve plus SED analysis

The SED of a binary can yield important information on the two stellar radii and the two T_eff values when the distance to the source is known accurately. The SED by itself is not generally sufficient to yield all four of these parameters, however. Additionally, the eclipsing light curve contains key information on quantities such as the ratio of the two T_eff values and the ratios of the stellar radii to the semi-major axis. The latter requires the system mass in order to translate R₁/a and R₂/a into information about the physical stellar radii, however. In turn, radial velocities can provide information about the masses. In any case, supplemental information about the radii and effective temperatures can indeed facilitate the inference of a number of the stellar properties from the SED data.

Under the assumption that the two stars in the binary have evolved coevally (i.e., without a prior episode of mass transfer), the SED can be fit with just four basic parameters: M₁, M₂, system age, and an interstellar extinction, A_V. If the metallicity is not known from the spectral observations, it can be added as a fifth parameter.

The mass and age yield the radius and T_eff of the two stars via stellar evolution tracks, and the SED can then be fit, provided that there is some other supplemental information such as the ratios of radii or T_eff from a light-curve analysis, for instance, from PHOEBE2. Since PHOEBE2 does not have a built-in SED fitter, the light-curve fitting with PHOEBE2 and the SED fitting could be iterated with separate codes. We instead chose to write a combined SED fitter with a simplified light-curve emulator that can be used to jointly fit for both the SED and the eclipsing light curve. The code follows the description by Jayaraman et al. (2024). Briefly, the simplified light-curve emulator uses (i) spherical stars with limb darkening, and (ii) sets of sines and cosines to represent the out-of-eclipse behavior, including ellipsoidal light variations, illumination effects, and Doppler boosting (see, e.g. Carter et al. 2011).

The SED portion of the fitting code uses the MESA Isochrones and Stellar Tracks (MIST) evolution tracks (Choi et al. 2016; Dotter 2016; Paxton et al. 2011, 2015, 2019). For any given mass, age, and metallicity, the tracks provide the stellar radius and T_eff. The stellar atmosphere models were taken from Castelli et al. (2003). The SED fitting portion of the code has been extensively tested on a number of multistellar systems, including binaries (Jayaraman et al. 2022b), triples (Rappaport et al. 2022), and even a sextuple star system (Powell et al. 2021). For EL CMi, metallicities of [Fe/H] = 0.0 and –0.3 were tested in keeping with what was learned from the stellar spectra (see Sect. 2.2). There was not much difference in the results depending on the metallicity.

This code was used to jointly fit the light curve, SED, and RV data for EL CMi, in agreement with the PHOEBE2 result. We found no viable solutions in which the two stars evolve in a coeval manner. Thus, there must have been a prior episode of mass transfer onto the current primary star (the pulsator). We therefore changed the fitted parameters to the following six quantities: primary and secondary mass (M₁ and M₂), radius and T_eff of the secondary star (R₂ and T_eff, 2), the age of the primary star (since its rejuvenating accretion episode, τ), and A_V.

The fitted light curve is shown in Fig. 7. Figure 8 shows a comparison of the observed to theoretical SED derived above as well as the fitted light curve. The agreement is satisfactory. Only the point in the far ultraviolet is somewhat off. This may be due to the wide bandwidth of the GALEX FUV filter and the steep continuum slope in this wavelength range, however (see Jayaraman et al. 2024).

Fig. 8. Fit to the SED of the EL CMi system. The orange dots with the error bars show the observational values. The blue line is the SED of the primary star, the red line shows the SED of the secondary, and the black line shows the sum of the two.

The parameters of the binary system are given in two columns of Table 4 along with some other derived parameters in Table 5. The second column of the results parameters in Table 4 lists the case without prior restrictions on the radius of the secondary star, except that it cannot overfill its Roche lobe. The third results column in Table 4 lists the restricted case when the secondary star always just fills its Roche lobe. This is motivated by evolutionary studies (see Sect. 4), which indicated that after the episode of rapid mass transfer, the donor star remains close to filling its Roche lobe for 1.3 Myr after the bulk of the mass transfer has ceased.

A comparison of the three results columns in Table 4 paints a consistent picture of the stellar parameters for EL CMi. The main difference in the three results columns arises in the radii, depending on whether the secondary is assumed to fill its Roche lobe. In the former case, the primary is 24%±11% larger than the secondary, while in the latter case, the two stars are more nearly equal in size, with the secondary being slightly larger.

4. System evolution

As the data analysis suggests that the radius of the secondary is too large to be explained as a product of a single-star evolution, it is logical to suspect that EL CMi is a product of a system that has passed through a mass-transfer phase. To model the system as the result of a prior episode of mass transfer, we used the MESA code (Modules for Experiments in Stellar Astrophysics, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023, version 23.05.1) with the MESA-BINARY module under the assumption of an implicit mass-transfer rate.

To find a model that fits the parameters of the system, in particular the masses, radii, effective temperatures, and the orbital period, we built a grid of binary evolution tracks for initial periods spanning a range from 0.6 to 10 d and for initial masses in the range of 1.0–3.0 M_⊙, with a constant step size of 0.2 d and 0.2 M_⊙ respectively. Additionally, we used the results from Section 3 and we imposed the constraints that the initial total mass of the system should not be lower than 3.0 M_⊙, but no higher than 3.5 M_⊙ to account for some small mass loss. After obtaining an initial fit, we refined the grid resolution around the best-matching model by reducing the steps to 0.02 M_⊙ in mass and 0.05 d in orbital period.

Four different values of metallicity, Z = 0.005, 0.010, 0.014, and 0.020, corresponding to [M/H] = −0.47, −0.16, −6.5 × 10⁻³, and 0.16, respectively, were used and scaled using the linear scaling relations by Choi et al. (2016) to determine the initial abundances of hydrogen and helium. We used the AGSS09 (Asplund et al. 2009) initial chemical composition of the stellar matter and the OPAL opacity tables supplemented with data from Ferguson et al. (2005) for lower temperatures. We treated convective instability using the Ledoux criterion, combined with the mixing length theory based on the model by Henyey et al. (1965). In regions that are stable according to the Ledoux criterion but unstable by the Schwarzschild criterion, we used semiconvective mixing with a scaling factor of α_sc = 0.1, following the formalism of Langer et al. (1985). To mix regions in which the mean molecular weight was inverted, such as those that formed during mass accretion, we applied thermohaline mixing using the formalism of Kippenhahn et al. (1980) with an α_th = 1 coefficient. As these models include neither rotation nor rotationally induced mixing, we used a minimum diffusive mixing coefficient of D = 10 cm² s⁻¹, which also allowed us to smooth out numerical noise or discontinuities in internal profiles.

The values of free stellar structure parameters were fixed according to the following theoretical models and asteroseismic calibrations. Given the primary mass, the convective-core overshooting parameter was set to f_ov = 0.02 for both components, in agreement with the values predicted by Claret & Torres (2016). Furthermore, the mixing length parameter and the fraction of mass lost from the system during mass transfer were set to α_MLT = 1.5 and β = 0.2, in agreement with asteroseismic determinations for δ Scuti stars in post-mass-transfer binaries (e.g. Miszuda et al. 2022).

From each set of binary evolutionary tracks, we selected models that agreed with the measured current-epoch orbital period, and we calculated a normalized χ² discriminant using observed and model masses, radii, and effective temperatures for these models. This resulted in a set of 40 models that fit the observations within 1σ and cluster around P_orb, 0 ∈ (1.0, 1.05) d, M_don, 0 ∈ (1.9, 2.06) M_⊙ and M_acc, 0 ∈ (0.9, 1.0) M_⊙ for metallicities Z = 0.014 and Z = 0.02. A visualization of the best-fitting model is presented in Figure 9.

Fig. 9. Hertzsprung-Russel diagram and time evolution of the stellar masses, radii, effective temperatures, Roche-lobe filling factor, and the orbital period for the binary evolution track that best reproduces the observed parameters of the system. This track was computed for an initial orbital period Porb, 0 = 1.05 d, primary and secondary masses Mdon, 0 = 1.94 M⊙ and Macc, 0 = 0.94 M⊙, metallicity Z = 0.014, overshooting parameter fov = 0.02, mixing-length parameter αMLT = 1.5, and mass-transfer efficiency β = 0.2. The evolutionary tracks of the accretor and the donor are shown with orange and blue lines, respectively, and the same colors are used to indicate the observed ranges of the corresponding parameters, which are marked with shaded bands in the plots. The model at the time corresponding to the observed orbital period is marked with vertical dashed lines in the lower panels and with circular markers in the HR diagram.

According to these models, the system is approximately 1.7 ± 0.3 Gyr old and is still undergoing mass transfer. The mass-transfer rate is currently fairly low at Ṁ ≈ 10⁻¹⁰ M_⊙ yr⁻¹. The fastest phase of mass transfer began when the system was ≈1 Gyr old, and it peaked at Ṁ ≈ 10⁻⁶ M_⊙ yr⁻¹. At this time, the donor star was in the later part of its main-sequence evolution, with a central hydrogen abundance of X_c = 0.26, while the accretor had barely evolved from the zero-age main sequence. This suggests that the binary went through type A mass transfer (MT) (Kippenhahn & Weigert 1967). Currently, after rejuvenation of the accretor and mass-ratio reversal, the donor is close to exhausting its central hydrogen (X_c < 0.1) and has retained a sufficiently thick hydrogen envelope to start ascending the red giant branch soon. On the other hand, the accretor is currently evolving through the main sequence, with a central hydrogen abundance of X_c = 0.45. These parameters are summarized in Table 6.

Despite the decline in the mass-transfer rate following the initially faster phase, the donor remains close to Roche-lobe filling over an extended timescale, with R_don/R_RL ≲ 1. Although the donor is initially driven out of thermal equilibrium, it quickly recovers equilibrium after the mass ratio reverses. From this point onward, the donor becomes the less massive star, and its radius evolution, governed by hydrogen burning, begins to regulate the mass-transfer rate. As a result, the subsequent mass transfer proceeds on the nuclear timescale of the donor, allowing the stable low-rate Roche-lobe overflow configuration to persist until the central hydrogen is exhausted. During this phase, the donor maintains sufficient pressure support to avoid contraction, and the system evolves stably under slow mass exchange (e.g., Kolb & Ritter 1990).

Initially, the accretor temporarily departs from thermal equilibrium due to the rapid accretion. This results in an increased luminosity. The accretor regains equilibrium shortly after, however, and continues its main-sequence evolution with a steadily growing mass and radius. When it reaches the end of the mass-transfer phase, the accretor will attain a mass of 2.16 M_⊙, cross the Hertzsprung gap, and become a red giant.

As mass transfer continues, the donor gradually loses its hydrogen-rich envelope. When the mass transfer ceases, the envelope is no longer massive enough (M = 0.2 M_⊙) to maintain efficient shell hydrogen burning or to provide the pressure required to compress the core toward helium ignition. At this point, the central temperature only reaches log T_c ≈ 7.52, and the helium core mass is limited to M_core ≈ 0.17 M_⊙. These values are both well below the thresholds required for helium ignition. Consequently, instead of ascending the RGB up to helium ignition, the star will evolve toward the pre-helium white dwarf (pre-He WD) phase. Despite its relatively high initial mass, the donor ultimately avoids helium ignition and becomes a low-mass helium white dwarf star.

At an age of ≈2.9 Gyr, the accretor star evolves into an AGB star and starts mass transfer to the helium white dwarf. At this point, the AGB star has a helium core mass of 0.488 M_⊙ and a carbon/oxygen core of 0.280 M_⊙. Since the binary mass ratio is high, the mass transfer is dynamically unstable, and the binary system will enter a common envelope (CE) phase. By comparing the binding energy of the envelope with the orbital energy of the system, we find that the CE cannot be ejected even for a CE ejection efficiency of 1.0. Consequently, the binary system ultimately merges into a single star.

It is unclear whether the ongoing accretion in the current evolutionary state of the system affects the triaxial pulsations of the primary (the accretor), but we do not expect strong effects. Because the primary is underfilling its Roche lobe and its current accretion rate is low, we do not expect substantial deviations from the hydrostatic equilibrium that is assumed in the triaxial pulsation model. The accretion power is L_ac ≈ GMṀ/R ≈ 0.003 L_⊙ (Ṁ/10⁻¹⁰ M_⊙/yr), which is much lower than the luminosity of the primary. Its structure and pulsations are therefore unlikely to be greatly modified.

5. Discussion and conclusions

The first objects called “tidally tilted pulsators” were explained as stars whose pulsational axis was pulled into the orbital plane and pointed at their close companion star at each point in the orbit. The first such pulsator, HD 74423, had a single dipole oscillation mode whose amplitude was also enhanced in one hemisphere. This motivated the term “single-sided pulsator” (Handler et al. 2020). The second such object, CO Cam, showed multiperiodic axisymmetric pulsations with a higher amplitude in the direction of the companion (Kurtz et al. 2020), and the third object, TIC 63328020, appeared to be distinctly different with its tesseral tidally tilted dipole mode (Rappaport et al. 2021).

The idea that stars can pulsate triaxially was based on the observational studies of TIC 184743498 by Zhang et al. (2024) and of TIC 435850195 by Jayaraman et al. (2024), which demonstrated that both objects have very many tidally tilted dipole modes. The geometries of these binary systems, and hence the inclination of the tilted stellar pulsational axis to the line of sight, however, made an explanation of these objects in terms of simple tidally tilted pulsators impossible because the observed amplitude and phase modulations of the individual pulsations during the orbit are predictable depending on the type of the pulsation modes and on the orbital geometry of the system (Reed et al. 2005). The behavior of the pulsations of these latter two stars did not conform to a simple picture of tidal tilting, but could be explained with pulsation around three different axes.

The theoretical foundation for this interpretation was provided by Fuller et al. (2025), who showed that the different azimuthal orders of a dipole frequency multiplet are tidally coupled in these systems. Thus, they align with the three principal axes of the triaxially deformed pulsating star. Based on this theory, theoretical predictions of the manifestations of nonradial quadrupole oscillation modes were provided. We presented the first observational detection of such a quadrupole mode.

The primary component of EL CMi is the third object in which pulsation modes that only occur in triaxially distorted stars reported. It is therefore logical to suspect that the previously discovered tidally tilted or single-sided pulsators may be just special cases of triaxial pulsators. After all, all of these objects are distorted into triaxial shapes, and their pulsation modes conform to the same framework. The single mode of HD 74423 is consistent with a tidally enhanced Y_10x mode, whereas the four modes of CO Cam appear to be of the same type. The behavior of the single mode of TIC 63328020, which was originally classified as a tidally tilted l = 1, m = 1 (Y_11x) mode, is also consistent with a Y_10y mode. The numerous mode doublets and one triplet in TZ Dra (Kahraman Aliçavuş et al. 2022) are also consistent with Y_10x and Y_10y modes, respectively, and the pulsation modes of the subdwarf B star HD 265435 (Jayaraman et al. 2022a) conform to the predictions of the behavior of modes in triaxially distorted stars by Fuller et al. (2025) as well.

The evolutionary state of the tidally tilted or triaxial pulsators with a main-sequence star as the primary components show for EL CMi, TIC 63328020 and TZ Dra that they have undergone or are undergoing mass transfer, whereas HD 74423, CO Cam, TIC 184743498, and TIC 435850195 are detached systems. The appearance of these special types of pulsation therefore appears to be not connected to a particular evolutionary phase of close binary systems.

The pulsating star in EL CMi is a very simple and clear example that demonstrates that the main features of the theoretical description of Fuller modes are consistent with observations. It is reasonable to expect that stars with richer, more complicated pulsation spectra will be discovered in the future. These pulsations in essence come with an observational identification of the underlying modes and occur in stars that can be characterized with a high degree of confidence and accuracy because they are located in binary and multiple systems. These objects might therefore be used for detailed asteroseismic studies based on oscillations that occur around three different symmetry axes. In other words, these objects hold the promise of enabling us to study the interior structure of distant stars in three dimensions.

Data availability

We make all files needed to recreate our MESA-BINARY results publicly available at Zenodo: https://zenodo.org/records/15745684

Acknowledgments

GH and AM both thank the Polish National Center for Science (NCN) for supporting this study through grant 2021/43/B/ST9/02972 and Filiz Kahraman Aliçavuş for her comments on the spectrum of the target. DJ acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant “Nebulosas planetarias como clave para comprender la evolución de estrellas binarias” and the European Regional Development Fund (ERDF) with reference PID-2022-136653NA-I00 (DOI: 10.13039/501100011033). DJ also acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant “Revolucionando el conocimiento de la evolución de estrellas poco masivas” and the European Union Next Generation EU/PRTR with reference CNS2023-143910 (DOI: 10.13039/501100011033). H-LC is supported by the National Key R&D Program of China (grant Nos. 2021YFA1600403), the National Natural Science Foundation of China (grant Nos. 12288102, 12333008 and 12422305). VK acknowledges support from NSF grant AST-2206814. 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 (https://dx.doi.org/10.17909/t9-r086-e880), 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 (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. This research has also made use of the VizieR catalogue access tool, CDS, Strasbourg, France. Software: This paper made use of the following codes/packages: LIGHTKURVE (Lightkurve Collaboration 2018), IRAF (Tody 1986) SPECTRUM (Gray & Corbally 1994) PERIOD04 (Lenz & Breger 2005) MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023).

References

All Tables

All Figures

	Fig. 1. Phase-folded TESS Sector 34 PDCSAP light curve of EL CMi without the pulsation.
In the text

	Fig. 2. Échelle diagram for EL CMi from Sectors 34 (black), 61 (red), and 88 (blue). Top: All data. Bottom: Primary eclipses removed from the light curve. The sizes of the plot symbols relate to the logarithm of the pulsation amplitude; the cross marks the expected location of the centroid of mode ν2.
In the text

Fig. 3. Fourier transforms of the Sector 34 TESS data for EL CMi with subsequent pre-whitening steps (first, only the orbital harmonics, then the orbital harmonics plus 7, 18, and 19 pulsation frequencies indicated by the blue arrows, respectively). The inset in the top panel shows the spectral window function of this dataset. The level of the S/N curve is applicable for all of these amplitude spectra, but is plotted only once for clarity.

In the text

	Fig. 4. Same as Fig. 3, but for the data without the primary eclipses. The inset in the top panel shows the spectral window function of this dataset.
In the text

	Fig. 5. Run of the pulsation amplitude and phase for the three main modes of EL CMi. Black: Sector 34 data. Red: Sector 61 data. Blue: Sector 88 data. Upper panel: Orbital light curve. Second and third panels: All data. Lower two panels: Reconstructions from the data without the primary eclipses. The frequencies are for the multiplet centroids.
In the text

	Fig. 6. Orbital NOT radial velocities (black points) with the simple Keplerian (red line) and the PHOEBE (blue line) fits that also include the predicted Rossiter-McLaughlin signature (Rossiter 1924; McLaughlin 1924). The uncertainties on the measurements are about the size of the plot symbols.
In the text

	Fig. 7. Upper panel: Fit to the phase-binned TESS PDC light curve from Sector 34 (black points) compared to the fit derived with PHOEBE2 (blue line) and the combined SED + light-curve fitting approach (red line, data and fit shifted by 10% in intensity for clarity). Lower panel: Residuals from these fits. Blue: PHOEBE2, and red: SED + light-curve fitting, shifted by 4% in intensity for clarity.
In the text

	Fig. 8. Fit to the SED of the EL CMi system. The orange dots with the error bars show the observational values. The blue line is the SED of the primary star, the red line shows the SED of the secondary, and the black line shows the sum of the two.
In the text

Fig. 9. Hertzsprung-Russel diagram and time evolution of the stellar masses, radii, effective temperatures, Roche-lobe filling factor, and the orbital period for the binary evolution track that best reproduces the observed parameters of the system. This track was computed for an initial orbital period Porb, 0 = 1.05 d, primary and secondary masses Mdon, 0 = 1.94 M⊙ and Macc, 0 = 0.94 M⊙, metallicity Z = 0.014, overshooting parameter fov = 0.02, mixing-length parameter αMLT = 1.5, and mass-transfer efficiency β = 0.2. The evolutionary tracks of the accretor and the donor are shown with orange and blue lines, respectively, and the same colors are used to indicate the observed ranges of the corresponding parameters, which are marked with shaded bands in the plots. The model at the time corresponding to the observed orbital period is marked with vertical dashed lines in the lower panels and with circular markers in the HR diagram.

In the text