TIC 184743498: The First Tri-Axial Stellar Pulsator

We have discovered a $\delta$ Scuti pulsator in a tight binary (P = 1.053 d) with nine pulsation modes whose frequencies are between 38 and 56 d$^{-1}$. Each of these modes exhibits amplitude modulations and $\pi$-rad phase shifts twice per orbital cycle. Five of these modes exhibit amplitude and phase shifts that are readily explained by dipole pulsations along an axis that is aligned with the binary's tidal axis. The novelty of the system lies in the remaining four pulsation modes, which we show are dipole pulsations along an axis that is perpendicular to both the tidal axis and the binary's orbital angular momentum axis. There are additionally two pulsation modes whose amplitudes and phases do not change significantly with orbital phase; they are explained as dipole modes along an axis aligned with the orbital/rotation axis. Hence, we propose that TIC 184743498 is a tri-axial pulsator, the first of its kind.


INTRODUCTION
Studies of the photometric variability of stars due to stellar pulsations have been carried out for over a century.For instance, the variability of the bright early F-type star, Delta Scuti (δ Sct), which lends its name to an entire class of pulsating stars, was discovered in radial velocities by Campbell & Wright (1900).They found a range of about 10 km s −1 in 7 measurements made from 1899 to 1900.Fath (1935) and Colacevich (1935) noted sinusoidal photometric variability in δ Sct stars with maximum light at minimum radial velocity over a period of 0.193 d, and concluded that the star is not a spectroscopic binary.Baglin et al. (1973), Breger (1979Breger ( , 2000)), and Kurtz (2022) summarize more recent work on this class of star.
δ Sct stars, along with almost all other classes of pulsating variables, are observed to have nonradial modes.Unlike the symmetrical radial pulsations of the dominant modes of Cepheid variables (as described in, e.g., Eddington 1926), nonradial modes have a pulsation axis.Originally, it was implicitly assumed that a star's ⋆ E-mail: vzhang25@andover.edupulsation axis coincides with its rotation axis, as it is about this axis that most stars deviate from spherical symmetry.It is still the case in modern asteroseismology (Aerts, Christensen-Dalsgaard & Kurtz 2010) that the starting point for modeling and inference is the assumption of the alignment of the pulsation and rotation axes.
However, stars can be distorted from spherical symmetry by other effects, primarily the presence of strong, global magnetic fields, and by tides arising from close companions.Therefore they may have pulsation axes not aligned with the rotational axis.The first stars found to pulsate about an axis other than the rotation axis were the rapidly oscillating Ap (roAp) stars (Kurtz 1982), which have primarily dipole magnetic fields that are inclined to the rotation axis.The combined effects of the rotation and magnetic fields lead to a pulsation axis that is inclined to the rotation axis and possibly offset from the stellar center.As the star rotates, the pulsation mode is then seen from varying 'latitudinal view angles' 1 , Θ, lead-1 Consider a star with a spin axis defined by ŝ and a pulsation axis, p lying at an obliquity angle β with respect to the spin axis, and is viewed by a distant observer along a direction v, such that v • ŝ ≡ cos i, and i is the inclination angle.We define an angle Θ, the latitudinal view of the pulsation axis, ing to modulation of both the observed pulsation amplitude and phase.This phenomenon is known as oblique pulsation.Kurtz (1982) introduced the oblique pulsator model for the roAp stars, which has been refined over the last 40 years for asteroseismic inference.Holdsworth et al. (2021) provides the most upto-date analysis of a large sample of roAp stars that were observed during the first two years of the Transiting Exoplanet Survey Satellite's (TESS) mission (Ricker et al. 2015).A natural extension of the oblique pulsator model is to stars in close binary systems, wherein the pulsation axis could align with the tidal axis (the line joining the two stars).Since δ Sct stars are common upper main-sequence pulsators that are often found in close binary systems, it was no surprise that such "tidally tilted pulsators" (TTPs) were found in TESS data.
The first such system that exhibited TTPs, HD 74423, is a pair of nearly identical chemically peculiar λ Boo stars that nearly fill their Roche lobes in a 1.58-d binary (Handler et al. 2020).Only one of the stars is pulsating and in only one mode, which is unusual for a δ Sct star.The pulsation mode in HD 74423 is a highly distorted dipole mode that is, remarkably, largely confined to one hemisphere of the star.This discovery provided the impetus for the theoretical groundwork for understanding the interaction of pulsation and tidal distortion in close binary stars (Fuller et al. 2020).Such pulsations can provide significantly more information about the mode geometry than can be gleaned from 'normal' nonradial modes, and allows for the identification of the mode degree, ℓ, and azimuthal order, m (see, e.g., Reed, Brondel & Kawaler 2005).
The second TTP found was CO Cam, a marginal Am star in a 1.27-d binary with an undetected companion.It pulsates in at least four tidally tilted modes (Kurtz et al. 2020).Then the third discovery was TIC 63328020, which showed, for the first time, a tidally tilted sectoral (|m| = ℓ) dipole mode (Rappaport et al. 2021).Previously, all identified modes in roAp stars had been zonal (m = 0) distorted dipole and quadrupole modes; this discovery allowed for further study of non-axisymmetric tidally tilted pulsations.The fourth discovery was HD 265435, a subdwarf B star which has been stripped of its H-rich envelope; this star has over 30 modes, of which at least 25 are tidally tilted (Jayaraman et al. 2022).
There has also been recent interest in less extreme examples of this phenomenon (the "tidally perturbed pulsators," see, e.g., Southworth et al. 2020;Van Reeth et al. 2023;Johnston et al. 2023).Because our understanding of mode excitation, mode selection, and mode trapping in all pulsating stars remains rudimentary, discoveries of tidally tilted (and perturbed) pulsators with more, and different, kinds of pulsation modes provides further constraints on developing theory and asteroseismic inference.
However, the oblique pulsator model has been unable to explain certain observations of roAp stars.For instance, Kurtz et al. (2011) concluded that the two pulsation modes of the roAp star KIC 10192926 had different pulsation axes, both of which were inclined to the rotation axis.For another roAp star HD 6532, Kurtz & Holdsworth (2020) found that the single pulsation mode appeared to have a different pulsation axis when observed through the red TESS filter, compared to the pulsation axis derived from groundbased B-band observations.Because these two filters sample difas −p• v.If p revolves about ŝ by an angle ϕ due to either the rotation of the star or the orbital motion in the case of a tidally tilted pulsation axis, then it is straightforward to show that: cos Θ = cos i cos β + sin i sin β cos ϕ.If β = 0, i.e., no obliquity, then Θ = i, a constant.ferent atmospheric depths, perhaps the mode in this star has a complex three-dimensional geometry; this result may imply that different pulsation modes in a single star may be trapped in mode cavities with different pulsation axes.This effect remains unstudied.
In this work, we report the discovery of TIC 184743498, a close binary system whose primary star exhibits δ Sct pulsations over the frequency range 40−55 cycles d −1 .This star exhibits nonradial modes with pulsation axes aligned with the rotation axis, the tidal axis, and a third axis perpendicular to both of those.Thus, it is a tri-axial pulsator, the first of its kind.First, Section 2 describes how this source was discovered.Then, Section 3 discusses the data we used in addition to the photometry from TESS.In Section 4, we use the available observational data on this system to make robust estimates of the stellar parameters of the two binary stars, as well as a tertiary star, in TIC 184743498.The pulsation spectrum is discussed in Section 5, and the amplitudes and phases of these pulsations are reconstructed as a function of orbital phase in Section 6.The expected density of radial pulsation modes for TIC 184743498 is estimated in Section 7. We present in Section 8 a rudimentary perturbation model that can explain in a natural way how this star could have three different orthogonal pulsation axes.We summarize our results in Section 9, and look at future applications of our findings.

