Physical properties of the eclipsing binary KIC 9851944 and analysis of its tidally-perturbed p- and g-mode pulsations

Stars that are both pulsating and eclipsing offer an important opportunity to better understand many of the physical phenomena that occur in stars, because it is possible to measure the pulsation frequencies of stars for which the masses and radii are known precisely and accurately. KIC 9851944 is a double-lined detached eclipsing binary containing two F-stars which show both pressure and gravity mode pulsations. We present an analysis of new high-resolution spectroscopy of the system and high quality light curves from the Kepler and TESS space missions. We determine the masses and effective temperatures of the stars to 0.6% precision, and their radii to 1.0% and 1.5% precision. The secondary component is cooler, but larger and more massive than the primary so is more evolved; both lie inside the {\delta} Scuti and {\gamma} Doradus instability strips. We measure a total of 133 significant pulsation frequencies in the light curve, including 14 multiplets that each contain between 3 and 19 frequencies. We find evidence for tidal perturbations to some of the p- and g-modes, attribute a subset of the frequencies to either the primary or secondary star, and measure a buoyancy radius and near-core rotational frequency for the primary component. KIC 9851944 is mildly metal-rich and MIST isochrones from the MESA evolutionary code agree well with the observed properties of the system for an age of 1.25 Gyr.


INTRODUCTION
Binary and multiple systems make up the vast majority of all observed medium-and high-mass (M > 1.15 M ⊙ ) stellar objects in the galaxy (Duchêne & Kraus 2013).For double-lined spectroscopic binary systems (SB2s), where light emitted from both components is visible in the stellar spectrum, the radial velocity (RV) variations of each component throughout the orbit can be measured.Combining this with the analysis of the light variability due to geometrical effects when the components also eclipse one another, i.e., eclipsing binaries (EBs), can lead to estimations of the mass and radius to within 1% (Torres et al. 2010), making these objects our best source of precise fundamental stellar parameters.Since the methods used to derive the properties of the components in double-lined eclipsing binaries (DLEBs) do not rely on theoretical stellar models, these objects are critical for testing and verifying stellar evolution theory (Pols et al. 1997;Pourbaix 2000;Lastennet & Valls-Gabaud 2002;de Mink et al. 2007).
Such measurements for the masses and radii of stars allow for a precise determination of evolutionary status and age (e.g., Higl & Weiss 2017) and the calibration of interior physics.As examples, the mass dependence of convective core overshoot-⋆ E-mail: z.jennings@keele.ac.uk ing (Claret & Torres 2016, 2017, 2018, 2019), and the amount of core mixing in massive stars has been probed using DLEBs (Pavlovski et al. 2018;Tkachenko et al. 2020).
It has been known since the 1970s that stellar pulsations occur in the components of EBs (e.g., Tempesti 1971, for the system AB Cassiopeiae).Oscillation signals measured in the light curves of pulsating stars can be studied using asteroseismology.This is an alternative way by which accurate stellar parameters can be derived (Aerts et al. 2010) and is a powerful tool for probing the internal structures of stars (Aerts et al. 2017;Chaplin & Miglio 2013).The analysis of space-based photometry of pulsating stars has revealed information about important phenomena such as internal rotation, core overshooting and angular momentum transport (e.g., Saio et al. 2015;Lovekin & Guzik 2017).
It is advantageous to study pulsating stars in EBs because analysing the intrinsic variability due to oscillations of the stellar interior and atmosphere as well as the light variability during eclipses leads to independent constraints on theoretical models (Guo et al. 2019;Liakos 2021;Miszuda et al. 2022).The parameters derived from EBs can also be used to constrain the seismic models (e.g., Sekaran et al. 2020), and this can be useful, for example, in mode identification for δ Scuti stars which is difficult (Streamer et al. 2018;Murphy et al. 2021).The complimentary nature of these two methods of analysis makes pulsating stars in EBs (particularly DLEBs) the best objects to use to refine our knowledge of stellar structure and evolution (Guo et al. 2016).Lampens (2021) and Southworth (2021) have reviewed the impact of space missions on the study of EBs with pulsating components, and on binary stars in general, respectively.
δ Scuti stars are a class of short-period (15 min to 8 h; Uytterhoeven et al. 2011;Aerts et al. 2010), multiperiodic pressure mode (p-mode) pulsators.They are dwarfs or subgiants located at the lower end of the classical instability strip in the Hertzsprung-Russell (HR) diagram (Breger 2000;Dupret et al. 2005;Murphy et al. 2019).Their spectral types are A2 to F5, corresponding to typical effective temperatures (T eff s) of 6500 -9500 K (Liakos 2021).Their masses range from 1.5 to 2.5 M ⊙ (Aerts et al. 2010;Yang et al. 2021), which places them in the transition region between lower-mass stars with radiative cores and thick outer convection zones, and massive stars with convective cores and thin outer convective zones (Yang et al. 2021).They thus provide an opportunity to study the structure and evolution of stars in this transition region (Bowman & Kurtz 2018).Low-order, non-radial modes are generally observed for these stars and these modes are driven by the κ mechanism acting in the partial ionization zone of He II (Pamyatnykh 1999;Breger 2000;Antoci et al. 2014;Murphy et al. 2020).Higher-order non-radial pulsations have been observed in some δ Scuti stars such as τ Peg (Kennelly et al. 1998), where the κ mechanism may not be sufficient to explain the modes, and may instead be attributed to the turbulent pressure in the hydrogen convective zone (Antoci et al. 2014;Grassitelli et al. 2015).
The γ Doradus class of pulsators are also located at the lower end of the classical instability strip in the HR diagram; they exist near the red edge of the δ Scuti instability strip (Yang et al. 2021).Such stars pulsate with typical periods between 8 and 80 h (Handler 1999) in high order gravity modes (g-modes) (Kaye et al. 1999), believed to be excited by the interaction of convection and pulsations (e.g.Guzik et al. 2000;Grigahcène et al. 2005;Dupret et al. 2005).Typical masses are between 1.3 and 1.8 M ⊙ (Hong et al. 2022;Aerts et al. 2023) with spectral types of F5 to A7.The high radial order of their pulsations puts them in the asymptotic regime, meaning their oscillations are equally spaced in period for non-rotating, non-magnetic, chemically homogeneous stars (Tassoul 1980).Deviations from a homogeneous spacing emerge due to chemical gradients and rotation (Sekaran et al. 2020), where the former are influenced by the effects of diffusive mixing; mixing reduces the steepness of chemical gradients.Period spacing diagrams can therefore be used to derive information on chemical composition gradients, internal rotation rates and diffusive mixing processes (e.g.Bouabid et al. 2013;Bedding et al. 2015;Saio et al. 2015;Van Reeth et al. 2016;Ouazzani et al. 2017;Li et al. 2019;Miglio et al. 2008;Moravveji et al. 2015;Sekaran et al. 2021).
Space missions such as Kepler (Koch et al. 2010;Borucki et al. 2010), CoRoT (Baglin et al. 2006) and TESS (Ricker et al. 2015), have delivered a large amount of photometric data with precisions unachievable from the ground.The unprecedented precision has allowed for the detection of extremely low-amplitude frequencies (Murphy et al. 2013;Bowman & Kurtz 2018) while long sequences of continuous observations (Lehmann et al. 2013) means that longer-period pulsations, such as those typical of γ Doradus stars, can be studied.The overlap of the δ Scuti and γ Doradus instability strips supports the existence of δ Scuti/γ Doradus hybrids (Breger & Beichbuchner 1996;Handler & Shobbrook 2002;Yang et al. 2021).These were expected to be rare, based on early calculations by Dupret et al. (2005), but the lower detection thresholds provided by space missions has led to the discovery that such behaviour is indeed common (Uytterhoeven et al. 2011;Grigahcène et al. 2010;Bradley et al. 2015;Balona et al. 2015;Guo et al. 2019).Later calculations by Xiong et al. (2016) conform better with these findings.
Hybrids have great potential for asteroseismology (Schmid & Aerts 2016).The p-modes probe the stellar envelope while the g-modes carry information about the near-core regions (Yang et al. 2021;Grigahcène et al. 2010;Saio et al. 2015;Kurtz et al. 2014).The behaviour also means that information about two different driving mechanisms can be obtained (Hong et al. 2022).
KIC 9851944 is an EB showing δ Scuti/γ Doradus hybrid pulsation signatures.Therefore, the object is an ideal candidate for testing our understanding of stellar structure and evolution given the large amount of constraints that can be obtained due to the advantages associated with hybrid behaviour as well as binarity.KIC 9851944 is included in the Kepler Eclipsing Binary Catalogue (KEBC; Prša et al. 2011;Kirk et al. 2016), as well as a study by Matson et al. (2017) who presented RVs for 41 Kepler EBs.Gies et al. (2012Gies et al. ( , 2015) ) studied the eclipse times for KIC 9851944 and found no evidence of apsidal motion or a third body; the object was also included in a catalogue of precise eclipse times of 1279 Kepler EBs by Conroy et al. (2014).Guo et al. (2016) combined the analysis of Kepler photometry with medium-resolution spectra (R = λ /∆λ ∼ 6000) to determine the atmospheric and physical properties of KIC 9851944; we list these results in Table 8.Evolutionary modelling based on these properties shows the post-main-sequence secondary to be more evolved than the main-sequence (MS) primary.The authors concluded that both components show δ Scuti type pulsations, which they interpreted as p-modes and p and g mixed modes, and attempted to identify the modes by comparison with theoretically computed frequencies; the range of theoretically predicted unstable modes agreed roughly with observations but the authors conclude that mode identification is still difficult in δ Scuti stars, even with constrained mass, radius and T eff .
This work aims to be complementary to the work by Guo et al. (2016) on KIC 9851944; we additionally include observations by TESS in our photometric analysis and combine this with the analysis of high-resolution (R = 60000) spectroscopic observations.Section 2 outlines the photometric and spectroscopic observations.We determine the orbital ephemeris based on the photometric observations in Section 3. In Section 4 we analyse RVs derived from the spectroscopic observations and in Section 5 we present the spectroscopic analysis.The analysis of the photometric light curves is given in section 6 and in Section 7 we present the physical properties of the system.An investigation of the pulsations is given in Section 8.The discussion and conclusion are given in Sections 9 and 10, respectively.

