The Transiting Exoplanet Survey Satellite TESS has enabled the discovery of numerous tidally tilted pulsators TTPs which are pulsating stars in close binaries where the presence of a tidal bulge has the effect of tilting the primary stars pulsation axes into the orbital plane. Recently the modeling framework developed to analyze TTPs has been applied to the emerging class of triaxial pulsators which exhibit nonradial pulsations about three perpendicular axes. In this work we report on the identification of the secondever discovered triaxial pulsator with 16 robustly detected pulsation multiplets of which 14 are dipole doublets separated by 2 SUBorbSUB. We jointly fit the spectral energy distribution and TESS light curve of the star and find that the primary is slightly evolved off the zeroage main sequence while the less massive secondary still lies on the zeroage main sequence. Of the 14 doublets we associate eight with Y SUB10x SUB modes and six with novel Y SUB10y SUB modes. We exclude the existence of Y SUB11x SUB modes in this star and show that the observed pulsation modes must be Y SUB10y SUB. We also present a toy model for the triaxial pulsation framework in the context of this star. The techniques presented here can be utilized to rapidly analyze and confirm future triaxial pulsator candidates.
['Asteroseismology', 'Tidal distortion', 'Detached binary stars', '73', '1697', '375', 'Astrophysics - Solar and Stellar Astrophysics']
TIC 435850195 The Second Triaxial Tidally Tilted Pulsator
{'TIC 435850195: The Second Tri-Axial, Tidally Tilted Pulsator': 'Rahul Jayaraman , 1 Saul A. Rappaport , 1 Brian Powell , 2 Gerald Handler , 3 Mark Omohundro, 4 Robert Gagliano, 5 Veselin Kostov , 2, 6 Jim Fuller , 7 Donald W. Kurtz , 8, 9 Valencia Zhang, 10 and 1  \nGeorge Ricker  \n1 Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA  \n2 NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  \n3  \nNicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, PL-00-716 Warszawa, Poland 4 Citizen Scientist, c/o Zooniverse, Dept., of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK 5 Amateur Astronomer, Glendale, AZ 85308, USA  \n6 SETI Institute, 189 Bernardo Avenue, Suite 200, Mountain View, CA 94043, USA  \n7 TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA  \n8 Centre for Space Research, North-West University, Mahikeng 2745, South Africa  \n9 Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK 10 Phillips Academy, Andover, MA 01810, USA  \n(Dated: September 9, 2024)  \nSubmitted to ApJ', 'ABSTRACT': "The Transiting Exoplanet Survey Satellite (TESS) has enabled the discovery of numerous tidally tilted pulsators (TTPs), which are pulsating stars in close binaries where the presence of a tidal bulge has the effect of tilting the primary star's pulsation axes into the orbital plane. Recently, the modeling framework developed to analyze TTPs has been applied to the emerging class of tri-axial pulsators, which exhibit nonradial pulsations about three perpendicular axes. In this work, we report on the identification of the second-ever discovered tri-axial pulsator, with sixteen robustly-detected pulsation multiplets, of which fourteen are dipole doublets separated by 2 ν orb . We jointly fit the spectral energy distribution (SED) and TESS light curve of the star, and find that the primary is slightly evolved off the zero-age main sequence, while the less massive secondary still lies on the zero-age main sequence. Of the fourteen doublets, we associate eight with Y 10 x modes and six with novel Y 10 y modes. We exclude the existence of Y 11 x modes in this star and show that the observed pulsation modes must be Y 10 y . We also present a toy model for the tri-axial pulsation framework in the context of this star. The techniques presented here can be utilized to rapidly analyze and confirm future tri-axial pulsator candidates.", '1. INTRODUCTION': "The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) has revealed diverse classes of pulsating stars, including the novel class of 'tidally tilted pulsators' (TTPs), wherein pulsating stars in tight binaries (with periods typically less than ∼ 2 d) have their pulsation axes tilted into the orbital plane by the tidal bulge induced by their companions. Specifically, in these systems, the pulsation axis is aligned with the tidal bulge,  \nrather than with the spin axis of the star. To date, there have been a handful of tidally tilted pulsators discovered (Handler et al. 2020; Kurtz et al. 2020; Rappaport et al. 2021; Jayaraman et al. 2022a; Zhang et al. 2024; Jennings et al. 2024b), and several other candidates reported in the literature (see, e.g., Kahraman Ali¸cavu¸s et al. 2024).  \nA key property of tidally tilted systems is that the observer is able to study the star through a wide range of latitudinal angles (from 0 · to 360 · ) with respect to the pulsation axis. Tidally tilted modes exhibit amplitude and phase modulations around the orbit, thereby leading to splitting of the modes, seen in periodograms of their light curves. Unlike typical pulsators, where the  \nobserver's viewing angle remains constant with respect to the pulsation axis, tidally tilted pulsations yield a unique perspective on pulsating stars as they orbit their respective companions.  \n'Tidally perturbed' pulsators, which also exhibit oscillation modes whose amplitudes are modulated over the course of the orbit, could also be interpreted as less extreme manifestations of the tidal tilting phenomenon (see, e.g., Southworth et al. 2020; Van Reeth et al. 2022, 2023; Johnston et al. 2023; Jennings et al. 2024a). However, there may also be other factors at play in these stars causing the observed changes in the pulsation modes, such as asynchronous rotation (see, e.g., the Introduction of Bowman et al. 2019 and references therein) or a tidal amplification mechanism, described in Fuller et al. (2020).  \nSignificant modeling of the first three tidally tilted pulsators that were discovered was performed by Fuller et al. (2020); a unique attribute of the tidal tilting phenomenon is that the pulsation modes are able to be identified with specific ℓ and | m | values, providing valuable insights into both the star's evolutionary state and its interior structure (as was shown in Jayaraman et al. 2022a). Initially, the framework for understanding TTPs was based upon the oblique pulsator model, developed by Kurtz (1982) for the rapidly oscillating Ap (roAp) stars; in this scheme, the gravitational field of the companion star in TTPs plays an analogous role to the magnetic field intrinsic to roAp stars. However, it remains an open question why only certain tight binaries experience tidal tilting or, for that matter, why certain modes in a given star can be tidally tilted, while others in the same star are not.  \nRecently, Zhang et al. (2024) (hereafter Z24) reported the unique discovery of TIC 184743498, the first 'tri-axial pulsator' (TAP). This star was found to pulsate along three different axes, and Section 8 of Z24 presented a perturbative model for how modes coupled by the tidal bulge could naturally produce pulsations along three different axes. This work presents the second TAP discovered with TESS. Section 2 presents the details of the TESS data used in the analysis, and our initial processing of the light curve. Section 3 presents our joint fit to the system's light curve and spectral energy distribution (SED), and enumerates the system parameters resulting from that analysis. Section 4 presents an asteroseismic analysis of the system, and Section 5 presents a toy model of the tidal tilting framework for this system, while exploring (and refuting) alternative explanations for this star's behavior.", '2. OBSERVATIONS AND DATA PROCESSING': "Figure 1. A ∼ 3.5 d portion of the TESS light curve of TIC 435850195 from the Quick-Look Pipeline (Huang et al. 2020a,b), sampled at 200 s cadence. The amplitudes and phases of most of the pulsation modes are seen to vary systematically with the orbital phase. There are also prominent ellipsoidal light variations (discussed further in Section 3.1.2) and prominent primary eclipses.  \n<!-- image -->  \nA team of citizen scientists, including MO and RG, have been conducting a visual survey of light curves from the TESS full-frame images (FFIs); more information about this effort can be found in Kristiansen et al. (2022). One major goal of this survey is to find unique stars in eclipsing binary (EB) systems. To this end, Powell et al. (2022) generated a set of millions of EB light curves from the TESS full-frame images using the eleanor pipeline (Feinstein et al. 2019) and a machine learning approach, detailed in Section 2 of Powell et al. (2021). TIC435850195 was initially identified as part of this effort, during a review of EB light curves identified by the neural network from TESS Sector 56. This star was observed for the first time in Sector 56 at 200 s cadence. Further details about the system are given in Table 1.  \nFollow-up analyses utilized the 200 s-cadence light curve from the MIT Quick-Look Pipeline (QLP; Huang et al. 2020a,b). The QLP light curve was downloaded using lightkurve (Lightkurve Collaboration et al. 2018), and we used the orbital period calculated by the Gaia mission (Gaia Collaboration et al. 2016, 2023), reported in Table 1 (Siopis et al., in prep), to analyze the eclipses. Orbital period estimates from both the TESS data and archival data from the All-Sky Automated Survey for Supernovae (Kochanek et al. 2017) were consistent with the Gaia value. Using the TESS data, we also calculated a value for the ephemeris t 0 , using the time of the first primary eclipse. We find that t 0 = BJD2459826.36448 (where BJD is 'Barycentric Julian Date' in the TDB system, as defined in Eastman et al. 2010).  \nAfter downloading the normalized light curve, we fit for 50 harmonics of the orbital period, along with a  \nFigure 2. The Discrete Fourier Transform (DFT; Kurtz 1985) of TIC 435850195, with sequential subtractions of the eclipsing binary signal in the raw light curve (visible in the top panel) and the δ Scuti pulsations (visible in the middle panel). Each set of periodic signals was fit for using leastsquares and then subtracted from the light curve; the bottom panel shows the DFT of the residual light curve, after both the EB signal and pulsations were fitted and removed. Note the differing scales between the top panel and the lower two.  \n<!-- image -->  \nconstant offset, to the light curve. This procedure was used to reconstruct the eclipsing binary light curve for subsequent analysis with our light curve fitting code. We used the residuals from this fit as the 'pulsational' light curve for our asteroseismic analysis in Section 4. With this 'pure' pulsational light curve, we then fit for the most significant frequencies, including their amplitudes and phases. High-amplitude singlet frequencies, and those frequencies that formed part of a multiplet, are enumerated, along with their best-fit amplitude and phase, in Table 4. The Discrete Fourier Transforms Kurtz 1985 of each light curve (raw QLP light curve, 'pure' pulsational light curve, and the residuals after all the frequencies are subtracted), are shown in Figure 2.  \nTable 1. Information about the TIC 435850195 system. The reference numbers in the third column correspond to: (1) Gaia Collaboration et al. (2016, 2023); (2) Stassun et al. (2019); (3) Katz et al. (2023); (4) Siopis et al. (in prep).  \nThe spectral energy distribution was downloaded from the Vizier SED viewer 1 (Ochsenbein et al. 2000). We used the photometric measurements of this system available from the following sources: Gaia (Gaia Collaboration et al. 2016, 2023), Pan-STARRS1 (Chambers et al. 2016), the Wide-field Infrared Survey Explorer (Wright et al. 2010), the 2-Micron All-Sky Survey (Cutri et al. 2003; Skrutskie et al. 2006), the Tycho-2 catalogue from the Hipparcos survey (Høg et al. 2000), the Sloan Digital Sky Survey (Ahumada et al. 2020), and the Galaxy Evolution Explorer ultraviolet source catalog (Bianchi et al. 2017).", '3. ESTIMATING BINARY SYSTEM PARAMETERS': "We simultaneously fit the orbital light curve and the composite SED of TIC 435850195 using a custom Markov Chain Monte Carlo (MCMC) fitting code. There are 14 (or 16) free parameters in the fit: the two masses ( M 1 , M 2 ), the system age, the orbital inclination angle i , 8 physically motivated Fourier components of the light curve (and an additional constant offset), the distance to the system, and the line-of-sight extinction A G . The aforementioned eight Fourier components correspond to the first four orbital harmonics, and represent a combination of the ellipsoidal light variations (ELVs), illumination effects, and any co-rotating starspots; physical explanations for these components are given in Table 2.  \nTo characterize this system, we make use of the MIST stellar evolution tracks (Dotter 2016; Choi et al. 2016), which yield the stellar radii ( R 1 , R 2 ), effective temperatures ( T eff , 1 ; T eff , 2 ), and luminosities ( L 1 , L 2 ) as a func-  \ntion of the stellar masses and system age. As a result, at each step of the MCMC, we are able to utilize the masses to calculate the orbital separation a (from the period and Kepler's third law), allowing us to compute the eclipse geometry. We can also use the stellar radii and effective temperatures to model the composite SED, using the model atmospheres of Castelli & Kurucz (2003). As part of our analysis, we assume that the two stars have evolved in a coeval manner since their formationin particular, we assume that there have been no prior episodes of mass transfer between them.  \nThe final possible free parameters are the distance to the source and the line-of-sight extinction A G . We first ran an MCMC with both the distance and the extinction as free parameters, to ensure agreement with the best-fit values from Gaia. Then, we fixed the distance to 521 pc, and the extinction to 0.42, and re-ran the fit. For each value in the SED, we corrected for extinction using the Cardelli et al. (1989) extinction law and this A G value.  \nIn the following subsections, we describe in more detail the simplified orbital light curve model and the SED fit. The MCMC was run with ten parallel walkers, for 175000 steps. We discarded the first 20 000 samples in each chain as burn-in.", '3.1. Light Curve Fitting': "We performed a fit to the folded and binned Fourierreconstructed light curve. This light curve was generated using the coefficients from the fit to 50 orbital harmonics described in Section 2, and was binned to 2 min. This approach was taken in order to minimize any leakage of signals from the pulsation frequencies into the 'pure' eclipsing binary light curve. Our simplified light curve fitting code incorporates three parts: (i) eclipses from two spherical stars 2 ; (ii) limb darkening, for which we utilize the approximation scheme described in Appendix A; and (iii) a series of four sine and four cosine terms to capture out-of-eclipse variability, including ellipsoidal light variations, the illumination effect, and starspots that are co-rotating with the orbit. The bolometric correction for the TESS band, which is a function of the stellar effective temperature, was approximated by assuming a blackbody spectrum and then calculating the ratio of the integrated flux between 6 000 ˚ Aand 10000 ˚ Ato the total integrated flux across all wavelengths. Bolometric  \nFigure 3. Our best-fit model to the binned, folded, Fourierreconstructed light curve of TIC 435850195, along with residuals. The residual structure around the times of eclipse could be partially explained by our use of fixed limb-darkening coefficients, rather than letting them be free parameters.  \n<!-- image -->  \ncorrection values were calculated for blackbodies having temperatures between 4 000 K to 10 000 K.", '3.1.1. Eclipse Modeling': "At each link in the MCMC chain, we use the stellar radii and luminosities for the eclipse modeling. We set the baseline flux as the sum of the two luminosities, multiplied by their respective bolometric corrections. We then calculated the sky-projected separation between the two stars' centers for the current values of the orbital inclination angle i and semi-major axis a , which then allowed us to determine the overlapping area between the two stars as a function of orbital phase.  \nDuring the portions of overlap, we utilized a onedimensional integral to calculate the limb darkening (LD) flux variations for the primary star. We used a quadratic LD model and the coefficients from Claret (2017). We numerically integrated the limb-darkening over the overlapping area, in a scheme that is a simplification of the standard approach from Mandel & Agol (2002). The relatively simple geometry of this system allows for a large number of the cases (III-XI) in Mandel & Agol to be omitted for the case of the primary eclipse. Further details of our approach to calculating the limb darkening can be found in Appendix A. For the secondary star, whose limb darkening has a significantly smaller effect on the light curve shape, we calculated an 'average' limb-darkening coefficient for the overlap region and then applied this across that entire area; this approach provided a reasonable fit (see residuals in Figure 3).  \nFor the primary, we used the limb-darkening coefficients for a star with T eff = 7500K, log g = 3 . 5, and  \nTable 2. Best-fit amplitudes to our semi-empirical model of the out-of-eclipse variability. Two of these parameters ( a 2 and b 1 ) correspond to physical attributes of the system, and are further discussed in Section 3.1.2. Other parameters account for co-rotating spots on the stellar surface.  \nmicroturbulent velocity ξ =2kms -1 . For the secondary, we used the coefficients for a star with T eff = 4250K, log g = 4 . 0, and microturbulent velocity ξ = 0kms -1 . We recognize that there exists a significant debate in the exoplanet community over whether limb-darkening coefficients should be fixed or free parameters when fitting transit light curves (see, e.g., Csizmadia et al. 2013; Espinoza & Jord'an 2015). Tests with LD coefficients as free parameters did not yield a significant improvement in the residuals, so we fixed them to the Claret values.  \nTwo orbital cycles of the Fourier-reconstructed light curve, with our best-fitting model overlaid (including the out-of-eclipse fit described in Section 3.1.2), are shown in Figure 3. Our model fits the light curve rather well, with a maximum residual of approximately 1 ppt.", '3.1.2. Out-of-Eclipse Variability': 'For the eclipse modeling, we assumed a flat out-ofeclipse baseline that was normalized to 1. However, in reality, light curves of tight binaries exhibit significant out-of-eclipse variability, for the reasons described earlier. As a result, we added in a semi-empirical model with nine free parameters that has the functional form, which is similar to Equation 8 from Carter et al. (2011):  \nc 0 + 4 ∑ n =1 a n cos nω orb t + b n sin nω orb t. (1)  \nThe values of these parameters and their physical significances are enumerated in Table 2. With these best-fit parameter values, we can attempt to constrain certain physical attributes of the system.  \nEllipsoidal Light Variations -Stars in tight binaries are distorted into ellipsoids by the tidal forces from their companion (a more detailed discussion can be found in  \nSection IV.2 of Kopal 1959). This affects the amount of flux that we observe, due to the fact that both the crosssectional area from which the flux appears to be emitted and the gravity darkening vary throughout the orbit with a frequency of 2 ν orb . These contribute to the dominant term in the ELV (cos 2 ωt ); terms with frequencies ν orb and 3 ν orb also contribute to the ELV, but to a much lesser degree. Equations 14-16 from Carter et al. (2011) relate the amplitude of each cosine term ( a 1 , a 2 , and a 3 in our notation) to physical parameters.  \nThe ELVs are analytically related to the physical properties of the system through the linear limb-darkening coefficient u , and the gravity-darkening exponent y (von Zeipel 1924). For TIC 435850195, we use the u and y values derived for the TESS band by Claret (2017); these are 0 . 4114 and 0 . 1245, respectively. The equation relating a 2 and the physical attributes of the system is  \na 2 ≃ -Z 1 (2) q ( R 1 a ) 3 sin 2 i. (2)  \nWe present below a simplified version of the expression for Z 1 given in Morris (1985). Our expression does not include the k 1 term in the original equation that accounts for precession induced by a third body.  \nZ 1 (2) = 45 + 3 u 20(3 -u ) (1 + y ) = 1 . 004 (3)  \nThe dominant source of uncertainty in this expression arises from the uncertainty in the coefficients u and y . We substitute our value of Z 1 (2) into Equation 2, along with the relevant best-fit parameters from Table 3. Using our calculated value of a = 6 . 96 R ⊙ , we find that a 2 ≃ 9813 ppm. This is ∼ 70% of the best-fit value for a 2 in Table 2, suggesting that the ELV contribution is the dominant term for the component of the out-of-eclipse variability with a frequency twice that of the orbit.  \nDoppler Boosting -Of particular interest is the b 1 sin ωt term, which corresponds to the Doppler boosting effect in tight binaries first detected by Maxted et al. (2000); a theoretical treatment of this effect is given in Loeb & Gaudi (2003). Doppler boosting is caused by the motion of the primary inducing three key effects, all of which are comparable in magnitude: An increased rate of photon arrivals, an increase in the net energy of the emitted photons, and slight relativistic beaming due to the motion of the star. Using synphot (STScI Development Team 2018), we calculated the Doppler boosting coefficient α TESS to be 2.89 for a 7800 K star (assuming a blackbody spectrum); for a Vega-type star, this coefficient increases to 2.92.  \nThe Doppler boosting amplitude is given by  \nA DB = α K c , (4)  \nwhere α is the coefficient calculated previously. If we were to use the Gaia value for K 1 (46 km s -1 ), we expect the coefficient b 1 of the sin ω orb t term to be +443 ppm. By comparison, Table 2 shows a value for b 1 = -83 ppm. This discrepancy can be explained by invoking co-rotating starspots, as these can strongly affect the amplitude of low-frequency orbital Fourier terms. We can gauge how large these might be by looking at terms a 4 , b 2 , b 3 , and b 4 , where there should be minimal to no contribution from physical effects such as ellipsoidal light variations, illumination effects, and Doppler boosting. At least two of these amplitudes ( a 4 and b 2 ) are large enough to interfere with the Doppler boosting term ( b 1 ), if they occur at the orbital frequency. Thus, we conclude that the Doppler boosting term is undetectable here.  \nAn interesting question is whether spots in this system can account for amplitudes of ≳ 1000 ppm at the lower orbital harmonics. The secondary, which is cool enough to have large spots, contributes ∼ 1% of the system light. To produce signals as large as 1000 ppm requires spots of 10% amplitude on the secondary, which is plausible.', '3.2. SED Fitting': 'Simultaneously with the light curve fit, we also fit the SED of the TIC435850195 system. The SED fitting technique is based upon that presented in Kurtz et al. (2020) and Rappaport et al. (2021, 2022); the MCMC implementation itself is derived from that used in Jayaraman et al. (2022b). For a given trial pair of masses and a system age, we obtain the radius, effective temperature, and luminosity from the MIST evolutionary tracks. The luminosity value is an input into the light curve calculation, as discussed previously; the radii and effective temperatures are used to estimate the emitted flux at each wavelength for which we have observations, using the Castelli & Kurucz model atmospheres. 3 The best fit SED is shown in Figure 4, and the best-fit stellar parameters are enumerated in Table 3.  \nWe find that the best-fit parameters yield an SED that agrees quite well with observations, but slightly underpredicts the two bluest points. However, we note that the GALEX FUV and NUV passbands are rather broad (see, e.g., shaded areas in Fig. 1 from Martin et al. 2005, and the horizontal error bars for the GALEX points in Fig. 4). Our analysis also shows that the secondary contributes a very small fraction of the overall system light, and is fainter by at least two orders of magnitude.', '3.3. System Parameter Results': "Table 3. Best-fit parameters from a joint MCMC fit to the SED and light curve. Uncertainties are given as 1σ confidence intervals, calculated from the posterior distributions shown in Figure 6. Parameters indicated with an * were directly varied as part of the MCMC run. The remaining parameters were uniquely determined for a given mass and age by its position on the MIST evolutionary tracks (Dotter 2016; Choi et al. 2016). The † indicates that we ran one MCMC with the distance and extinction as free parameters to ensure that our result was sensible, and then fixed them to the respective Gaia-determined values for future runs.  \nFigure 4. The best-fit SED. The yellow and green lines are the best-fit model atmospheres from Castelli & Kurucz (2003) for the primary and secondary, respectively, while the thicker purple line is the sum of the two models. At nearly all wavelengths, the primary star's emission dominates that of the secondary's by at least two orders of magnitude. We have added a horizontal error bar to show the width of the GALEX FUV and NUV bandpasses, and an upper limit for the measurement in the W4 band. The uncertainties of ∼ 10% on each flux point are too small to be seen at this scale.  \n<!-- image -->  \nThe final results of our joint fit to the orbital light curve and the composite SED are summarized in Table 3. The pulsating primary star has a mass of 1.8 M ⊙ and radius of 2.2 R ⊙ , and is somewhat evolved off the  \nFigure 5. Hertzsprung-Russell diagram showing evolutionary tracks for stars with masses between 0.7 and 1.9 M ⊙ from the MIST tracks (Dotter 2016; Choi et al. 2016). This diagram also shows the locations of the primary and secondary components, as stars and squares, respectively, for both the TIC184743498 (Z24) and TIC435850195 systems.  \n<!-- image -->  \nzero-age main sequence. The secondary is a K star of mass 0.66 M ⊙ that contributes under 1% of the system light. The orbital separation is close to 7 R ⊙ . Using the approximation for the effective Roche lobe (RL) radius r L from Eggleton (1983), with a mass ratio q of 0.37, we find r L = 0.59 a . With a ≃ 7 R ⊙ , we find that the primary star has r L = 4.13 R ⊙ ; thus, it fills slightly over 50% of its Roche lobe. The orbital inclination angle, a key parameter for this system, is 73.3 +0 . 6 -0 . 5 · .  \nWe note a potentially important discrepancy in what our system parameters predict for the K 1 velocity of the primary, compared to the measured value from Gaia, 46.13 km s -1 . Our best-fit system parameters indicate a K 1 value of 67 . 41 ± 4 . 23kms -1 , which differs from the Gaia value by 5.0 σ . To decrease this discrepancy to 2σ , either the inclination angle must be lowered to be between 35.6 · -46.9 · -which would make the eclipses disappear-or the mass of the secondary star must be lowered to M 2 ≲ 0 . 5 M ⊙ . Such a star, if near the ZAMS, would have an insufficient intrinsic luminosity to produce a secondary eclipse as deep as the one that is observed, which has a depth of approximately 1%; the eclipses produced for such a low-mass star would have a depth of ∼ 0.2%. We posit that the Gaia uncertainty on their measurement of K 1 is underestimated by a factor of 2-3.  \nOur best fit for the system age suggests that it is approximately 970 Myr old. The primary of this system appears to be slightly more evolved than the one in TIC184743498 (the first TAP to be discovered); however, its mass and other parameters are very similar to the primary in that system (cf. Table 3 in Z24). A comparison of the evolutionary states of the two systems is shown in  \nFigure 5, alongside evolutionary tracks for various stellar masses, assuming a solar metallicity. Both the TIC 184743498 and TIC 435850195 systems have R 1 /a close to 0.3 and, as we will see, have similar ranges for their δ Scuti pulsations. In TIC 184743498, however, the secondary is twice as massive as the one in TIC 435850195, and there also likely exists a tertiary component in that system that has comparable mass to the secondary.  \nPosterior distributions for a selection of the parameters tracked throughout the MCMC are shown in Figure 6.", '4. PULSATIONAL ANALYSIS': "To analyze the pulsations, we utilized Period04 (Lenz & Breger 2005) in order to perform a multifrequency fit for the amplitudes and phases at the time of a reference primary eclipse for all the pulsation modes.  \nWe first performed a standard frequency analysis, identifying the most significant pulsations. We do this by finding the highest peak in the Fourier amplitude spectrum, carrying out a linear least squares fit for the amplitude and phase of that frequency, and then subtracting it from the light curve. We repeat this process and sequentially remove the highest Fourier peaks until the noise floor is reached. From this initial list of frequencies, we constructed an 'echelle diagram (Figure 7), where the pulsation frequency ν puls is plotted against the ''echelle phase', which is defined as ( ν puls mod ν orb ) /ν orb .  \nWethen proceeded to identify frequencies in this 'echelle diagram that could be separated by integer multiples of the orbital frequency. Period04 is able to provide optimal light-curve fits for multiperiodic signals including harmonic, combination, and equally spaced frequencieswhich is essential for the present analysis. We used this option to test whether or not a given candidate multiplet of frequencies is consistent with the assumption of spacings that are related to the orbital frequencies; if not, their amplitudes would notably change using the force-fitted frequencies. Having identified all possible multiplets, we fixed all their frequencies, and calculated the amplitudes and phases for each pulsation frequency via least-squares fitting of each set of doublet frequencies to the light curve. All the frequencies in multiplets, and their amplitudes and phases, are enumerated in Table 4.  \nIf a pair of pulsations share the same 'echelle phase and are separated in frequency by 2 ν orb , we call that a 'dipole doublet.' A set of three pulsations separated by either ν orb or 2 ν orb (sets A and B in Figure 7) is referred to as either a dipole or a quadrupole triplet, respectively. Dipole doublet modes that have the same pulsation axis are circled in the same color in Figure 7.  \nTable 4. Multifrequency solution for all the observed pulsation multiplets, and the two highest-amplitude singlets. Multiplets are labeled according to their numbering in the 'echelle diagram (Fig. 7) and the amplitude-phase reconstructions (Figs. 8 and 9). We also provide mode identifications for the listed modes; the doublets and singlets are all Y 10 modes, while the mode identifications for the triplets are much more difficult to ascertain. Uncertainties for the frequencies range from 0.0002 d -1 for the strongest peaks, to 0.0015 d -1 for the weakest peaks that are still confidently detected.  \nFigure 6. Posterior distributions for stellar parameters from the joint MCMC fit to the light curve and the SED. The top panel shows the mass, radius, and temperature of the primary (in gray), as well as the inclination (in tan); the bottom panel shows the mass, radius, and temperature of the secondary (in magenta), and the age of the star (in turquoise). The median value is denoted by the solid vertical black lines, while 1σ confidence intervals are denoted by the dotted black lines. Note that even though we fit for the age indirectly in terms of the 'Equivalent Evolutionary Phase' (EEP), we show the posterior for the age in Gyr. Prior to the histogram calculation, the ages were normalized by the density of EEP samples in that particular age range.  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 7. An 'echelle diagram for TIC 435850195, showing the tidally tilted δ Scuti pulsations. The abscissa is labeled by the 'echelle phase, which is the pulsation frequency ν , modulo the orbital period, divided by the orbital period. The sizes of the points are scaled linearly according to their amplitudes in the periodogram of the light curve. From our analysis, we were able to extract 16 multiplets of tidally tilted pulsations. The doublet dipoles have been circled based on their pulsation axes. The numbered dipole doublets are shown in Figures 8 and 9; the triplets A and B are further discussed in Section 4.2.2.  \n<!-- image -->  \nTo identify the ℓ and m values for a given pulsation, we calculated its amplitude and phase as a function of the orbital phase. This was done following the formalism from Jayaraman et al. (2022a) for the special case of a doublet. The amplitude C and phase ϕ mult of a doublet, as a function of orbital phase Φ, can be written as  \nC 2 (Φ)= A 2 + B 2 +2 AB cos( ϕ A -ϕ B -2Φ) (5)  \nϕ mult = ATan2 [ A sin( ϕ A -Φ)+ B cos( ϕ B +Φ) A cos( ϕ A -Φ)+ B cos( ϕ B +Φ) ] (6)  \nHere, A and B are the amplitudes of the pulsations that comprise the multiplet, and ϕ A and ϕ B are their phases at a particular epoch time t 0 , which we have defined as the time of the first primary eclipse in the data (see Section 2). These amplitude and phase variations are shown as a function of orbital phase in Figures 8 and 9.  \nThe doublets with an amplitude maximum at the time of primary eclipse exhibit π phase shifts at orbital phases 0.25 and 0.75 (i.e., times of ELV maxima). These correspond to modes 1-6, as well as modes 13 and 14, in Figure 7. We identify these as Y 10 x , modes following the notation in Z24, where ℓ = 1, m = 0, and the x -axis corresponds to the tidal axis. In this case, we are looking  \ntoward the pulsation poles of these modes during eclipse, and the phase is expected to flip by π radians when we, the observers, switch to viewing the other pulsational hemisphere at orbital quadrature. These modes are similar to the ones presented in Figure 2 of Reed et al. 2005 (henceforth shortened to R05). 4 Their predictions for i = 75 · , and I p = 90 · -when the pulsation axis is pulled fully into alignment with the line of apsides with the binary-agree with our observations for these modes.  \nIn contrast, the other dominant set of modes from the 'echelle, numbered 7-12, exhibit amplitude minima at the times of primary and secondary eclipse, and π phase shifts at these times. Thus, these modes exhibit exactly the same behavior as the Y 10x modes, except they are shifted by 90 · in orbital phase. This similarity with the Y 10 x modes extends to the fact that they are also doublets that show no significant central peak. As was done in Z24, we tentatively identify these pulsations with  \nFigure 7 shows a number of other pulsation modes in addition to the dipole doublets discussed previously; these include singlets, as well as two potential tripletsone a dipole (B), and the other a quadrupole mode (A). Most of the pulsations that we see in TIC 435850195 can be explained within the tidal tilting framework, including some that may only be partially tilted ('tidally perturbed'). For the purposes of our analysis, we focus  \n<!-- image -->  \nFigure 8. The amplitude (left) and phase (right) variations of the six strongest Y 10x modes in TIC 435850195, as a function of orbital phase. These modes exhibit amplitude maxima at primary eclipse (phase 0, 1, and 2), and amplitude minima at the maxima of the ELVs (phases 0.25, 0.75, 1.25, and 1.75). These modes also undergo a π phase shift at minimum amplitude. The abrupt 2 π jumps in phase are an artifact of the cyclic nature of the phases and are not physically meaningful.  \n<!-- image -->  \nY 10 y modes. The description of what the observer sees as a function of orbital phase is exactly the same as for the Y 10 x modes, except shifted by 90 · in orbital phase. In other words, this is a Y 10 mode with a pulsation axis lying along the system's y axis, perpendicular to the axes of both the tidal bulge and the orbital momentum.  \nIf we were to adopt the formalism of R05, the amplitude variations of these modes, at first glance, appear most similar to Y 11x modes (cf. their Figure 2.2, with i = 75 · , and I p = 90 · ). However, we highlight that the corresponding R05 periodogram shows a detectable central peak that we do not observe in these modes of TIC 435850195. We suggest that the explanation given in Section 6 of Z24 reveals the true provenance of these modes, and identify these as Y 10y modes.  \nTo further disprove the hypothesis that the observed doublets 7-12 arise from Y 11x modes, we simulated the frequency spectrum expected from these modes as a function of orbital inclination angle i (shown in Figure 10). This figure, similar to Figure 3 from R05, shows that Y 11 modes must show a significant central peak for a binary with i ≲ 75 · . However, there exists no evidence for such a central peak in any of the doublets whose amplitude and phase variations are shown in Figure 9.  \nAs part of our verification procedure to confirm that the doublets do not arise from Y 11x modes, we also calculated the ratios of the central peak to the sidelobes that would be expected for these modes, as a function of inclination angle i . This calculation was based on the peaks shown in Fig. 10. In order to estimate the  \nnoise level, we fit a Rayleigh distribution to the white noise in the periodogram between 60 and 80 d -1 , which is a region without any significantly-detected frequencies. The 97% upper limit for this distribution, which has a mode of 18.74 ppm, was 50 ppm; below this value, we are unable to claim a detection of any periodicity in the light curve, even at marginal significance.  \nUsing this value, we set empirical upper limits on the presence of any central peak for all of modes 7-12, and thereby establish secure upper limits on the ratio of any putative central peak to its sidelobes (shown in Figure 11). For an inclination angle of ∼ 75 · , we find that the ratio of the central peak to the sidelobe should be slightly over 0.5; the limit for the observed modes, at the observed inclination angle, is under 0.4 for two of the modes, and under 0.1 for three of the modes. This further reinforces the fact that the modes we are seeing are, in fact, not Y 11x modes, and instead are the novel Y 10y modes discovered and reported by Z24.", '4.2. Other Observed Pulsations': '<!-- image -->  \nFigure 9. As Figure 8, but for the six highest-amplitude Y 10y modes in TIC 435850195. These modes have amplitude minima at primary eclipse (phase 0, 1, and 2), and amplitude maxima at the ELV maxima (phases 0.25, 0.75, 1.25, and 1.75). These modes also undergo a π phase shift at minimum amplitude. In every way, they are identical to the Y 10x modes, except for being shifted by 90 · in orbital phase. The phases for mode 12 have uncertainties of 0.243 and 0.125, respectively, which may partially explain why the calculated pulsation amplitude minima and maxima are slightly offset from the times of eclipse and ELV maxima.  \n<!-- image -->  \nFigure 10. Simulated Y 10y (left panel) and Y 11x modes (middle panel) for a star similar to TIC 435850195. These show that at the inclination angle of the system, roughly 75 · (indicated in blue), we expect to see a central peak that is approximately half the amplitude of the sidelobes for any Y 11x modes that may be present. In contrast, we observe only a frequency doublet of pulsations, regardless of inclination angle, for a Y 10y mode. The right panel shows a vertical zoom in on the periodogram for the Y 11x modes to highlight the Fourier peaks arising from eclipse mapping (also called spatial filtration), further discussed in Section 5.3. Note the differing x-axis scales between the middle and right panels.  \n<!-- image -->  \non the highest-amplitude singlets, and those modes that are in clear multiplets separated by multiples of ν orb .', '4.2.1. On the Existence of Radial Modes': "The 'echelle diagram from Figure 7 shows mostly dipoles and one quadrupole mode, but there are two prominent  \npulsation singlets: one 5 at 53.01285 d -1 (#16), and another one at 56.18624 d -1 (#15). At first glance, these  \nFigure 11. A comparison of the ratios of the central peak to the sidelobe amplitudes for Y 11x modes to the limits obtained from the Y 10y modes in TIC 435850195. The red lines correspond to the ratio of the 2σ detection threshold to the mean amplitude of each Y 10y doublet. The blue shaded area corresponds to the 1σ uncertainty for the inclination angle of the system, from our light curve fit.  \n<!-- image -->  \nmay appear to be radial modes. However, this interpretation is most likely incorrect.  \nWe first calculate the period ratio of these two oscillations, along with the large frequency spacing, to evaluate this hypothesis. Using the models of Stellingwerf (1979), we find that the ratios of the periods of the first, second, and third overtones to the fundamental period are 0.772, 0.629, and 0.527, respectively. The period ratio of these two modes is 0.944, so it is unlikely to be a low-order radial mode. For higher-order radial modes, the asymptotic (large) frequency spacing is given by  \n∆ ν = ∆ ν ⊙ √ ¯ ρ ¯ ρ ⊙ . (7)  \nHere, ∆ ν ⊙ is a constant, 135 µ Hz, and ¯ ρ is the mean stellar density. Utilizing our mass and radius values from Table 3 yields a stellar density of 0.249 g cm -3 . Substituting this into Equation 7 yields ∆ ν ≈ 57 µ Hz, or 4.9 d -1 , which is larger than the separation between these two singlet frequencies. Consequently, radial modes probably cannot explain these two observed singlets.  \nThese pulsations, being singlet modes, do not exhibit the amplitude and phase variability that we observe in the Y 10x and Y 10y modes. We thus interpret at least one of these modes, or possibly both, as non-radial Y 10z modes, similar to the framework presented in Z24, making TIC435850195 a novel tri-axial pulsator. We also suggest the possibility that one of these modes (but not both) could be a radial ( Y 0 ) mode instead of Y 10z .  \nTriplet A (as numbered in Figure 7) has each of its components separated by 2 ν orb (making it a quadrupole mode), while triplet B has a separation of just ν orb , making it a dipole mode. Using the formalism from Jayaraman et al. (2022a) to reconstruct the amplitude and phases of these triplets, as we did for the other modes in Subsection 4.1, we find that triplet A has two amplitude maxima per orbit; these exactly correspond to the ELV maxima. However, there are no π phase jumps, as would be expected for an m = 0 mode whose pulsation axis is the z -axis; moreover, given that the components are separated by 2 ν orb , this is likely to be a ℓ = 2 mode. While there is not an exact match from R05 for a quadrupole mode with a strong central peak and sidelobes separated by 2 ν orb , we note that there could be additional multiplet components that are obscured by the noise, making a conclusive mode identification difficult with the existing data set.  \nThere are two other possibilities that we could explore for these modes: Either it is an ℓ = 2 mode, with m = ± 2, or it is rotationally split; this latter hypothesis will be discussed further in Section 4.2.3. We note that R05 indicate that for ( ℓ, m ) = (2 , ± 2), there should be no π phase jumps under any condition. However, our inclination angle of i ∼ 73 . 5 · should yield a quintuplet spaced by ν orb for a fully tidally tilted mode ( I p = 90 · ). There are no clear matches from R05 for this scenario, even for tidally perturbed modes ( I p < 90 · ). We also note that it is possible for ℓ = 2 modes with different m = -2 , 0 , 2 to couple and produce what we refer to as Y 2+ and Y 2 -modes; however, these modes are predicted to have two π phase jumps for each orbital cycle, which we do not observe here. Excluding these possibilities suggests that this triplet is difficult to interpret with our current understanding of tidally tilted pulsations.  \nTriplet B has its components separated by ν orb ; each component frequency has a comparable amplitude. This case could correspond to the I p = 30 · case for a Y 10 mode, making this a tidally perturbed mode rather than a fully tilted mode. Given that the separation between the triplet components is exactly ν orb , we rule out the rotational splitting explanation (see Section 4.2.3). Another possibility is that this pulsation is tidally perturbed as part of a higher-order mode (e.g., with ℓ = 2), with the other frequencies being indistinguishable from the noise (cf. the I p = 30 · case, for i = 75 · in Fig. 3.4 in R05). This triplet could also arise from the coupling of modes with different ℓ values; however, that calculation is beyond the scope of this work, and will be further explored in Fuller et al. (in prep).  \nWe note that neither triplet is likely to be a Y 11x mode, based upon the ratio of the amplitudes of the  \npeaks shown in the simulated periodograms in Figure 10. Specifically, the central peak is the strongest component of the multiplet, when for Y 11x modes, it is usually the weakest, when observed at i ∼ 75 · .", '4.2.3. Ruling Out Rotational Splitting': "Rapidly-rotating stars can have their pulsations split by the rotational frequency, yielding multiplets whose components are given by the following expression:  \nν n,ℓ,m = ν n,ℓ + m (1 -C n,ℓ ) ν orb , (8)  \nwhere C n,ℓ is the Ledoux constant. For the observed 14 doublets that we identified as Y 10x and Y 10y modes, there are two compelling arguments against these arising from rotational splitting. First, the splitting of the doublets is equal to 2 ν orb to within an rms fractional uncertainty of ≲ 0 . 001. Therefore, the Ledoux coefficient would have to be ≲ 0 . 001, which is implausible-for δ Scuti stars, this value is between 0 and 0.2 (Goupil et al. 2000), with a specific value of 0.08 for the δ Scuti star KOI-976 (Ahlers et al. 2019). Second, the phase shifts between the doublet elements in the rotational-splitting scenario can, in principle, take on any value, including π . However, the likelihood that all of the doublets would have π phase shifts between their elements at exactly the same orbital phase is very low, if not zero. Thus, rotational splitting is categorically not the origin of the observed doublets.  \nFor triplet A, despite the difficulty in conclusively assigning a mode to it, we suggested that it may be an | m | = 2 mode. As a result, C n,ℓ must again approach 0 to within one part per thousand in order to explain this mode as rotational splitting. A similar argument for m = 1 holds for triplet B. As a result, we can eliminate the rotational splitting hypothesis as a viable explanation for what we observe in this star's pulsation spectrum.", '5. DISCUSSION': 'We have conclusively shown that TIC 435850195 pulsates along three different orthogonal axes. The system exhibits almost exclusively dipole modes, with one likely quadrupole mode. In this section, we introduce a toy model for mode coupling in the TIC 435850195 system, and further discuss the multiplets and the effects of eclipse mapping on the observed periodogram.', '5.1. Toy Model for Tidal Tilting': "As discussed in Z24 and in greater detail by Fuller et al. (in prep), the perturbation T to the potential from the tidal bulge in a close binary can be represented by  \nT ∝ x 2 r 2 , (9)  \nwhere ˆ x is the direction along the tidal bulge, ˆ z is the orbital angular momentum axis, and ˆ y is a direction perpendicular to both. In turn, this can be written in terms of spherical harmonics as  \nT ∝ Y 22 ,z + √ 2 3 Y 20 ,z + Y 2 -2 ,z . (10)  \nHere, the z subscript indicates that angles in these spherical harmonics are measured with respect to the z axis. As discussed briefly in Z24 and originally in Fuller et al. (2020), a perturbation analysis suggests possible mode couplings that the tidal bulge can induce. Equation (10) indicates that the Y 22 ,z term in the tidal perturbation can couple the mode Y 11 ,z with Y 1 -1 ,z (with ∆ m = +2), and vice versa, due to the Y 2 -2 ,z part of the tidal perturbation. The Y 20 ,z term in the tidal perturbation couples to none of the ℓ = 1 modes. Since the ℓ = 1 modes represent the vast majority of what we detect in TIC 435850195, we do not consider the coupling of ℓ = 2 modes-which will be discussed in Fuller et al. (in prep).  \nThe tidal perturbation thus gives new eigenmodes which are the sum and difference of the normal ℓ = 1 modes. Here, we compare the unperturbed dipole modes ( D 1 ± 1 ) to the tidally perturbed ones ( D ± ):  \nD 1 ± 1 ∝ Re { Y 1 ± 1 ,z e iω 1 t } ∝ sin θ cos( ω 1 t ± ϕ ) D + ∝ Re { ( Y 11 ,z + Y 1 -1 ,z ) e iω + t } ∝ sin θ sin ϕ sin ω + t D -∝ Re { ( Y 11 ,z -Y 1 -1 ,z ) e iω -t } ∝ sin θ cos ϕ cos ω -t  \nHere, ω + and ω -are the two perturbed eigenfrequencies, and Re {} represents the real part of the complex exponentials. The fundamental difference between the unperturbed ℓ = 1 modes and the new perturbed modes is that the unperturbed modes are waves that travel around the equator of the star, while the perturbed modes are standing waves, much like a Y 10 mode. In fact, the perturbed modes can be written explicitly as 'tidally tilted' modes with axes along the x and y axes:  \nD + ∝ y sin ω + t ∝ Y 10 y sin ω + t (11)  \nD -∝ x cos ω -t ∝ Y 10 x cos ω -t (12)  \nThese are exactly Y 10 modes with pulsation axes along the y and x directions, respectively. They describe the 14 tidally tilted dipole doublet modes we have detected in TIC 435850195, shown in Figures 7, 8, and 9.", '5.2. Other Possible Multiplet Groupings': "We note that the interpretation of certain groupings of peaks in the 'echelle (Fig. 7) can be subjective. Specifically, we discuss mode 16, which we choose to interpret as a singlet, and modes 1 and 4, which we choose to  \ninterpret as a doublet Y 10x mode, even though there is a very weak central component present.  \nFirst, we discuss the alignment of the frequency doublet (50.407261, 52.60155). This is split by 3 ν orb , which does not correspond to any clear mode identification in the tidal tilting framework. Comparison with the modes from R05 also does not yield a clear match for two significant peaks separated by 3 ν orb . This could either be an octupole mode, with the third component (also spaced by 3 ν orb ) hidden below the noise. An alternative (and perhaps more likely) explanation is that this is an ℓ = 3 mode, which might be able to produce doublets that are spaced by 3 ν orb in a triaxial tidally tilted pulsator (Fuller et al. in prep). We also note that highℓ modes are often strongly geometrically canceled (Dziembowski 1977), which could make this a chance alignment.  \nWhat we have tentatively called a singlet, Mode 16, also has another companion peak with a much smaller amplitude separated by -ν orb . However, we note that such a small peak does not induce any significant amplitude and phase variability in this mode; rather, this could be part of a triplet, with the peak at + ν orb buried in the noise. In that case, this could plausibly be interpreted as a dipole triplet. However, we would need data with a lower noise floor in order to extract this frequency and conclusively assign a mode identification to it, using the amplitude and phase variations throughout the orbit.  \nWe opted to interpret modes 1 and 4 as doublet Y 10x modes, despite the presence of a weak ( ≲ 10-20%) central peak in both of them. Comparing these to Figure 3 of R05 suggests that these, instead of being fully tidally tilted pulsations ( I p = 90 · ), may actually be tidally perturbed, or only partially tilted Y 10x modes. The presence of the weak central peak is a close match for the I p = 60 · and I p = 75 · cases, suggesting that these modes may not have their axis fully tilted into the orbital plane. If this were true, this star would have a unique combination of tidally perturbed modes, tidally tilted modes, and non-tidally tilted (singlet) pulsation modes, making the pulsation behavior far richer than initially thought.", '5.3. Further Asteroseismic Constraints': "This star exhibits 14 dipole doublet pulsations, 2 singlet pulsations, and 2 triplet pulsations-one dipole, and one quadrupole. We suggest that the triplet dipole mode may not be fully tidally tilted, and that the quadrupole mode is difficult to interpret as part of the tidal tilting framework. Some of these modes may be related to the coupling of ℓ = 2 modes, or coupling between modes that have different ℓ values-a derivation of which is beyond the scope of this work, and will be addressed by Fuller et al. (in prep). These modes may also be strongly affected  \nby the Coriolis force in the star, which has been neglected in our modeling of previously-discovered TTPs.  \nAn explanation that has sometimes been posited for the multiplets that are observed in 'echelle diagrams is 'eclipse mapping,' or 'spatial filtration' (see, e.g., Johnston et al. 2023, and references therein, including Gamarova et al. 2003). This phenomenon occurs when the occultation of the pulsating star obscures different parts of the surface of a star that is undergoing nonradial pulsations. This leads to considerably complicated variations in the amplitudes of pulsation multiplets (see, e.g., Fig. 10 of Johnston et al. 2023). However, we note that the periodogram of pulsations that have been spatially filtered will exhibit a long series of peaks split by ν orb that grow more pronounced as the inclination angle increases, as seen in the right panel of Figure 10. This effect is most pronounced for binary systems with i ∼ 90 · , and becomes nearly imperceptible (given the noise properties of the data) for i ∼ 75 · . Even the strongest peaks arising from eclipse mapping are only ≲ 5% of the Y 11x triplet. Thus, we can be confident that the observed multiplets are indeed caused by tidal phenomena, and are not simply a function of our observational perspective on the star.", '6. CONCLUSIONS': "In this work, we report the discovery of the second ever tri-axial pulsator, exhibiting Y 10x , Y 10y , and Y 10 z modes. Given the richness of the observed pulsational behavior, this star represents a unique laboratory through which we can investigate the effects of a companion's gravitational field on the pulsations of a star. We have also simultaneously fit the SED and the binary lightcurve for the system parameters using an MCMC algorithm. This analysis showed that this system consists of a slightly evolved primary δ Sct star, with a secondary K-type star that is still on the zero-age main sequence.  \nWe find that our best-fit inclination angle of 73 . 3 · ± 0 . 6 · is well outside the expected regime for the observed doublets to be Y 11x dipoles. This provides a useful way to rapidly confirm future tri-axial pulsators; we also provide a toy model framework in which this system and others can be interpreted in terms of tidal tilting. A targeted search for similar pulsators is underway, with a specific focus on ellipsoidal binaries (Green et al. 2023) and eclipsing binaries that lie in the δ Sct instability strip. Such a search does not solely rely on the presence of eclipses, allowing us to test our predictions even in the absence of eclipses, as well as for a large range of orbital inclination angles. As the TESS mission releases more light curves at 200 s cadence, we will be able to identify tri-axial pulsation in many more classes of stars.", 'ACKNOWLEDGEMENTS': 'RJ would like to thank Kevin Burdge for discussions about light curve fitting using Lcurve (Copperwheat et al. 2010), which we (unfortunately!) did not end up using as part of this work. RJ would also like to thank Michael Fausnaugh for information about generating barycentric corrections for TESS light curves.  \nV. B. K., S. R., and B. P. acknowledge financial support of the NASA Citizen Science Seed Funding Program, grant number 22-CSSFP22-0004. GH thanks the Polish National Center for Science (NCN) for supporting this study through grant 2021/43/B/ST9/02972.  \nThis 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 was obtained from MAST ( 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.  \nThis 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.  \nFacilities:  \nTESS, Gaia  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018, 2022), lightkurve (Lightkurve Collaboration et al. 2018), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), synphot (STScI Development Team 2018), TICA (Fausnaugh et al. 2020)', 'REFERENCES': 'Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2 von Zeipel, H. 1924, MNRAS, 84, 665,  \ndoi: 10.1093/mnras/84.9.665  \nWright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868  \nZhang, V., Rappaport, S., Jayaraman, R., et al. 2024, MNRAS, 528, 3378, doi: 10.1093/mnras/stae010', 'A. NUMERICAL INTEGRATION FOR LIMB-DARKENING': "As part of our light curve fitting code, we implemented a one-dimensional numerical integral to account for limb darkening. For each step in the light curve, we calculated the distance between the two stars' centers s , and the distance from the primary star's center to the limb of the secondary ( s -r 2 ). We then split the range from [ s -r 2 , r 1 ] into 75 steps, which corresponded to 75 circles concentric with the primary, with radii lying within that range; we denote these radii by r conc . For each of the circles, we calculated the distance from the center of the primary star to the chord between the circle's two points of intersection with the secondary's limb. This value is given by the expression  \nd = s 2 -r 2 2 + r 2 conc , x 2 s , (A1)  \nwhere x corresponds to a given concentric circle.  \nThis value enabled us to calculate the angle θ subtended by the chord between the two points of intersection, as θ/ 2 = cos -1 ( d/r conc , x ). The expression for the area of the portion of an infinitesimally thin circular ring that subtends an angle θ is given by rθ dr . We then integrated this against the limb-darkened flux L ( r ), from [ s -r 2 , r 1 ]. The formal integral for the limb-darkened flux (of the primary) within the overlapping area is given by  \n∫ r 1 s -r 2 2 r cos -1 ( d r ) L ( r r 1 ) dr (A2)  \nThe radius r was normalized to r 1 in the argument of the LD function (given below), as the limb darkening formula typically assumes a circle with radius scaled to 1.  \nL ( r ) = 1 -c 0 (1 -√ 1 -r 2 ) -c 1 (1 -√ 1 -r 2 ) 2  \nThis was then normalized with the total limb-darkened flux across the surface of the star, and then subtracted from 1, to determine the fraction of flux from the primary that was visible. To calculate the total limb-darkened flux, we numerically integrated the following expression, with ∆ r = 0 . 005:  \n∫ r 1 0 2 πr dr L ( r r 1 ) dr (A3)  \nFor the secondary, whose flux was 1% that of the primary, we instead calculated the limb darkening at ( r 1 + r 2 -s ) / 2 and applied that across the entire overlapping area. We note that this approximation scheme breaks down when the center of the secondary is occulted by the primary star; however, given that this scenario only occurs for a small number of phase bins in our system, and we do not observe any significant residuals in Figure 3 as a result of this approximation scheme, we opt to ignore this higher-order correction. Obscuration of the stellar center would correspond to cases VII and VIII in the scheme of Mandel & Agol (2002)."}