DISCOVERY
The first two TTPs discovered (Handler et al. 2020;Kurtz et al. 2020) were found serendipitously by visual examination of millions of TESS light curves from full frame images (Kristiansen et al. 2022).Since then, we have focused on TESS targeted stars with 120-s cadence calibrated light curves generated by the Science Processing Operations Center (SPOC; Jenkins et al. 2016).One of us (RJ) carried out a preliminary search for periodicities in all of the approximately 20,000 SPOC sources processed for each sector.A single summary sheet is generated for the approximately 1,000 stars from each sector that exhibit significant periodicity.The summary sheet includes plots of the raw light curve, the Fourier transform (FT) with various scalings, the inferred principal period, and a light curve folded about this period.The summary sheets also provide properties for each source taken from the TESS Input Catalog (TIC v8.0; Stassun et al. 2019), including the effective temperature, sky location, and estimated radius.Some of us then visually examine these summary sheets to search for interesting periodic stellar phenomena, including stellar pulsations in binary systems.
As a part of a continuing effort to identify TTPs, VZ visually inspected summary sheets for sources from the recent TESS Sectors 56 -62, specifically looking for stars in binary systems whose FT displayed peaks at higher frequencies (typically ≳ 15 d −1 ).During visual inspection of the summary sheets from Sector 62, she found an interesting set of pulsations in TIC 184743498.In particular, she noticed high frequency pulsations between 40 − 55 d −1 in a ∼ 1 d eclipsing close binary.Fig. 1 shows a small portion of the TESS light curve for TIC 184743498 covering about 3 d during Sector 62. Shallow eclipses, ellipsoidal light variations (ELV), and prominent pulsations are salient features of the raw light curve.
The top panel of Fig. 2 shows the Fourier amplitude spectrum of data from TESS Sectors 61 and 62 before any filtering.The FT clearly displays numerous orbital harmonics and a rich set of pulsations from TIC 184743498.The lower panels show the FT after the orbital harmonics have been removed, and then after the largest amplitude pulsations have been cleaned out.This 'cleaning' process for the Fourier transform is explained later in Section 5.

DATA AVAILABILITY
3.1 TESS Data TIC 184743498 was observed by TESS at 30-min cadence in Sector 8 (from 2019 May 16 to June 13), at 10-min cadence in Sectors 34 and 35 (2021 January 14 to March 6), 200-s cadence in Sector 61 (2023 January 18 to February 12), and 2-min cadence in Sector 62 (2023 February 12 to March 10).This target was selected for 2-min cadence observations because it was part of the TESS Candidate Target List (CTL v8.01), a collection of bright (T < 13) stars in the TESS field of view that could host exoplanets.We accessed the data using the Python package lightkurve (Lightkurve Collaboration et al. 2018).We used the Sector 61 and 62 data (200-s and 2-min cadence) for our principal pulsation analysis, but also utilised the Sector 34 and 35 data to confirm our findings based on the later sectors, and to better determine the orbital period.The Sector 61 and 62 data, when combined, span 51 d and contains 28 700 data points.

Other Data
We also made use of archival data from Gaia Data Release 3 (Gaia Collaboration et al. 2023), the Mikulski Archive for Space Telescopes (MAST) 2 , the All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014;Kochanek et al. 2017), and the online VizieR SED viewer (Ochsenbein, Bauer & Marcout 2000) 3 .The basic archival parameters of TIC 184743498 from Gaia and MAST are summarized in Table 1.We used the ASAS-SN data primarily to obtain an accurate and independent value of the orbital period over a baseline of ∼10 yr.The spectral energy distribution (SED) data were obtained from VizieR; these points were used in conjunction with other data to estimate the properties of the two stars in the binary in Sect 4.
Fortuitously, TIC 184743498 is one of the relatively rare double-line spectroscopic binaries in the Gaia data set (Katz et al. 2023;Blomme et al. 2023).The Gaia K1 and K2 velocities, also listed in Table 1, are 119.3 and 154.9 km s −1 , respectively.