Photometry
The Kepler space telescope (Koch et al. 2010) was launched in March 2009 with the primary aim of detecting Earth-like extrasolar planets around solar-like stars.However, the mission has also provided a huge amount of high-quality data on stars and stellar systems (Hong et al. 2022).The data were collected using two modes of observation (Gilliland et al. 2010)  solar arrays pointed toward the Sun during its Earth-trailing heliocentric orbit, the spacecraft was rotated by 90 • every three months (one quarter).Primary functions of the Kepler photometric analysis module of the Kepler science data processing pipeline are responsible for generating the Simple Aperture Photometric (SAP) flux time series of the observed objects, while secondary photometric analysis functions performed by the Pre-Search Data Conditioning (PDC) module support systematic error correction, giving rise to the PD-CSAP fluxes (see Jenkins 2017, for further details on the SAP and PDCSAP measurements).KIC 9851944 was observed by Kepler in six quarters (0,12,13,14,16,17) in short cadence mode between May 2009 and May 2013.
The Transiting Exoplanet Satellite Survey (TESS; Ricker et al. 2015) was also designed to discover extrasolar planets and also produced a large amount of high-quality data for many other celestial objects.More than 2 × 10 5 main-sequence dwarfs were observed at 2 min cadence during the all-sky survey, which was divided into overlapping sectors on the sky.KIC 9851944 was observed by TESS in sectors 14, 15, 41, 54, 55 and 56 between July 2019 and September 2022.The TESS science pipeline is based on the Kepler science pipeline so the architecture and algorithms are similar; the photometric analysis module computes the SAP flux time-series and the PDC module corrects for systematics to compute the PDCSAP fluxes (Jenkins et al. 2016).
The light curves from the Kepler quarters and TESS sectors mentioned above were downloaded from the Mikulski Archive for Space Telescopes (MAST) and are used in Section 3 to determine the ephemeris of the system, as well as Section 6 to obtain the final model of the light curve.Both SAP and PDCSAP measurements are available; the SAP and PDCSAP fluxes were similar, which we verified by inspecting the SAP fluxes with the PDCSAP fluxes overplotted after dividing them by their median flux values to put them on the same scale.Thus, we used the SAP measurements to avoid possible biases due to the additional processing applied to the PDC-SAP data.A second-order polynomial was fitted to fluxes that correspond to positions of quadrature, i.e., the maximum of the ellipsoidal brightening, to estimate systematic trends present in the light curves.
Subtracting the difference between this polynomial and the median flux of the light curve yields a model for the local median level of out-of-eclipse flux.This model was then divided out to remove systematic trends.The residual value of a smoothed version of the light curve subtracted from the observed light curve was calculated and observed fluxes that deviated by more than 1% were rejected.Fluxes were converted to magnitudes and errors were propagated following Prša (2018).The short cadence Kepler light curve from Quarter 0 and the two-minute cadence TESS light curve from Sector 15, after performing these operations, are shown in the top and middle panels of Fig. 1, respectively.KIC 9851944 was also observed by the Wide Angle Search for Planets (WASP) telescope (Butters et al. 2010) between May 2007 and July 2010.WASP consisted of two robotic observatories, one in the Nothern Hemisphere at Observatorio del Roque de los Muchachos on La Palma and the other at the South African Astronomical Observatory (SAAO), each with eight 20 mm telescopes on a single mount.Observations of KIC 9851944 collected by WASP between May 2007 and July 2007 as well as between June 2008 and July 2008 are shown in the bottom panels of Fig. 1 phase folded about the orbital period determined in Section 3. The data collected by WASP for KIC 9851944 were only used to constrain the times of primary minima when performing preliminary fits to the Kepler and TESS light curves in Section 3, to determine the orbital ephemeris.

Spectroscopy
A set of 33 spectroscopic observations were obtained for the target using the Hamilton échelle spectrograph (Vogt 1987) on the Shane 3 m telescope at Lick Observatory during two observing runs, one in July 2012 and the other in June 2013.The data were obtained using CCD chip no. 4, giving a wavelength coverage of 340 -900 nm over 89 échelle orders at a resolving power of R ≈ 60000.
The data were reduced using the standard pipeline for the spectrograph (e.g.Fischer et al. 2014).Flat-fields were obtained with a quartz lamp and divided from the spectra.The wavelength calibra- tion was obtained from exposures of a thorium-argon emission lamp taken roughly every hour during the night.Details for the normalization of the one-dimensional extracted spectra are given in Sections 4 and 5 because each of those components of the analysis used different approaches.
Table 1 gives the epochs of each of the 33 spectrosopic observations as well as the signal-to-noise ratio (S/N), estimated as the square root of the counts close to the peak of the best-exposed échelle order.Also given are the corresponding RVs of each component that are derived in Section 4.

ORBITAL EPHEMERIS
A first model of the Kepler and TESS light curves was obtained using the JKTEBOP code (Southworth 2013), which approximates stellar shapes in detached eclipsing binaries based on the simple and efficient NDE biaxial ellipsoidal model, so is computationally fast (Nelson & Davis 1972;Southworth et al. 2004).We fitted for the period P, epoch T 0 , surface brightness ratio J, sum of the fractional radii r A + r B , ratio of the radii k = r B r A , inclination i, the Poincaré elements e cos ω and e sin ω, and a light scale factor.We define star A to be the star eclipsed at the primary (deeper) minimum and star B to be its companion.
We also fitted for the linear coefficients u of the quadratic limb darkening (LD) law while quadratic terms were taken from Claret & Bloemen (2011) and Claret (2017) for the Kepler and TESS bands, respectively; we used the estimated atmospheric parameters reported in the KIC to choose the appropriate values for the quadratic coefficients.We performed fits to the WASP light curve but the lower quality of this data compared to the Kepler and TESS light curves means that including these results in the calculation of the overall preliminary eclipse model would lead to less well determined parameters.Thus, we simply include the epochs of primary minimum estimated from the WASP light curves as additional observational constraints on T 0 in the preliminary fits to the Kepler and TESS light curves.
The adopted values for the light curve parameters from this preliminary analysis were taken as the weighted averages of the results from the fits to the individual Kepler quarters and TESS sectors, where the reciprocal of the squared errors from the covriance matrix were used as weights.These values are given in Table A1, and Fig. 2 shows the fit to the Kepler Quarter 0 light curve; high-frequency variability due to pulsations of the δ Scuti type is clearly visible at all orbital phases in the residuals shown in the lower panel, with amplitudes of around ±0.2 mmag.
The value for r A + r B in Table A1 suggests that the components of KIC 9851944 are deformed beyond the limits of applicability of the biaxial ellipsoidal approximation (Popper & Etzel 1981) due to their close proximity and thus strong mutual deformations; the ellipsoidal variation (∼ 0.04 mag) observed in Fig. 1 and Fig. 2 is another indication of this.Furthermore, plotting the results from individual quarters and sectors against each other reveals strong degeneracies between the fitted parameters, specifically r A + r B , k, i, J and the LD coefficients.This is highlighted by the results giving a wide range of ratios of the light contributions of the two stars.Therefore, we present a detailed analysis of the Kepler and TESS light curves using a more sophisticated modelling code, as well as constraints on the light ratio ℓ B /ℓ A from our spectral analysis, to reliably determine the light curve solution in Section 6.
Among the reliably constrained parameters from this preliminary analysis are the period and epochs of primary minimum from individual quarters and sectors.Therefore, the analysis is useful for establishing the orbital ephemeris.A linear ephemeris was fit to the resulting values of T 0 against orbital cycle using the reciprocal of the squared errors from the covariance matrix as weights; these errors were rescaled during the fitting procedure to yield a reduced chi-squared value of χ 2 ν = 1.Fig. 3 shows the residuals of the fit and the red line represents the linear ephemeris.
There seems to be a trend in the O-C diagram compared with the linear ephemeris, with data points corresponding to TESS observations (>1000 cycles) all appearing above the red line.Thus, we also attempted to fit for a quadratic ephemeris, which is represented by the blue dashed line in Fig. 3.However, the corresponding quadratic term was a similar size to its errorbar.We therefore decided to stick with the linear ephemeris which was measured to be: Min I (BJD TDB ) = 2456308.304072(14)+ 2.163901775(98) E.