Orbital period determination
We determined the orbital period of TIC 184743498 using several different approaches.
First, using all five sectors of TESS data, we performed a Box Least Squares (BLS) analysis 4 (Kovács, Zucker & Mazeh 2002).This method resulted in a period of 1.053 234 d.We also performed a BLS transform of the nearly 7,000 ASAS-SN archival flux points spanning 11 yr.The period derived from that data set is 1.053 238 d.The Gaia archives provide a period of 1.053 245 ± 0.000 014 d, based on their RV solution to the orbit.These results are summarized in Table 2.
Additionally, we generated an eclipse timing variations (ETV) 8.753 ± 0.025 Notes.curve from the TESS data.The results are shown in Fig. 3.The reference period used to compute this ETV curve was 1.053 240 d.The residuals from a linear ephemeris are nicely fitted with a quadratic term which can be described as a fractional rate of change in the orbital period of Ṗorb /P orb ≃ −(1.030 ± 0.013) × 10 −5 yr −1 .This curvature in the ETV curve could represent either orbital decay or a portion of a long-period orbit about a third body.
Given that Gaia additionally reports an 'acceleration' solution for this object, this implies a non-linear component to the proper motion (Halbwachs et al. 2023).This is consistent with orbital motion about a third body rather than orbital decay, and we adopt the former hypothesis.
The mass of the third body necessary to induce this curvature in the ETV curve can be estimated to be in the range of ∼0.7-1.3M⊙ if the outer orbit has a period in the range ∼2000-6000 d.Such an orbital period would be consistent with (i) the fact that Gaia detects an accelerated astrometric motion, (ii) the duration of the Gaia observations, and (iii) the parabolic shape of the ETV curve, as opposed to a sinusoidal shape.With this possible mass range for the third star, we have to allow for the fact that it may make a non-negligible contribution to the system light.

ESTIMATING THE STELLAR PARAMETERS
We determined the stellar parameters of the stars in TIC 184743498 with an SED fitting code that utilises the K velocities from Gaia (see Table 1) and parameter ratios from the light curve modeling code LIGHTCURVEFACTORY (see, e.g.Borkovits et al. 2019, 2020, andreferences therein).
The SED code is based on a Markov Chain Monte Carlo (MCMC; see, e.g., Ford 2005) fitting routine with only five adjustable parameters: (i) the primary mass, M1, (ii) the secondary mass M2, (iii) the mass of the tertiary star in a wide outer orbit Notes.(a) Gaia SB2C orbit (Katz et al. 2023;Blomme et al. 2023).Gaia also reports an acceleration solution (Halbwachs et al. 2023) which indicates non-linear proper motion.This is consistent with the non-linear ETV behavior we find with TESS.(b) ASAS-SN (Shappee et al. 2014;Kochanek et al. 2017).(c) See Fig. 3.
M3, (iv) the age of the system, τ , and (v) the interstellar extinction, AV .The code models the following 'input' information: (i) 21 SED points, (ii) the K1 and K2 values with uncertainties from Gaia, (iii) two ratios (T eff,1 /T eff,2 and (R1 + R2)/a) and the orbital inclination angle, i, from a fit to the orbital light curve with LIGHTCURVE-FACTORY (Borkovits et al. 2019(Borkovits et al. , 2020)), and (iv) a set of stellar evolution tracks from MIST (Dotter 2016;Choi et al. 2016;Paxton et al. 2011;Paxton et al. 2015;Paxton et al. 2019), where a solar composition is assumed.We utilise the Castelli & Kurucz (2003) model stellar atmospheres for 4 000 < T eff < 10 000 K. The distance is known with exquisite accuracy from Gaia parallax measurements.
The only constraints we have on the tertiary star are (i) from the total light of the system at each of the 21 SED points, and (ii) the fact that its mass is likely to be in the range of 0.7 to 1.3 M⊙ (based on the ETV curve, see Fig. 3 and the discussion of the figure).For the latter, we simply adopted a uniform prior on M3 over the range 0.5-1.6M⊙ to be conservative.
Our approach follows that of Kurtz et al. (2020), Rappaport et al. (2021), andRappaport et al. (2022).The MCMC code evaluates the following five parameters: M1, M2, M3, AV , and the system age via the MIST equivalent evolutionary phase (EEP) of the primary star.The use of EEPs as a preferred sampling parameterrather than directly sampling the stellar age-is described in detail in Kurtz et al. (2020).This procedure for the SED fitting has been utilised extensively in the study of numerous compact triples (Rappaport et al. 2022, Rappaport et al. 2023).The validity of the results from those SED fits has been checked in a number of cases where there are several additional constraints on the system parameters.
At each step of the MCMC routine, the procedure is as follows.The three masses and age of the binary (via the EEP) yield the stellar radii and T eff values through the stellar evolution tracks.This, of course, assumes that the two EB stars have evolved in a coeval fashion and, in particular, that they have not previously exchanged any mass.This, coupled with the model atmospheres and current value of AV , allow for a model of the measured SED points.The χ 2 value from the SED fit is then registered.Additionally, the known orbital period and the two masses provide the semi-major axis of the system.When combined with the orbital inclination angle (provided by LIGHTCURVEFACTORY), this enables us to compute K1 and K2, to compare with the measured Gaia values.The χ 2 values from the fit to the K velocities are also recorded.Finally, we check how well the model ratios of radii and T eff match those provided by LIGHTCURVEFACTORY.
The light curve that we fit with LIGHTCURVEFACTORY is shown in Figure 4, and is based on Sectors 61 and 62 of the TESS data.Instead of the usual folded, or folded, binned, and averaged data, we use a light curve that is reconstructed from the first 45 orbital harmonics5 that we fit for in the TESS data using a sim-  ilar technique to the one described later in Sect. 5.In principle, this synthetic light curve is the same as a folded, binned, and averaged light curve, except that (i) we avoid one of the pulsations that coincides with the 51st orbital harmonic and would otherwise appear in the folded light curve, and (ii) this also eliminates some of the high frequency noise that is not relevant to the shape of the light curve.The LIGHTCURVEFACTORY fit to the light curve of TIC 184743498 is shown in Figure 4 in red.The three key paramcleaning the data prior to the FT (see Sect. 5).The 45 harmonics are used here to extract the maximum high-frequency content for the light curve without running into the 51st orbital harmonic which coincides with a pulsation frequency (see Table 4).eters extracted from this fit are: T eff,1 /T eff,2 = 1.241 ± 0.010, (R1 + R2)/a = 0.500 ± 0.025, and i = 65 Figure 5 shows the results of the SED fitting.We plot the measured SED points in orange, the model SED of the primary star in red, that of the secondary in blue, and of the tertiary in green.The total flux is plotted in black.We find the following stellar parameters: M1 = 1.83 ± 0.07 M⊙, M2 = 1.37 ± 0.046 M⊙, R1 = 1.72 ± 0.06 R⊙, and R2 = 1.35 ± 0.06 R⊙.The parameters of the third star in the long outer-period orbit are less constrained, but are comparable to the properties of the secondary eclipsing binary (EB) star.The remaining fitted parameters can be found in Table 3.
We show in Fig. 6 the correlations in the MCMC posteriors for the mass vs. radius of the three stars in the TIC 184743498 system.The mass and radius of the primary are fairly well localized.By contrast, the secondary star in the EB shows a high degree of correlation between M2 and R2.This results from the fact that the secondary is lower in mass and, at the age of the system, it is still firmly on the zero-age main sequence (ZAMS).The same is true for the tertiary star (in the wide outer orbit), but the degree of correlation between M3 and R3 appears even more extreme due to the high degree of uncertainty in M3.The top panel of Fig. 6 shows the 1D posterior distributions for the stellar masses of each star.