RADIAL VELOCITY ANALYSIS
The spectral range between 4400-4800 Å is a suitable region to carry out the RV extraction procedure because there are many wellresolved spectral lines and the region does not contain any wide Balmer lines.Thirteen échelle orders within this range were corrected for cosmic ray spikes by first resampling onto a homogeneous wavelength grid before using a median filter to smooth them, then resampling the orders back to the original grid and subtracting them from the observed spectrum (Blanco-Cuaresma et al. 2014a).Pixels where the residual deviated by a threshold percentage, typically between 5 and 10%, were masked and their values were corrected for using linear interpolation.The appropriate value for the threshold was investigated for each order individually based on visual inspections of resulting flagged pixels.The blaze signature was removed using the method of Xu et al. (2019), where alpha shape fitting (see, Edelsbrunner et al. 1983) combined with local polynomial regression is used to estimate the continuum without dividing out important spectral features.
Small errors in the continuum estimation are magnified where the blaze function approaches zero.Following Xu et al. (2019), a better estimation at the edges can be obtained by taking a weighted average in the overlapping regions and applying weights such that pixels further away from the edge within their own order are favoured.For this to be applied correctly, the orders were re-sampled onto grids with a common wavelength step so that the pixel averages were calculated for corresponding wavelength values between neighbouring orders.A useful result from the edge correction procedure is that neighbouring orders share a common edge and can be merged after truncating the orders in the centre of the overlapping regions.The following analysis was carried out after merging the 13 spectral orders with wavelengths between 4400-4800 Å using this process.
Preliminary estimations for the projected rotational velocity, v sin i, of the components of KIC 9851944 were obtained by cross correlating the observed spectra with a template broadened to a range of v sin i values between 0 and 150 km s −1 in steps of 10 km s −1 .Observations taken less than 0.125 orbital phases from the eclipses were excluded to minimise issues arising due to lines from each component blending or being eclipsed.Interpolating the peak heights of the correlation functions and taking the maximum yields v sin i A = 46.3± 0.4 km s −1 and v sin i B = 57.2± 0.5 km s −1 for the primary and secondary component, respectively.These values for v sin i were then used and applied to synthesise primary and secondary templates using the ISPEC code (Blanco-Cuaresma et al. 2014b).In the first iteration, the T eff and log(g) values were taken from the Kepler Input Catalogue (KIC).For the final iteration of the RV extraction, the atmospheric parameters of the templates corresponded to those derived from our atmospheric analysis of the disentangled spectra (see Section 5.1.)Our own implementation of the TODCOR two-dimensional crosscorrelation algorithm (Zucker & Mazeh 1994) was used to extract the RVs.Systematic errors arise when using the cross-correlation technique to extract RVs because neighbouring peaks and sidelobes in the doubled-peaked CCF disturb each other, and this is related to blending between the spectral lines of the two components (Andersen 1975;Latham et al. 1996).The result is that the estimated RVs may be shifted relative to the true value, where the size and direction of the shift depends on the phase (see Fig. 4), i.e, the relative Doppler shifts between the components (Latham et al. 1996).
This systematic error can be minimised by correcting the RVs for this effect (e.g., Latham et al. 1996;Torres et al. 1997;Torres et al. 2000;Torres & Ribas 2002;Southworth & Clausen 2007).A pre- liminary orbit was obtained by making an initial fit to the RVs extracted using TODCOR.A synthetic SB2 model was then produced with orbital parameters corresponding to the initial fit by Doppler shifting and adding the template spectra weighted by their relative light contributions as estimated by the resulting value for ℓ B /ℓ A from the initial TODCOR run.The TODCOR algorithm was then used to extract the known RVs from the synthetic SB2 model at orbital phases corresponding to the observations.The difference between the calculated RV and the known RV of the model gives an estimation for the magnitude and direction of the systematic error at that orbital phase.The correction was then added to the actual RVs measured from the observed spectra.This procedure was iterated multiple times using templates with different atmospheric parameters; the top panel of Fig. 4 shows the resulting orbital fit to the corrected velocities that were derived in our final iteration, where we used templates with atmospheric parameters corresponding to those derived in Section 5.1, and the results are presented in Table 2.The middle panel of Fig. 4 shows the residuals of the fit and the lower panel, split into two for the primary and secondary, gives the corrections that were applied to the extracted RVs.RVs with corrections larger than 2.5 km s −1 were excluded and are not shown because these RVs correspond to phases of conjunction, where line blending effects are most severe, and negligible information is contained on the velocity semiamplitudes.The fit was constrained to a circular orbit because attempts to fit for eccentricity yielded values consistent with zero and the study by Guo et al. (2016) suggests that this system has circularised.The uncertainties on the RV measurements were rescaled by a constant factor during the fitting procedure to yield a χ 2 r of unity.The root mean square (RMS) of the residuals of the fit compared to the RVs for the primary and the secondary are 1.1 km/s and 1.4 km/s, respectively.
The differences between our final results for K A and K B presented in Table 2 and those derived from the initial run (both after applying corrections), which used templates with atmospheric parameters corresponding to those from the KIC, was ∼ 0.1%.The corresponding difference in ℓ B /ℓ A was ∼ 9%.This shows that the TODCOR light ratio is more sensitive to the atmospheric parameters of the templates than the derived RVs.We added these differences in quadrature, to the uncertainties derived from the covariance matrix of the fit for the orbital parameters, and to the standard error in the mean value of ℓ B /ℓ A derived from different observations, in our calculation of the parameter error bars in Table 2.
Applying the corrections resulted in an increase in the velocity semiamplitudes K A and K B by 0.2% and 0.3%, respectively.This is a small increase, suggesting that TODCOR is less sensitive to blending effects than the one-dimensional cross correlation technique.However, the 0.2% and 0.3% increase in the velocity semiamplitudes translates to a 0.6% and 0.9% increase in the derived masses, which is significant considering that we aim to achieve precisions smaller than these values.This suggests that the corrections are necessary.Furthermore, the corrections are clearly systematic, and depend on the relative Doppler shifts of the components, as shown in the bottom panels of Fig. 4.

Atmospheric parameters
The method of spectral disentangling (SPD) makes it possible to separate the spectra of individual components from a time-series of binary or multiple star spectra (Simon & Sturm 1994).SPD has other advantages: the optimal orbital parameters are determined at the same time, and the disentangled spectra have a higher S/N than the original observed spectra.The disentangled spectra contain information on the atmospheric parameters (T eff , surface gravity, v sin i and metallicity) of each component.The disentangled spectra are still in the common continuum of the binary system, but can be renormalised if the light ratio between the components is known from other sources, e.g. from the analysis of light curves (Hensberge et al. 2000;Pavlovski & Hensberge 2005, 2010).Alternatively, the disentangled spectra can be fitted with synthetic spectra to determine the light ratio between the components, and this can be useful for partially-eclipsing systems where there is a degeneracy between the radius ratio and light ratio in the system (Tamajo et al. 2011;Tkachenko 2015;Pavlovski et al. 2023).In well-determined cases, the precision of the light ratio can approach 1% (Pavlovski et al. 2022).
Since we planned to apply the method of SPD to extract the individual spectra of the components, normalisation of the observed spectra was of critical importance.We used a different approach than that in Section 4, where we extracted RVs.Here, we used the dedicated code described in Kolbas et al. (2015).First, the blaze function of échelle orders was fitted with a high-order polynomial function.Then, the normalised échelle orders were merged.When overlapping regions of successive échelle orders are sufficiently long, the very ends were cut off because of their low S/N.Échelle orders containing broad Balmer lines, in which it is not possible to define the blaze function with enough precision, were treated in a special way.For these échelle orders, the blaze function was interpolated from adjacent orders.This produces more reliable normalisation in spectral regions with broad Balmer lines than the usual pipeline procedures.For recent applications of this approach, please see Pavlovski et al. (2018Pavlovski et al. ( , 2023)); Lester et al. (2019Lester et al. ( , 2022)); Wang et al. (2020Wang et al. ( , 2023)).
SPD was performed in the Fourier domain with the prescription by Hadrava (1995).The FDBINARY code, developed in Ilijic et al. (2004), was applied to the time-series of normalised échelle spectra of KIC 9851944.FDBINARY uses the fast Fourier transform approach which allows flexibility in selection of spectral segments for SPD whilst still keeping the original spectral resolution.The orbital parameters, specifically K A and K B , were fixed to the solution reported in Table 2, thus SPD was used in pure separation mode (Pavlovski & Hensberge 2010).At this point the reconstructed individual spectra of the components were still in the common continuum of the binary system.A portion of separated spectra for both components in the spectral range 4900-5290 Å is shown in Fig. 5.
Overall, a slight difference in line depth between the two components is seen.The more massive component has deeper lines and slightly faster rotation.The Balmer lines are broadened by Stark broadening, and generally are not sensitive to the rotational broadening.Thus, we first optimised portions of the disentangled spectra free of the Balmer lines, primarily to discern the v sin i values.We then performed optimal fitting of disentangled spectra centred on the Hβ and Hα lines, with fixed v sin i.
It is well-known that the hydrogen Balmer lines are excellent diagnostic tools for the determination of the T eff s for stars with T eff < 8000 K because the degeneracy with the surface gravity van-   ishes (Gray 2005).Moreover, we can use log(g) determined for the components since this quantity is determined with high precision from the analysis of DLEBs.In the case of KIC 9851944, the log(g) values are determined with uncertainties of about 0.01 dex, for both components.Therefore, preference is given to the determination of the T eff for the components in KIC 9851944 from line profile fitting of Balmer lines, with fixed surface gravity.This is advantageous over the calculation of the T eff from metal lines because their strength depends on the metallicity of the stellar atmosphere.
Without any appreciable changes in the light ratio in the course of the orbital cycle, i.e. no spectra were observed in eclipse, ambiguity exists in the reconstruction of the components' spectra and only separation of the spectra still in the common continuum of the binary system is possible.The components' spectra are correctly reconstructed but with scaling factors, i.e. the shapes of the spectral lines are correct, but not their strength.In such a case, the light ratio can be determined from fitting the separated spectra with synthetic (theoretical) spectra, in the course of determination of the atmospheric parameters.The optimal fitting was performed in constrained mode, as defined in Tamajo et al. (2011), with the condition that the sum of the light dilution factors is equal to unity.The STARFIT code (Kolbas et al. 2015) uses a genetic algorithm to search for the best fit within a grid of synthetic spectra pre-calculated using the UCLSYN code (Smalley et al. 2001).The uncertainties were calculated using the MCMC approach described in Pavlovski et al. (2018).This task was straightforward due to the fixed surface gravities and v sin i values: only the light ratio and T eff s were optimised.
The results of the optimal fitting for all three spectral segments are given in Table 3.We adopt the weighted average of the results from Hα and Hβ for T eff , which is 6964 ± 43 K and 6840 ± 37 K for the primary and secondary, respectively; these values are used to obtain the light curve solution in Section 6 and are presented formally in Section 7. Our T eff values indicate that the KIC T eff (6204 K) for this object is underestimated, consistent with prior studies (e.g., Molenda-Żakowicz et al. 2011;Lehmann et al. 2011;Tkachenko et al. 2013;Niemczura et al. 2015) reporting similar underestimations for stars with T eff ∼ 7000 K.
With the T eff determined from the Balmer lines, the depths of the metal lines cannot be matched -the metal lines for both components are slightly deeper than the synthetic spectra for given parameters.This could be explained with the metallicity effect, which affects metal lines more than broad Hβ and Hα lines for which Stark broadening dominates.The effects of metallicity explain the somewhat higher T eff (by about 200 K) found from the metal lines.From a grid of calculated synthetic spectra for given atmospheric parameters derived from Balmer lines, we found that the metallicity for both components is [M/H] ∼ 0.15 dex.Since the components of KIC 9851944 are at the cool end of where Am stars are found, this motivated us to disentangle the spectra between 3853 -4067 Å and compare the Ca K lines to synthetic spectra; the findings did not suggest any Am peculiarity although the Ca K line for the primary appears slightly weakened.Due to the relatively high v sin i for both components and thus severe line blending, as well as the low S/N of the individual spectra, we did not attempt to determine individual elemental abundances.
The light ratio from the optimal fitting of the wings of Hβ and Hγ lines, spectral segments containing only metal lines, and the Mg I triplet at around 5180 Å, are 1.315 ± 0.018, 1.304 ± 0.025 and 1.278 ± 0.033 (Table 3).All three values for the light ratio are consistent within their 1σ uncertainties, with the light ratio determined from the wings of Hβ line being the most precise of the three.graphically reconstructed spectra of the components.Their spectra cover the wavelength range 3930-4610 Å at a medium spectral resolution of R = 6000.The analysis by Guo et al. (2016) was similar to ours, but with an important difference that they fitted complete separated spectra, whilst we concentrate on the wings of Balmer lines.
Our results for the Balmer lines corroborate their findings to within 1σ .We also find good agreement for the v sin i values, which are also within 1σ for both components.Guo et al. (2016) obtained a mean light ratio 1.34 ± 0.03 from spectra centred at 4275 Å -this is slightly larger than but still consistent with our own results.

Direct fitting for the light ratio
We were initially unable to constrain the light ratio of the system from the photometric analysis alone (see Section 3), so estimated it using several independent methods.The first method was the TOD-COR light ratio reported in Table 2, the second was the optimal fitting of the disentangled spectra (ODS), and a third method is developed here.
We minimize the sum of the squared residuals between synthetic composite spectra and the observed spectra, where the synthetic composite spectra are built by adding Doppler-shifted synthetic spectra weighted by light fractions according to some value of ℓ B /ℓ A , and with atmospheric parameters from the ODS analysis.We search in a grid of ℓ B /ℓ A values between 0.8 and 2 with ten uniformly spaced samples; these bounds were determined after an initial trial run with bounds of 0.5 to 3. The minimum of a polynomial fit to the sum of the square residuals against trial values for ℓ B /ℓ A then yields a best estimate for its value.The only free parameters of the minimization for each trial ℓ B /ℓ A were the coefficients of a sixth-order polynomial that was used to normalize the observed spectra against the synthetic spectrum.The applied Doppler shifts were fixed according to the corresponding RV values derived in Section 4. The minimization was carried out using the Nelder-Mead method as implemented in the Scipy python package MINIMIZE (Virtanen et al. 2020).
We repeated the process for all spectral orders showing sufficient well-resolved lines in regions unaffected by tellurics from Earth's atmosphere.To save computing time, only the two observations closest to positions of quadrature were used.This approach resulted in  30 spectral orders showing adequate fits with well defined minima in the sum of the square residuals.Fig. 7 shows the result of this method applied to order 66, which demonstrates the effectiveness of optimizing the normalization of the observed spectra in the fitting routine for an order where line blending is significant.The resulting values for ℓ B /ℓ A after applying this method to the selected orders are shown in Fig. 8, where each result is the average of the result from each of the two observations used.Thus, after taking the average of the results in Fig. 8, which corresponds to the blue line, we have three independent estimations for the light ratio of the system, i.e., TODCOR, ODS, and direct fitting for the light ratio.These estimations are given in Table 4, where we take the average and its standard deviation from the ODS method.For the rest of this study, we adopt ℓ B /ℓ A = 1.300 as the weighted mean of those values, with an errorbar of ± 0.036 which is their standard deviation.This value is used to ascertain the correct light curve solution is obtained in the next section.