PULSATIONS
We analyzed the pulsations of TIC 184743498 by applying a discrete Fourier transform program to the raw TESS data.We performed FTs following various stages of cleaning: We first removed the orbital harmonics from the data, and then removed the more prominent pulsation mode frequencies.As mentioned in Sect.3.1, we use Sectors 61 and 62 for pulsation analysis.
We started by subtracting the first 30 orbital harmonics from the data, the highest frequency of which is near 28.5 d −1 .This was done via a simultaneous linear least-squares fit to 30 sines and 30 cosines to represent the orbital modulations.We purposely stopped the cleaning at 30 orbital harmonics so as not to inadvertently remove any natural pulsation frequencies that happen to lie near an orbital harmonic or are, in fact, tidally excited at exactly that harmonic.The amplitudes of the orbital harmonics above this frequency are small enough so that they are not likely to be confused with stellar pulsations.Furthermore, possible residual orbital harmonics will have a phase in an echelle diagram (see below) close to 0 or 1, and will therefore be easily recognized.The middle panel of Fig. 2 shows the resultant FT in the frequency range from 0 to 80 d −1 after 30 orbital harmonics have been removed.We see a rich spectrum of pulsations within the range of 38-56 d −1 .The bottom panel of Fig. 2 shows the transform after removing the 13 most prominent frequencies (enumerated in Table 4).
To better visualize the organization of the pulsations in TIC 184743498, we generated an echelle diagram.Such a diagram plots the frequency of a pulsation on the y-axis against the echelle phase on the x-axis; this so-called "echelle phase" is the pulsation frequency modulo the orbital frequency, normalized by the orbital frequency.In the creation of an echelle diagram, we set a threshold for the minimum Fourier amplitude (≳ 7σ), such that only highly significant peaks are represented.Fig. 7 displays the echelle diagram for TIC 184743498 following the removal of 30 orbital harmonics.The amplitude threshold is 0.05 ppt (∼7 σ).We limit the echelle diagram to the relevant frequency range where the vast majority of peaks are located, and 65.2 size each dot by linearly scaling the Fourier amplitude to emphasize the more prominent peaks.We see 11 doublets (also referred to as "multiplets"), which are a vertically aligned pair of filled circles separated by twice the orbital frequency.These doublets are encircled in red ellipses and labeled by increasing frequency.These results are also tabulated in Table 4.We encircle two singlets in blue ellipses.
For each multiplet, we find the frequencies of the observed peaks (marked by ellipses in the echelle diagram) and then calculate the inferred (unseen) central frequency at which the star is actually pulsating.All the measured and inferred frequencies of the multiplets and singlets are summarized in Table 4.