ANALYSIS OF THE LIGHT CURVE
The components of KIC 9841944 are close to each other and thus have a significant tidal deformation.We therefore sought to model the light curve using a code that is based on Roche geometry.We selected the Wilson-Devinney (WD) code (Wilson & Devinney 1971;Wilson 1979) for this, and used the 2004 version of the code driven using the JKTWD wrapper (Southworth et al. 2011).The user guide which accompanies the WD code (Wilson & Van Hamme 2004) includes a description of all input and output quantities discussed below.
The WD code is computationally expensive and is not suited to the analysis of the full 500 000 short-cadence datapoints in one step.We therefore used the orbital ephemeris determined in Section 3 to convert the datapoints to orbital phase, and then binned them into a much smaller number of points.We chose a bin size of 0.001 orbital phases during the eclipses and 0.005 outside the eclipses, resulting in a total of 456 phase-binned datapoints suitable for analysis with the WD code.
Through a process of trying a large number of different fits with a range of fixed and fitted parameters, we arrived at a good solution to the light curve.We adopt this as the default solution, plot it in Fig. 9, and give the fitted parameters in Table 6.It was obtained in Mode = 0 with the following fitted parameters: the light contributions of the two stars; their potentials; their gravity darkening coefficients; the linear coefficients of the logarithmic LD law; and the orbital inclination.We fixed the mass ratio at the spectroscopic value, the orbital eccentricity to zero, the albedos to 1.1, the logarithmic coefficients of LD to values from Van Hamme (1993), third light to zero, and the rotation rates to synchronous.We adopted the maximum numerical precision values of N1 = N2 = 60 and N1L = N2L = 60, the simple treatment of reflection, and the Cousins R-band as a proxy for the Kepler passband.This solution gives a light ratio in excellent agreement with the spectroscopic values in Table 4.The fractional radii in Table 6 are volume-equivalent values calculated by the LC component of the WD code.
The main confounding factors in fitting the light curve were the albedo and gravity darkening.Albedo values around 1.0 give a good fit to the data, but there is a wide local minimum of χ 2 at albedos in the region of 0.0 that gives a worse fit but is frequently found by the steepest-descent minimization method implemented in the WD code.After extensive experimentation, we found that fixing the albedos to values slightly above 1.0 yielded the best fits and did not risk descending into the local minimum at lower abedo values.Fitting for gravity darkening turned out to be a crucial step in obtaining a good model, and the best-fitting values were somewhat varied but generally around 0.8 for the primary star and 1.1 for the secondary star.We interpret this as a systematic issue caused by residual pulsations in the light curve and account for it in the uncertainties.
For the determination of the uncertainties of the fitted parameters we ran a wide range of solutions with a large number of different possible approaches to modelling the light curve (see Southworth 2020).We calculated the errorbar for each parameter by adding in quadrature the contribution from every model choice, which in turn was taken to be the amount that parameter changed by versus the default solution.The following different approaches were explored.
(i) We varied the amount of binning prior to the WD solution, finding that it had a negligible effect on the fitted parameters.
(iii) We tried the detailed reflection effect and found almost identical results.
(iv) Attempts to fit for mass ratio returned a value very close to the spectroscopic one and almost no change in the other parameters.
(v) Fitting for the rotation rates of the stars yielded values close to and consistent with synchronous rotation and little change in the other parameters.
(vi) Fitting for albedo gave a poorer fit but again very little change in the other parameters.
(vii) Allowing for third light gave very similar parameters and a very small amount of third light consistent with zero.
(viii) Fitting for one of the temperatures of the stars instead of the two light contributions directly (namely Mode = 2) changed the fractional radii by 0.2%.
(ix) Fixing the limb darkening coefficients to theoretical values changed the fractional radii by 0.5%.
(x) Using the square-root instead of logarithmic limb darkening law had a negligible effect.
The best-fitting parameters were highly robust against all these experiments.The only significantly discrepant fit (neglecting our original exploratory ones) was when we used the Cousins I-band instead of the R-band.However, we are able to reject this fit as it is not consistent with the spectroscopic light ratio.The calculated uncertainties are given in Table 6.
Fig. 9 shows that the best fit is extremely good, with an r.m.s. of 0.20 mmag, but that there are systematics remaining in the residuals.We attribute these to the WD numerical integration for points during the eclipses, and residual pulsations for points outside eclipse.Our approach of phase-binning the data gives a light curve practially without Poisson noise, so makes any imperfections in the fit easily noticable.We note that the systematics in the residuals of the fits we found to the data are too small to show up in plots of unbinned data, so may well be present in previous work on this object.

PHYSICAL PROPERTIES
We determined the physical properties of KIC 9851944 from the spectroscopic and photometric results obtained above.For this we used the JKTABSDIM code (Southworth et al. 2005), modified to use the IAU system of nominal solar values (Prša et al. 2016) plus the NIST 2018 values for the Newtonian gravitational constant and the Stefan-Boltzmann constant.Errorbars were propagated via a perturbation analysis.The results are given in Table .7.
We determined the distance to the system using optical BV magnitudes from APASS (Henden et al. 2012), near-IR JHK s magnitudes from 2MASS (Cutri et al. 2003) converted to the Johnson system using the transformations from Carpenter (2001), and surface brightness relations from Kervella et al. (2004).The interstellar reddening was determined by requiring the optical and near-IR distances to match.We found a final distance of 935 ± 12 pc, which is significantly shorter than the distance of 999 ± 12 pc from the Gaia DR3 Table 5. Summary of the parameters for the WD2004 solutions of the light curves of the systems.Detailed descriptions of the control parameters can be found in the WD code user guide (Wilson & Van Hamme 2004).A and B refer to the primary and secondary stars, respectively.Uncertainties are only quoted when they have been robustly assessed by comparison of a full set of alternative solutions.parallax (Gaia Collaboration 2016, 2021).We have no explanation for this at present, but note that the Gaia DR2 and DR3 parallaxes of this object differ by nearly 2σ so might be affected by its binarity.

Frequency analysis
Following the binary modelling, we continue with the asteroseismic analysis of the target.Because the observed pulsations have much smaller amplitudes than the binary signal, the quality of the TESS and WASP data are insufficient for the asteroseismic analysis, and we limit ourselves to using the residual Kepler light curve.This is the merged light curve of all available Kepler short-cadence data after subtracting the best-fitting binary model, hereafter referred to as the pulsation light curve.To minimize the impact of outliers and instrumental effects on the asteroseismic analysis of small-amplitude pulsations, we apply additional processing to the data.Firstly, we remove those parts where coronal mass ejections (CMEs) or thermal and pointing changes of the spacecraft, such as at the start of a quarter or after a safe-mode event, have a visible impact on the quality of the light curve.Secondly, we apply preliminary iterative pre-whitening (as described by, e.g., Van Reeth et al. 2023) to build a tentative mathematical model of the 20 most dominant pulsations using a sum of sine waves where a i , f i and φ i are the amplitude, frequency, and phase of the i th sine wave, respectively, and t 0 is the average time stamp of all data points in the pulsation light curve.We then identify individual outliers in the residuals using 5σ clipping, and remove these data points from the pulsation light curve as well.
Next, we use iterative pre-whitening to measure the pulsation frequencies from the resulting cleaned pulsation light curve.Hereby we iteratively fit additional sine waves to the time series, with frequencies that correspond to the dominant amplitudes in the Lomb-Scargle periodogram (Scargle 1982) of the residual pulsation light curve.However, in the particular case of KIC 9851944, there is a significant amount of red noise present in the data, most likely caused by residual instrumental effects, even after cleaning the pulsation light curve.To ensure that we detect all significant pulsation frequencies (with a signal-to-noise ratio S/N ≥ 5.6; see e.g.Baran et al. 2015), we customise our approach.Firstly, we split the frequency range between 0 d −1 and the Nyquist frequency (734.21d −1 ) in overlapping parts and apply the iterative pre-whitening to these individually: from 0 d −1 to 2 d −1 , from 1 d −1 to 6 d −1 , from 4 d −1 to 11 d −1 , from 9 d −1 to 21 d −1 , from 19 d −1 to 51 d −1 , from 49 d −1 to 201 d −1 , and from 199 d −1 to the Nyquist frequency.In this stage, we measure all frequencies with S/N ≥ 4.0.Secondly, we merge the different frequency lists, keeping only those frequencies that are dominant in a 2.5 f res -window, where f res is the frequency resolution of the light curve (Loumos & Deeming 1978) with a value of 0.00208 d −1 , and that have S/N ≥ 5.6.Thirdly, we non-linearly optimize this filtered frequency list using the least-squares minimization with the trust region reflective method from the lmfit python package (Newville et al. 2019).After the optimization, we redetermined the S/N that are associated with the different frequencies, again only keeping those with S/N ≥ 5.6.This leaves us with a final list of 133 measured frequencies in both the p-and g-mode regimes, considerably more than the 89 frequencies reported by Guo et al. (2016).Notably, we formally detect one frequency at 391.5095±0.0002d −1 with an amplitude of 6.3±1.2 µmag.However, a detailed analysis of this sine wave reveals that it originates from the noise properties of the short-cadence light curve.The signal has a maximal amplitude during the first parts of quarters 13 and 14, but is not detectable in most other parts of the light curve.Hence, we discard this frequency and limit ourselves to the remaining frequencies, with values below 25 d −1 .These are illustrated in Fig. 10 and listed