RECONSTRUCTING THE PULSATION AMPLITUDES AND PHASES
The amplitudes and phases of the elements of each of the 11 pulsation doublets (see Fig. 7) carry all the information about a particular pulsation mode that can be extracted from the TESS data.However, the specific pulsation modes cannot be entirely inferred merely by inspecting the echelle diagram.In this regard, it is extremely useful to reconstruct the behavior of the pulsation amplitudes and phases as a function of orbital phase.To do this, we analytically recon- . Correlations between the mass and radius of the three stars in the TIC 184743498 system, as well as posterior distributions of the masses.The mass and radius of the secondary EB star and the tertiary are well correlated because these stars are still on the zero-age main sequence, where there is a one-to-one correspondence between mass and the other stellar parameters.Red ellipses mark the pulsation multiplets that we have identified.These are numbered according to increasing central frequency to match the notation in Table 4.
struct the amplitudes and phases of 9 of the doublets in the echelle diagram (all except for doublets 6 and 7, whose phase and amplitude behavior are not robustly characterized).We follow the formalism for reconstructing multiplets provided in Eqns.( 2)-( 6) of Jayaraman et al. (2022).
To be conservative, we force fit a quintuplet of components to each doublet, at frequencies ν0, ν0±ν orb and ν0±2ν orb .If at one or more of the 5 frequencies there is no significant FT amplitude, then the only contribution to the reconstruction will be a small amount of added noise.However, there may be a weak (albeit real) signal there which did not surpass the echelle threshold, and we cannot neglect its contribution.Note that we use the same reconstruction process for all the pulsation modes, regardless of spacing or status as a singlet or doublet.Finally, we chose the same reference eclipse time, t0, as was used in the pulsation analysis in Section 5.
Fig. 8 shows the amplitude-phase plots for the central fre-at the ELV peaks and have π phase jumps near the eclipses.The pulsation modes in Fig. 8 can all be described by dipole modes of the form Y10, namely with ℓ = 1 and m = 0, if the pulsation axis has been tidally tilted into the orbital plane and follows the tidal axis as it orbits with the binary.At low inclination angles, even lower than the 65 • for this system, the doublet nature of such a pulsation mode along the tidal axis will be preserved (see, e.g., Reed, Brondel & Kawaler 2005).This means no central peak will appear in the echelle diagram, and indeed this is the case as we look at doublets for ν1, ν2, ν3, ν4, and ν9 in Fig. 7.We hereafter refer to these as Y10,x modes.
By contrast, the reconstructions in Fig. 9 all have pulsation maxima and π phase shifts at orbital phases that are 90 • displaced from those in Fig. 8. Thus, one might be tempted to explain these pulsations as dipole modes of the form Y11 (ℓ = 1, and |m| = 1), again for a tidally tilted pulsation axis.However, for this system's relatively low orbital inclination angle, strong central frequency components must appear in the echelle diagram for Y11 modes, yet we do not see such a central component for any of the doublets in Fig. 7.This is equivalent to noting that a Y11 mode tilted along the tidal axis which, in turn, is at an inclination of 65 • , could not have zero pulsation amplitude at the eclipses -as we see in Fig. 9.
In Fig. 10, we show explicitly how the FT of a tidally tilted Y11 mode varies with orbital inclination angle.As is clear from the figure, a significant central peak begins to appear at i ≲ 80 • , and is nearly equal to the split sidelobes by i = 60 • .Thus, it would be implausible to have such a mode produce the types of doublets we see in Fig. 7, which have neither a strong central peak nor nontrivial pulsations at the eclipses.
We therefore propose that the four modes shown in Fig. 9 are actually Y10,y modes with the pulsation axis along the "y-axis", which is perpendicular to the tidal axis (x) and the orbital angular momentum axis (z).This hypothesis resolves the issue of the missing central component of the doublets, and zero pulsation amplitude at the eclipses, while being fully consistent with the pulsation amplitudes and phase shifts with respect to orbital phase.Thus, this represents the first discovery of a stellar pulsator with a pulsation axis along the y direction.
In the same spirit as the above hypothesis, we also propose that the two singlets found in this source and listed in Table 4 are Y10 modes with a pulsation axis along z.Given that their two frequencies are 52.5411 and 55.6664 d −1 , with a difference of ∼ 3 d −1 , the interval here is too narrow to host radial modes with such frequencies.This hypothesis is further explained in Section 7.
Finally, we note that the modes with frequencies ν6 and ν7, which do not appear in Figs. 8 and 9, each have markedly different amplitudes for their two doublet peaks.However, the reconstructed pulsation amplitude and phase behavior with orbital phase of ν6 very closely resembles those in Fig. 9, i.e., Y10,y modes, while ν7 would be associated with those modes in Fig. 8, i.e., Y10,x.Thus, these two doublets are likely distorted Y10 tidally tilted modes.

DENSITY OF RADIAL MODELS
In Sect.5, we reported two singlets spaced in frequency by only ∼3 d −1 .In the presence of tidally tilted pulsations, such singlets are most readily explained by radial modes, as their pulsational amplitudes and phases will not be modulated over the orbit.We investigate whether radial modes are a plausible interpretation for those two frequencies.
To this end, we computed single-star pulsational models using the latest version of the Warsaw-New Jersey stellar evolution and pulsation code (e.g., Pamyatnykh et al. 1998), for a solar chemical composition (Asplund et al. 2004) and a rotational velocity of 100 km s −1 at the ZAMS.We evaluated the frequency differences of consecutive radial overtones of modes in the domain of 50 -60 d −1 along each evolutionary track and determined the locations of models that would reproduce the observed frequency difference ν13 − ν11 at the radial mode frequency closest to ν13. Figure 11 shows the loci of those models, as well as of the stellar components of the TIC 184743498 system, in a theoretical Hertzsprung-Russell Diagram (HRD).Inspection of Fig. 11 clearly demonstrates that (i) only the primary star (star 1) is expected to exhibit δ Scuti type pulsations, and (ii) frequencies ν11 and ν13 cannot both be radial modes, as all stellar models that produce consecutive radial modes that are closely spaced in frequency are too evolved.Models producing ν11 and ν13 as non-consecutive radial modes would be even more evolved.
Additionally, we look at the expected pulsation spectrum of dipole modes of a single star with the parameters we derived for the primary star in the TIC 184743498 system.The corresponding pulsation frequencies were computed analogously to those of the radial modes with the Warsaw-New Jersey code, for a model with log T eff = 3.9206, and log L = 1.1302 (cf.Table 3).Nonradial mode rotational frequency splittings were computed to first order in stellar rotation frequency.The schematic frequency spectrum so derived is shown in Fig. 12.
Evidently, the expected spectrum of axisymmetric dipole modes (Y10) is not dense enough to explain all the observed frequencies.As the underlying model is fairly unevolved, no avoided crossings (mixed modes) are present in the theoretical frequency spectrum; the modes shown here are of radial overtone 5 to 7.However, by including sectoral dipole modes (Y1±1), the density of the observed apparent dipole frequency spectrum can be explained.
The above arguments lead us to the conclusion that at least one of the ν11 and ν13 singlets, and possibly both, must indeed be a Y10 mode with a pulsation axis along z, i.e., Y10,z modes.