Tidal perturbation analysis
As already demonstrated by Guo et al. (2016) and illustrated in Fig. 11, many of the detected frequencies form orbital-frequency spaced multiplets, also called "tidally split multiplets".In our work, we identify each multiplet by looping over all measured frequencies f i in order of decreasing amplitude, and consider frequencies f j to be part of a multiplet around it if where n is an integer chosen to minimise the left-hand side of the Equation, and ∆ f is a chosen frequency tolerance of 0.002 d −1 , that is ≈ f res .In each iteration, the considered frequencies f j have amplitudes a j that are smaller than the amplitude a i , associated with f i , and are not yet associated with a different multiplet.
We identify 13 multiplets that consist of three or more com-ponents, as listed in Table B in Appendix B. One of these, with the dominant frequency at 1.386342(4) d −1 , consists of 19 frequencies that match integer multiples of the binary orbital frequency f orb within ∆ f , and is discussed in detail below in Sect.8.3.In their work, Guo et al. (2016) explained the other multiplets as consequences of (i) rotational frequency splitting, and (ii) partial occultations of the pulsation mode geometries during the eclipses (e.g., Reed et al. 2001Reed et al. , 2005;;Gamarova et al. 2003;Rodríguez et al. 2004;Gamarova et al. 2005).However, it is now known that pulsation modes are often tidally perturbed (e.g., Samadi Ghadim et al. 2018;Steindl et al. 2021;Van Reeth et al. 2023) or tilted (e.g., Fuller et al. 2020;Handler et al. 2020;Kurtz et al. 2020;Van Reeth et al. 2022) in close binaries.Hence, to determine the true origin of the detected multiplets, we evaluate their dependence on binary orbital phase in detail.Each multiplet with three or more components is studied individually by removing all measured variability that is not associated with that multiplet from the pulsation light curve, that is, we subtract all sine waves with frequencies that are not in the multiplet from the pulsation light curve.We then fold the residual light curve with the orbital period, split the data in 50 bins, and fit the dominant sine wave of the multiplet to the data, optimising the amplitude and phase for each orbital-phase bin.The reason for this approach is twofold: (i) the frequency spacings within the detected multiplets are not always exact multiples of f orb , and (ii) there can be frequencies missing within the detected multiplets.As a result, the analytical reconstruction from Jayaraman et al. ( 2022) is not well suited for this star.
The results are illustrated in Fig. 12 for the dominant pulsation, and shown in Appendix C in Figs.C1 to C12 for all pulsations.In each figure, the relevant frequency multiplet is plotted in the lefthand panel, with a white star marking the dominant frequency.On the right-hand side, the middle and bottom panels show the orbitphase dependence of the pulsation amplitude and phase, respectively.For reference, the top right panel shows the orbital-phasefolded light curve of the binary.From these figures we can draw several conclusions.Firstly, because most pulsation amplitudes and phases vary significantly at all orbital phases, and the scales of the observed pulsation phase modulations are of the order of 0.5 to 1.5 rad, we can conclude that the observed pulsations are tidally perturbed.While tidally tilted pulsations are expected to have pul-sation phase modulations of 0 rad, π rad or 2π rad, which can be smeared out (Fuller et al. 2020), tidal perturbations can result in much smaller pulsation phase modulations (e.g., Van Reeth et al. 2023).Moreover, the observed tidally split multiplets are not perfectly equidistant.This is in contradiction with current theoretical predictions made for tidally perturbed (e.g., Smeyers 2005) and tidally tilted pulsations (Fuller et al. 2020), indicating that aspects that are currently not included in these theoretical frameworks, such as the Coriolis force, also play a role.Secondly, based on the different morphology of the curves for the different pulsations, we can conclude that the pulsations have different mode geometries (Van Reeth et al. 2023).For example, while most observed pulsations are modulated twice per orbit, the amplitude of the pulsation with frequency f = 5.097165(3) d −1 only reaches one maximum per orbital cycle.Moreover, for some pulsations the pulsation phase decreases as a function of the orbital phase when the observed amplitudes are maximal, while for others the pulsation phase increases.There is no detectable correlation between these effects and the pulsation frequency, but the tidal modulations are most easily observed for pulsations with higher (average) amplitudes.This suggests that most if not all pulsations in this target are tidally perturbed, but that the S/N of the lower-amplitude pulsations is too low for a detection.
Finally, we can confirm that some pulsations belong to the primary, while others could belong to either the primary or the secondary component, in agreement with the inferences made by Guo et al. (2016) based on asteroseismic models.As seen in Fig. 12, the observable amplitude of the p mode with frequency 10.39971 d −1 drops during the primary eclipse, indicating that it belongs to the primary.The higher observed amplitudes just before and after the primary eclipse signifies that the pulsation has a higher amplitude on the side of the primary that is facing the secondary component, similar to what has been observed for g-mode pulsations in V456 Cyg (Van Reeth et al. 2022), andKIC 3228863 andKIC 12785282 (Van Reeth et al. 2023).By contrast, the observable amplitudes of the p modes with frequencies 11.52234 d −1 and 11.89048 d −1 peak during the primary eclipse, as shown in Figs.C6  and C7.This can either indicate that these pulsations belong to the secondary component, or it can be caused by reduced geometric mode cancellation during the eclipse, depending on the geometry of these pulsation modes.The origin of the observed pulsations can be investigated further using detailed asteroseismic modelling, which would allow us to calculate the probability that specific pulsations belong to one or the other component.However, such a modelling study lies outside the scope of the current work.

Orbital harmonic frequencies
In addition to the tidal perturbation of self-excited pulsations, and in agreement with Guo et al. (2016), we also detect an orbital harmonic frequency comb, illustrated in Fig. 13, which can indicate tidally excited oscillations.However, we do not recover the same orbital harmonics as Guo et al. (2016).They reported the detection of 8 f orb , 22 f orb , 46 f orb , and 50 f orb .While our detected orbital harmonic frequency comb consists of 19 frequencies, we only recover 22 f orb .This discrepancy is likely caused by differences between the reduced short-and long-cadence Kepler light curves, and between the photometric binary models.Thus, while we confirm the presence of orbital-phase dependent variability in the pulsation light curve, as seen in the bottom right panel of Fig. 13, its exact observed characteristics have to be treated with caution.
As pointed out by Guo et al. (2016), the detection of an orbital harmonic frequency comb is somewhat unexpected for synchronised binaries with circular orbits, though it has also been detected for other such systems (da Silva et al. 2014).Because the orbital eccentricity of the binary is zero, the equilibrium tides that are responsible for deforming the star and perturbing the pulsations are considerably larger than the dynamical tides that excite oscillations.Hence, this can indicate that this system has a slightly eccentric orbit or that one or both of the components is asynchronously rotating, within the uncertainty margins of our measurements.

Gravity-mode period-spacing pattern
Finally, we analyse the observed g-mode pulsations in detail.Li et al. (2020) reported the detection of a short period-spacing pattern of prograde sectoral quadrupole modes, that is with (k, m) = (0, 2), and assigned them to the primary component of the system, based on the stellar masses determined by Guo et al. (2016).In this work, we confirm the pattern detection by Li et al. (2020) using the methodology from Van Reeth et al. (2015), as shown in Fig. 14.Only the g mode with frequency ∼ 0.461 d −1 was not recovered because of the higher signal-to-noise ratio cutoff that we used, that is, S/N ≥ 5.6 instead of 4.0.Moreover, the dominant g mode in the pattern, with a frequency f of 2.239718(6) d −1 , was found to be tidally perturbed and exhibit spatial filtering during the primary eclipse, as shown in Fig. C1 and with the corresponding multiplet listed in Table B.These observations confirm that the g modes belong to the primary component of the system: the observed tidal perturbations exhibit a dip in the amplitude and a saw-tooth-like "glitch" in the pulsation phase during the primary eclipse, which can only be explained if the pulsation belongs to the primary (Van Reeth et al. 2022, 2023).
To investigate the potentially asynchronous rotation of the pulsating star, we modelled our detected period-spacing pattern by fitting an asymptotic period-spacing series, following the method described by Van Reeth et al. (2016).We confirmed the pulsation mode identification found by Li et al. (2020), (k, m) = (0, 2), and found that the star has a buoyancy radius Π 0 of 4370 +690 −660 s and a near-core rotation frequency f rot of 0.49 +0.05 −0.06 d −1 .These values agree within 1σ with the results from Li et al. (2020), who found Π 0 and f rot values of 3500 ± 500 s and 0.41 ± 0.05 d −1 , respectively.Moreover, both sets of values are consistent with synchronous rotation.However, because the detected pattern is so short, the uncertainties on these f rot values are relatively large.Hence, there is still a possibility that the primary is asynchronously rotating, but insufficiently strongly to be detected with the available data. the period spacing between consecutive modes in the detected pattern, as a function of the pulsation period.Because there is an undetected mode between the fourth and fifth detected pulsation modes, the fourth period spacing in the pattern is not shown.
The error margins are smaller than the symbol sizes.

DISCUSSION
We compared the observed properties of the components of KIC 9851944 to isochrones calculated by MIST (Mesa Isochrones and Stellar Tracks) using the MESA code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015;;Dotter 2016;Choi et al. 2016;Paxton et al. 2019).We searched for the best MIST isochrone using two methods.For method 1, we interpolated radius, T eff and luminosity as functions of mass.We then estimated these parameters using the interpolants at our measured masses.The objective function to minimize is then the sum of the quadrature distances between the interpolated and observed locations of the components with age and metallicity as free parameters.For method 2, we included mass in the calculation of the objective function, i.e., masses are not constrained to the observed values.For the minimization, we used SciPy's implementation of the stochastic differential evolution algorithm (Storn & Price 1997;Virtanen et al. 2020).Our grid of isochrones spanned from −4.0 to 0.5 dex in [Fe/H] and 0.5 to 10.3 dex in log 10 (age).The grid spacing in [Fe/H] was 0.05 dex between −0.5 and 0.5 dex and 0.25 dex outside that range up to ±2 dex, beyond which the grid spacing was 0.5 dex.A grid spacing of 0.05 dex in log 10 (age) was used through the full grid.
After the best matching coeval isochrone was found, we removed the corresponding metallicity value from our isochrone grid and repeated the procedure to explore the effect of varying [Fe/H] on the predicted age of the system.We did this twice, leading to first, second, and third best matching coeval isochrones, each corresponding to a different metallicity.For method 1, the metallicity values of the first, second and third best matching isochrones were 0.05, 0.1 and 0.0 dex, respectively.For method 2, these values correspond to 0.05, 0.0 and 0.1 dex, i.e., the same set, except the second and third best matches are swapped.In all cases, the predicted age of the system is 1.259 ± 0.073 Gyr.The error on this age estimate is taken to be half the average grid spacing either side of the best fitting isochrone.
Table 7. Model parameters of the best fitting isochrone with [Fe/H] = 0.05 dex and an age of 1.259 ± 0.073 Gyr.The results obtained from both comparison methods are given.Also given are the differences ∆ between the observations and model for each parameter; we quote this difference in units of the uncertainty associated to the observations in brackets.