TIDAL MODE-COUPLING MODELS
As suggested in Section 5, the pulsations of TIC 184743498 can be best explained by three different pulsation axes.In particular, it appears that there are five dipole Y10 modes with pulsation axis aligned with the tidal (or x axis), while four are dipole Y10 modes with a pulsation axis aligned along the y axis (i.e., in the orbital plane, but orthogonal to the x axis).Here we justify this hypothesis with a straightforward physical argument.
Rather than assuming that tidal distortion aligns pulsations with the tidal axis, we show that tidal coupling between modes naturally produces modes with three different pulsation axes.Consider an ℓ = 1 triplet, now using the rotation/orbital axis as a reference axis.We formally define the x−axis in the direction of the companion, the z−axis in the direction of the spin/orbital angular momentum, with the y−axis lying in the orbital plane.The tidal distortion is dominated by the ℓ = m = 2 component of the tidal potential, which couples the ℓ = m = 1 mode with the ℓ = 1, m = −1 mode.Writing out the coupled eigensystem for this mode triplet The reconstructed pulsation amplitude and phase variations as a function of orbital phase for the five doublets that we have classified as Y 10,x pulsation modes with a pulsation axis along the tidal axis (ν 1 through ν 4 and ν 9 ); the frequencies are listed in Table 4.The reconstructed pulsations have been vertically offset from each other by 0.01 for clarity, and each pulsation amplitude has been doubled from its actual value to make it more visible.The phase reconstructions are offset vertically from one another by 2 cycles.We label each reconstructed mode with its frequency numbering.The black curve superposed on the plots is the reconstructed orbital modulation.See text for details of the reconstructions.The reconstructed pulsation amplitude and phase variations as a function of orbital phase for the four doublets that we have classified as Y 10,y pulsation modes with a pulsation axis along the y direction (perpendicular to the tidal and angular momentum axes).The frequencies ν 5 , ν 8 , ν 10 , and ν 12 are listed in Table 4.The reconstructed pulsations have been vertically offset from each other by 0.01 for clarity, and each pulsation amplitude has been doubled from its actual value to make it more visible.The phase reconstructions are offset vertically from one another by 2 cycles.We label each reconstructed mode with its frequency numbering.The black curve superposed on the plots is the reconstructed orbital modulation.See text for details of the reconstructions.(i.e., Eqn. 8 of Fuller et al. 2020), we find Here, ω is the frequency of a mode of the coupled system, ω 2 α is the unperturbed frequency of the mode triplet, and δV1−1 and δT1−1 Since the m = 0 mode remains uncoupled in this simple scenario, its eigenfrequency is slightly perturbed and its eigenfunction remains unchanged.It is an m = 0 mode about the z−axis, i.e., the spin/orbital axis.The m = ±1 modes are coupled, however, and the eigensystem of Eqn. 1 for those modes reduces to where δω 2 tide = δV1−1 − ω 2 α δT1−1 . (5) Here, we have assumed small perturbations such that δω 2 1 , δω 2 −1 , δV1−1 ≪ ω 2 α , and δT1−1 ≪ 1, and we have dropped second-order terms in these small quantities.
However, in the limit of strong tidal coupling such that δω 2 tide ≫ δω 2 1 − δω 2 −1 , solving the eigensystem yields the new mode eigenfrequencies with corresponding eigenvectors Hence, the new mode eigenfunctions are equal superpositions of Y1,1 and Y1,−1, with angular flux perturbation patterns Some algebra shows that the spatial/time dependence of these two modes are δF+ ∝ sin θ cos ϕ cos(ω+t) and Note also that the m = 0 mode has δF0 ∝ cos θ cos(ωαt) The extra subscripts on the Y10 spherical harmonics (x, y, z) refer to the axis with respect to which the spherical harmonic is defined.Note that neither of the modes in Eqn. 9 or Eqn. 10 propagates around the equator like uncoupled m = ±1 modes.Instead, both modes are standing modes, aligned with the x− and y−axes.The m = 0 mode is a standing mode aligned with the z−axis.Thus, we find that strong tidal coupling naturally transforms an ℓ = 1 dipole triplet of states to a set of Y10 modes around three orthogonal axes.Hence, tidally coupled ℓ = 1 modes can behave identically to the tri-axial pulsation behavior necessary to explain the observations of TIC 184743498.Tidal coupling will have a similar effect on ℓ = 2 modes, coupling the Y21 and Y2−1 modes into a perturbed doublet, and the Y22, Y20, and Y2−2 modes into a triplet.This will produce more complicated spatial patterns than those discussed above, with different amplitude/phase modulation over the orbit.We plan to examine these signatures in future work and simultaneously search for these signatures in other stars in the TESS data.
As mentioned above, tidal coupling also couples modes of different ℓ, which is not accounted for above.For weak tidal distortion, the ℓ = 2 component of the tidal distortion dominates, such that coupling between modes differing by ∆ℓ = 2 and ∆ℓ = 0 is strongest.For stronger tidal distortion, the ℓ = 3 component of the tidal distortion becomes more important, enabling coupling between modes differing by ∆ℓ = 3 and ∆ℓ = 1.In this case, the star and mode eigenfunctions will be asymmetric across the x − yplane, such that pulsations can be strongly trapped on either side of the star (i.e., the side facing toward or away from the companion).This leads to the observed "tidal trapping" or "single-sided pulsator" phenomenon (e.g., Handler et al. 2020;Kurtz et al. 2020).
In general, tidal coupling could lead to complex behavior of tidally perturbed pulsations, and the effects of strong tidal distortion, centrifugal, and Coriolis forces should all be taken into account.However, the simple case shown above demonstrates a straightforward mechanism through which tidal coupling could induce pulsations to naturally align with the x, y, and z−axes of the star.More thorough calculations incorporating the effects listed above, and using larger networks of coupled modes, will be needed for a full understanding of pulsations of tidally distorted stars.

SUMMARY, CONCLUSIONS, AND FUTURE DIRECTIONS
In this work, we report the discovery and analysis of a tight eclipsing binary with at least nine tidally tilted pulsation modes, TIC 184743498.
Evidence from our own ETV curve, as well as from Gaia's discovery of an astrometric acceleration solution, both indicate that there is a third star in the system in a wide orbit of thousands of days.We have analysed the available archival data for this system to investigate the stellar parameters and the evolutionary state of the system using an SED fitting code.This code incorporates measured RVs from Gaia, stellar evolution models, and other parameters extracted from the light curve with the modeling code LIGHTCURVEFACTORY.We find that the pulsating primary star has M1 ≃ 1.83 M⊙, R1 ≃ 1.72 R⊙, and T eff,1 ≃ 8500 K, while the secondary has M2 ≃ 1.37 M⊙, R2 ≃ 1.35 R⊙, and T eff,2 ≃ 6870 K.The primary star has evolved somewhat away from the ZAMS at an age of 460 Myr, while the secondary star is still on the ZAMS.The uncertainties on the properties of the inferred third star are large, but it is likely not too different from the secondary star in the inner binary, and apparently contributes ∼10% of the system light.
We then analyzed the pulsations of TIC 184743498 in detail.The system exhibits nine pulsation modes, 5 of which have pulsation amplitude maxima at the binary eclipses and phase shifts of π at the ELV peaks.We conclude that the first 5 pulsations can be explained by dipole modes Y10 about a pulsation axis which has been tilted into the orbital plane and which lies along the tidal axis.We call these 'Y10,x' modes.
Four other pulsations exhibit amplitude maxima at the ELV peaks and π phase shifts at the eclipses.We have shown that while the latter modes may appear to be Y11 modes based on their amplitudes and phases, this explanation is untenable because at the relatively low orbital inclination angle of this system (65 • ), such doublets would have an unmistakable central peak in their echelle diagram.Such a central peak is a visual description that the mode amplitude is nonzero at t0 (i.e., at the eclipses) for these modes; however, the observed mode amplitudes do, in fact, vanish at the eclipses.The lack of such a central peak shows that these modes are actually Y10 about an axis we define as y, which is orthogonal to the tidal axis and the angular momentum axis.We term these 'Y10,y' modes.
We further found that one or both of two remaining singlet pulsation modes are also Y10 modes, but about the angular momentum axis z, and we investigated this hypothesis with single-star pulsation models.We found that the expected pulsation spectrum of radial modes with a star of the parameters derived from the SED fit is significantly less dense than the observed pulsation frequencies.Thus, the singlet pulsation modes cannot both be radial pulsations, but rather at least one of them must be a Y10,z mode, i.e., about the z-axis.
To study this hypothesis of three different pulsation axes, we explored a natural explanation in which strong tidal coupling can convert an otherwise uncoupled triplet of ℓ = 1 dipole states to all Y10 modes, but about three orthogonal axes.These axes are (i) the tidal axis, (ii) the orbital angular momentum (i.e., stellar rotation) axis, and (iii) the direction in the orbital plane perpendicular to the tidal axis.This is exactly what we are seeing in the total set of 9 dipole modes and two singlet modes observed in TIC 184743498.Additionally, we show that the transformed original Y1±1 modes, which have circulating pulsation modulations, are reduced to standing waves.
Every newly discovered tidally tilted pulsator presents distinctly unique behavior in its pulsations, and challenges previously understood ideas about tidally tilted pulsations.TIC 184743498 contributes to this diversity as the first identified tri-axial pulsator.
More generally, the discovery reported here of a tri-axial set of pulsation axes has a wider application to the broad subject of stellar pulsations.We have seen that simple modes (such as doubletdipoles) are easier to understand in the context of tidal tilting than complex multi-element modes; the latter of which are, in contrast, far easier to spot in echelle diagrams.Finding the numerous doublets in TIC 184743498, and their interpretation, have allowed us to identify and interpret the modes of nearly all the prominent pulsations in this star.
With this in mind, we hope to search for more TTP candidate stars with simpler modes such as doublets rather than those with complex multi-element modes, such as HD 265435 (Jayaraman et al. 2022).Future observations of TTPs in different photometric bands may reveal new insights into mode cavities/propagation at different stellar layers, which, in turn, can provide an exquisite new view into stellar interiors.We also highlight the fact that toy models such as the one presented in this work should allow us to develop a preliminary understanding of previously-unexplainable phenomena, such as complicated mode couplings in stellar interiors.
In this work we have also shown that, in systems with tidally tilted modes, part of the pulsation geometry can be constrained, in some cases rather tightly, by the determination of the orbital and system parameters, especially the orbital inclination angle.It was just this tilt of the orbit that allowed us to realize that we were dealing with tri-axial pulsations in TIC 184743498, otherwise, the pulsational multiplet structures would be different.This orbital constraint can, and should, be applied to other TTP systems.
This work also motivates us to reexamine other previously identified TTP candidates to see if any of their multiplets can be fit into this paradigm of tidally induced doublets.In this regard, the simple model for tidal tilting presented here needs expansion to quadrupole modes.Once the signatures of tidally tilted quadrupole modes, under the paradigm developed in this work, are understood, we will be able to search for these specific features in other pulsating stars in binaries.Given a small set of stars with precisely measured stellar parameters, and nearly complete sets of identified pulsation modes, detailed stellar modeling of these modes will be warranted.
Finally, we note that the 200-second cadence data that will continue to come from TESS's all-sky survey will be invaluable in helping us determine the relationship between the stellar pulsation axes and the strength of the tidal perturbation.These data will also allow us to find all manner of new types of pulsations, and perhaps answer some of the lingering questions in the field.