Parameter
3.214 7.8 (5.3σ ) 3.096 3.8 (2.6σ ) We note that this age estimation is in excellent agreement with the estimation by Guo et al. (2016) of 1.25 Gyr.
The three best-fitting isochrones are shown in Fig. 15 with the observed locations of the objects over-plotted in the M − log(T eff ), M − log(R), M − log(L) and HR diagrams; also shown are their isochrone neighbours, i.e., ± the grid spacing in age.For both model comparison methods, the isochrone with [Fe/H] = 0.05 dex yielded the model with the smallest quadrature distance to the observed quantities that were included in the objective functions so we adopt the corresponding model predictions from now on.Those values and their differences relative to the observed values are given in Table 7 for both model comparison methods.The locations of these values in the planes of Fig. 15 are also shown; note the relative scales of the axes, particularly for T eff where the scale of the axis is much smaller than the scales of the other axes.
Table 7 confirms that the luminosity predictions from our best matching isochrone resulting from both comparison methods are accurate for both components and this is also clear in the lower, left panel of Fig. 15.The other parameters are also in good agreement for star A, particularly regarding those resulting from method 2 where the agreement is excellent.However, for star B, the radius is underestimated while the T eff is overestimated, and these discrepancies compensate for each other to yield the accurate luminosity prediction.
A possible explanation is that the secondary star is approaching the terminal age main sequence (TAMS), i.e., an evolved stage where the sensitivity of the models increases.Thus, a full evolutionary modelling analysis, such as that carried out by Guo et al. (2016), might yield better model predictions because metallicity and other parameters relating to, e.g., overshooting, are included as free parameters so the model is more flexible.This would also provide the means for a more detailed discussion of model comparison with observations but is beyond the scope of this paper.Here, we simply searched in grids of pre-computed models, which are limited by the size of the grid steps in metallicity and age, as well as fixed input physics, to estimate the age of the system.The parameter ranges that the isochrone neighbours span in the planes shown in Fig. 15 (i.e., about twice the uncertainty), and the agreement with the value determined by Guo et al. (2016), is evidence that this age estimation is accurate.
The evolutionary tracks corresponding to the masses of the models in Table 7 are shown in the T eff -R, and HR diagrams in Figs.16 and 17, respectively, as well as the observed locations of the components.While Guo et al. (2016) find the secondary to be in the hydrogen-shell burning phase for the same age estimate as ours, we find that the secondary has not yet exhausted the hydrogen in its core, and is approaching the TAMS.The scenario found here is more likely to be observed on a statistical basis since the evolution up the HR diagram after the blue loop occurs on a very short timescale.In any case, the secondary is more evolved and thus larger and more massive than the primary, but cooler.
Also shown in those figures are is the blue and red edges of the instability strips for low-order p-(dashed lines), and g-modes (dotted lines) in δ Scuti and γ Doradus stars, respectively, from Xiong et al. (2016).Both stars are well within both instability domains which is in agreement with the findings in Section 8 that both components might be pulsating.Guo et al. (2016) find that the secondary is slightly hotter than the blue edge of the γ Doradus instability domain calculated by Dupret et al. (2005).This reflects the fact that the calculations by Xiong et al. (2016) (used here) predict a much broader overlap between the δ Scuti and γ Doradus instability domains so stars with a larger range in T eff are expected to pulsate simultaneously in p-and g-modes, i.e., hybrids are expected to be more common.We also plotted a sub-solar ZAMS ([Fe/H] = -0.25 dex) as well as the ZAMS corresponding to our best matching isochrone ([Fe/H = 0.05 dex] in Fig. 16 and Fig. 17; this shows how the ZAMS is raised to higher luminosities and radii for higher metallicity. We present the results previously derived by Guo et al. (2016) in Table 8.The 13 spectroscopic observations that the previous authors used to derive their RVs and atmospheric parameters have resolutions of R = 6000.While it was already discussed in section 5.1 that our atmospheric parameters agree with those derived by Guo et al. (2016), the higher resolution (R = 60000) of the spectra, and larger number of observations (33) used in this work yielded a much higher quality RV curve.This is reflected by comparing the reported masses between the two studies; the precision in the mass estimates by Guo et al. (2016) are 4.0% and 3.9% for the primary and secondary, respectively, while we attain precisions of 0.57% and 0.59%.The two sets of mass measurements agree to within 0.2σ for M A and 1.1σ for M B .In contrast, our radius measurements differ from those of Guo et al. (2016) by 6.6σ for star A and 3.5σ for star B. However, Guo et al. (2016) note the discrepancy between their values for the radius ratio derived from two spectroscopic methods (k = 1.22 ± 0.05 and k = 1.27 ± 0.29) and the value derived from their light curve modelling (1.41 ± 0.018), where they tentatively adopt the radii associated to the latter.The radius ratio derived here agrees with the values derived by Guo et al. (2016) using their spectroscopic methods within 0.8 σ and 0.3 σ , but shows the same discrepancy compared to the value derived from their light curve modelling, for which Guo et al. (2016) note there exists a family of comparable solutions due to the partial nature of the eclipses.
We have assumed a circularised and synchronized orbit for KIC 9851944 because our RV-and light-curve solutions were consistent with a circular orbit, and the values for the component of the synchronous velocity along our line of site are in excellent agreement with the values for v sin i derived in section 5.1 (see Table 3).Further justification was provided by Guo et al. (2016) where they examined the eclipse times by Conroy et al. (2014) and Gies et al. (2015), finding that the median deviation of the phase difference between the primary and secondary eclipses from that of a circular orbit suggests e ≥ 0.0001.Furthermore, the circularization timescale for a binary system like KIC 9851944 is 600 Myr, and the synchronization timescale is an order of magnitude shorter (Zahn 1977;Khaliullin & Khaliullina 2010;Guo et al. 2016); these timescales are shorter than the age of the system reported by both studies.
Our estimations of the near-core and surface rotations reported in section 5.1 and 8 suggest KIC 9851944 is rotating rigidly and syn- chronously.This is similar to the findings by e.g.Guo et al. (2017); Guo & Li (2019), that the short-period EBs KIC 9592855 and KIC 7385478 both contain a γ Doradus pulsator that is tidally synchronised at the surface as well as at the near-core regions.Such measurements allow for a calibration of the time-scales for synchronization at the surface compared to the core, and thus the time-scales associated to angular momentum transport, which further aids in the discrimination among different theories.
In contrast to the above, our detection of the orbital harmonic frequency comb for KIC 9851944 either suggests asynchronous rotation or a non-zero eccentricity that has not been detected due to observational error because these can be indicative of tidally excited oscillations.It is unclear whether the very small value for the lower limit on the eccentricity, i.e., e > 0.0001, reported above for this system would be enough to induce tidal excitation of modes at the amplitudes observed here.
The pulsation analysis here complements the study by Guo et al. (2016) who report splittings to some of the pulsation modes.We confirm the detection of tidally split multiplets, and explain their origin; we present evidence to suggest that these are due to perturbations to the pulsation mode cavaties, i.e., tidally perturbed pulsations, by investigating their phase and amplitude dependencies with orbital phase.Guo et al. (2016) attempted to interpret the modes by comparing the observations with theoretically computed frequencies, from which they concluded that the observations can be explained by loworder p modes in the primary and the secondary, or the g mode and  7.The evolutionary tracks are shown in green and purple for star A and star B, respectively, and their observed locations are indicated by the blue and black markers.The blue (blue lines) and red (red lines) edges of the instability domains for low-order p-and g-modes calculated by Xiong et al. (2016) for δ Scuti (dashed lines) and γ Doradus (dotted lines) stars are indicated.The thin, black line is the ZAMS corresponding to the models in  mixed modes of the secondary (Guo et al. 2016).We confirm that some of the p modes belong to the primary, and others could belong to either the primary or the secondary from modulation of the amplitudes during eclipses.Our evidence suggests the primary is the hybrid, and this is because the saw-tooth-like "glitch" in pulsation phase of the tidally perturbed g mode (see Fig. C1) can only be explained if the pulsation belongs to the primary (Van Reeth et al. 2022, 2023).Guo et al. (2016) report that mode identification was inconclusive, reflecting the difficulty in identifying the p-modes in δ Scuti stars.

CONCLUSIONS
We have determined the physical properties of KIC 9851944, a short-period detached eclipsing binary containing two F-type stars, both of which pulsate.Our analysis is based on 33 échelle spectra plus light curves from the Kepler and TESS missions.We measure masses and T eff s to 0.6%, radii to 1.0% and 1.5%, and 133 frequencies due to p-and g-mode pulsations.We find no evidence of a third component, apsidal motion, or eccentricity.We estimate the age of the system to be ∼1.26Gyr by comparison of the measured properties to the MIST model isochrones.
We investigated the systemic errors associated to using the cross correlation technique to extract RVs, which arise due to blending between the spectral lines of the components in a binary system.We find that the effect is small when using TODCOR but correcting for the issue is still necessary because the resulting shifts are clearly systematic in nature and have a non-negligible impact on the results.We used three independent spectroscopic methods to determine the light ratio for the system.The results are in agreement with each other and consistent with the value obtained from modelling the light curve which supports the reliability of our light curve solution.We compared our results to those reported by Guo et al. (2016) and find that we have improved the precision of the measured masses significantly, but the precision in the radius estimates are comparable.Both these outcomes are expected since we use much higher-resolution spectra but the same photometry (Kepler).
By analysing the residuals of the light curve model, we confirm the detection of tidally perturbed p-mode pulsations, possibly in both components of KIC 9851944.A short period spacing pattern was detected among the g modes and was assigned to the primary component, where perturbations to one of the g-modes was detected.Thus, the primary component is a δ Scuti/γ Doradus hybrid.If pulsation mode identification can be performed (as in, e.g., Bedding et al. 2020) on the p modes, KIC 9851944 will become a well-equipped laboratory for stellar physics.lines with synthetics, and sharing his thoughts on the possibility that the components of KIC 9851944 are Am stars.
We are grateful to Kelsey Clubb for assisting in obtaining and reducing the spectroscopic observations from the Shane telescope.TVR gratefully acknowledges funding from the KU Leuven Research Council (grant C16/18/005: PARADISE).We gratefully acknowledge financial support from the Science and Technology Facilities Council.This research has made use of the SIMBAD and da Silva R., Maceroni C., Gandolfi D., Lehmann H., Hatzes A. P., 2014, A&A, 565, A55 de Mink S. E., Pols O. R., Hilditch R. W., 2007, A&A, 467, 1181 Table B1.Parameter values associated with the iteratively prewhitened frequencies, which were obtained as described in Sec.8.1.The frequencies are grouped, with the first set consisting of those labelled as independent frequencies, and the following groups consisting of combination frequencies, sorted according to the dominant parent frequency.In the last column, we list the identified combinations.Combination frequencies that are found within the frequency resolution f res but not within 3σ , are indicated with * in the first column and with ≈ in the final column.In the middle and bottom right panels, we see the orbital-phase dependent modulations of the pulsation amplitudes and phases, respectively.These are calculated in 50 data bins, with the 1σ uncertainty range indicated by the dashed lines.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Cuts of the Kepler Quarter 0 (top panel) and TESS Sector 15 (middle panel) processed light curves.Also shown in the bottom panels are the phasefolded light curves observed by WASP between May 2007 and July 2007 (left), and between June 2008 and July 2008 (right).

Figure 2 .
Figure 2. JKTEBOP model fit to the phase folded Kepler light curve from quarter 0.

Figure 3 .
Figure 3. O-C plot from the linear fit (red line) to the times of primary minimum (black dots), where the grey shaded region represents 1 σ uncertainty associated to the corresponding calculated values..The blue dashed line represents an attempt to additionally fit for a quadratic term.

Figure 4 .
Figure4.Top: orbital fit to the corrected RVs using the merged approach.Bottom: velocity corrections for the corresponding RVs.

Figure 5 .
Figure 5. Portion of disentangled spectra in the spectral range 4900-5290 Å.The similarity between the spectra of both components is obvious.

Figure 6 .
Figure6.Optimal fitting of the wings for Hβ lines in disentangled spectra of binary system KIC 9851944.The spectral lines of metals were masked during the fitting, so were not included in calculations of the merit function.The disentangled spectra are shown in black, and the optimal fits in red.

Figure 8 .
Figure 8. Resulting ℓ B /ℓ A values from the grid search method averaged over the orders.

Figure 9 .
Figure9.The best-fitting WD model (green line) to the Kepler phase-binned light curve of KIC 9851944 (red filled circles).The residuals of the fit are plotted in the lower panel using a greatly enlarged y-axis to bring out the detail.

Figure 11 .
Figure 11.Échelle diagram of the prewhitened frequencies of KIC 9851944, folded with the orbital frequency f orb .Identified multiplets are indicated in different colours ranging from purple to dark yellow, with a white star marking the dominant frequency of each multiplet.Single frequencies are shown in black.The marker sizes indicate the associated amplitudes.

Figure 12 .
Figure 12.Tidally perturbed pulsation with frequency f = 10.399706(2)d −1 .Left: Associated frequency multiplet, as shown in Fig. 11.The white star marks the dominant frequency of the multiplet.Top right: Orbital-phase folded light curve.Middle right: Orbital-phase dependent modulations of the pulsation amplitude, calculated in 50 data bins, with the 1σ uncertainty range indicated by the dashed lines.Bottom right: Orbital-phase dependent modulations of the pulsation phase, calculated in 50 data bins, with the 1σ uncertainty range indicated by the dashed lines.

Figure 13 .
Figure 13.Tidally excited oscillations in KIC 9851944.Left: Associated frequency multiplet, as shown in Fig. 11.The white star marks the dominant frequency of the multiplet.Top right: Orbital-phase folded light curve.Bottom right: Orbital-phase folded residuals of the light curve after subtraction of the binary model.Individual data points are shown in black.The overplotted purple line shows the average variability as a function of orbital phase, evaluated in 50 equal bins.

Figure 14 .
Figure14.Detected period-spacing pattern of g modes with (k,m) = (0,2) that belong to the primary component of KIC 9851944.Top: part of the Lomb-Scargle periodogram of the pulsation light curve (black) with the pulsation periods of the modes that form the pattern (red dashed lines).Bottom: the period spacing between consecutive modes in the detected pattern, as a function of the pulsation period.Because there is an undetected mode between the fourth and fifth detected pulsation modes, the fourth period spacing in the pattern is not shown.

Figure 15 .
Figure 15.The three best fitting isochrones with the observed locations of the objects over-plotted in the M − log(T eff ) (top-left), M − log(R) (top-right), M − log(L) (bottom-left) and HR (bottom-right) planes.Also plotted are the best fitting isochrone neighbours, i.e., ± the grid spacing in age.The locations of the isochrones within those planes corresponding to the model with the smallest quadrature distance to the measured quantities that were included in the objective function are indicated for those from both model comparison methods.

Figure 16 .
Figure16.The T eff -R plane showing the evolutionary tracks corresponding to the models presented in Table7.The evolutionary tracks are shown in green and purple for star A and star B, respectively, and their observed locations are indicated by the blue and black markers.The blue (blue lines) and red (red lines) edges of the instability domains for low-order p-and g-modes calculated byXiong et al. (2016) for δ Scuti (dashed lines) and γ Doradus (dotted lines) stars are indicated.The thin, black line is the ZAMS corresponding to the models in Table7 ([Fe/H] = 0.05 dex) and the dotted black line is the ZAMS for a metallicity of [Fe/H] = −0.25 dex.The thick, black line represents the best fitting isochrone which is shown as a solid grey line in Fig.15.Transparent grey dotted lines show the evolutionary tracks of [Fe/H] = 0.05 dex stars with other labelled masses.

Figure 17 .
Figure 17.Same as Fig.16 but in the Hertzsprung-Russell diagram.
APPENDIX C: DETECTED TIDALLY PERTURBED PULSATIONSIn this Appendix we illustrate all detected tidally perturbed pulsations of KIC 9851944.In each figure, the left panel shows a zoom-in of the associated frequency multiplet, as plotted in Fig.11.In the top right panel, we see the orbital-phase folded light curve of the binary.

Table 1 .
RVs and S/N corresponding to the spectroscopic observations taken at times given in the BJD column.An asterisk next to the BJD value means that observation was not used to derive the orbital parameters in Section 4.

Table 3 .
Determination of the atmospheric parameters from disentangled spectra for the components of KIC 9851944.The surface gravity for each component was fixed to the values determined from light curve and RVs solution as listed in Table7.

Table 4 .
Light ratios measured and adopted for KIC 9851944.

Table 6 .
(Prša et al. 2016es measured for the four systems analysed in this work.The units labelled with a 'N' are given in terms of the nominal solar quantities defined in IAU 2015 Resolution B3(Prša et al. 2016).
in Table B in Appendix B.

Table 7 (
[Fe/H] = 0.05 dex) and the dotted black line is the ZAMS for a metallicity of [Fe/H] = −0.25 dex.The thick, black line represents the best fitting isochrone which is shown as a solid grey line in Fig.15.Transparent grey dotted lines show the evolutionary tracks of [Fe/H] = 0.05 dex stars with other labelled masses.

Table 8 .
Guo et al. (2016)ed results for KIC 9851944 byGuo et al. (2016).Also given are the discrepancies ∆ of our results compared to those results, given as a percentage of the previous result as well as units of sigma σ .