Figure 1 .
Figure 1.A portion of the Sector 62 TIC 184743498 light curve showing the eclipses, ellipsoidal light variations, and prominent pulsations.The TESS cadence for this sector was 120 s.

Figure 2 .
Figure 2. The Fourier amplitude spectrum of TIC 184743498.The panels highlight the amplitude spectrum after various stages of cleaning.The top panel shows the spectrum of the raw TESS light curve; the middle panel represents the spectrum following the removal of 30 orbital harmonics, while the bottom results from a further cleaning of the 13 pulsation modes discussed in the paper and listed in Table4.Note the changes in amplitude scale between the upper and lower panels.

Figure 3 .
Figure 3. ETV curve based on 107 primary eclipse times from TESS Sectors 8, 34, 35, 61, and 62.The reference period for this analysis was 1.053 240 0 d.The non-linear behavior is significant at the 80 σ level.

Figure 4 .
Figure 4. Model fit to the TESS sector 61 and 62 data.The black curve is the light curve reconstructed from the first 45 orbital harmonics.The red curve is a fitted model using LIGHTCURVEFACTORY.See text for details and references.

Figure 5 .
Figure5.21 SED points for TIC 184743498 (orange points) are fitted to stellar models for the three stellar components in this system.The contribution of the primary star is shown as a red curve, secondary is shown in blue, and tertiary is shown in green.The total model flux is given by the black curve.See text for the constraints that were used in the fit.

Figure 7 .
Figure 7. Echelle diagram following a cleaning of the data that removed the first 30 orbital harmonics.This diagram shows pulsation peaks as a function of their frequency and echelle phase, i.e., the normalized pulsation frequency modulo the orbital frequency.Points are sized linearly by their Fourier amplitude.Red ellipses mark the pulsation multiplets that we have identified.These are numbered according to increasing central frequency to match the notation in Table4.

Figure 8 .
Figure8.Left and right panels: The reconstructed pulsation amplitude and phase variations as a function of orbital phase for the five doublets that we have classified as Y 10,x pulsation modes with a pulsation axis along the tidal axis (ν 1 through ν 4 and ν 9 ); the frequencies are listed in Table4.The reconstructed pulsations have been vertically offset from each other by 0.01 for clarity, and each pulsation amplitude has been doubled from its actual value to make it more visible.The phase reconstructions are offset vertically from one another by 2 cycles.We label each reconstructed mode with its frequency numbering.The black curve superposed on the plots is the reconstructed orbital modulation.See text for details of the reconstructions.

Figure 9 .
Figure9.Left and right panels: The reconstructed pulsation amplitude and phase variations as a function of orbital phase for the four doublets that we have classified as Y 10,y pulsation modes with a pulsation axis along the y direction (perpendicular to the tidal and angular momentum axes).The frequencies ν 5 , ν 8 , ν 10 , and ν 12 are listed in Table4.The reconstructed pulsations have been vertically offset from each other by 0.01 for clarity, and each pulsation amplitude has been doubled from its actual value to make it more visible.The phase reconstructions are offset vertically from one another by 2 cycles.We label each reconstructed mode with its frequency numbering.The black curve superposed on the plots is the reconstructed orbital modulation.See text for details of the reconstructions.

Figure 10 .
Figure 10.Fourier transforms of simulated pulsation for a mode Y 11 with a pulsation axis along the tidal axis of the binary (x-axis).The orbital inclination angle for each simulation is labeled in color.The curves for different inclinations are separated both vertically and horizontally for clarity.A central peak starts to appear at i ≲ 80 • , and is nearly equal to the split sidelobes by i = 60 • .

Figure 11 .
Figure 11.Theoretical HRD with the locations of the three components of TIC 184743498 (stars 1, 2, and 3) indicated.The evolutionary tracks are the same as were used in the SED fit shown in Fig. 5.The dashed black lines are the observed boundaries of the δ Scuti star instability strip (Murphy et al. 2019).The green line connects models that best represent the observed pulsation frequencies ν 11 and ν 13 with radial modes.

Figure 12 .
Figure 12.Schematic pulsation frequency spectrum of TIC 184743498.The black lines denote the observed centroid mode frequencies, whereas the blue lines indicate the theoretical model radial-mode frequencies.Red lines represent theoretical dipole mode frequencies, with those of 'amplitude' 0.6 corresponding to m = 0 and those of 'amplitude' 0.4 to |m| = 1.

Table 1 .
Properties of the TIC 184743498 System

Table 2 .
Orbital Period of TIC 184743498

Table 3 .
Derived Parameters for the TIC 184743498 System