© The Authors 2024

1. Introduction

Eclipsing binaries (EBs) are a sub-class of stellar binaries for which the orbital plane is along the observer’s line of sight. They are crucial to the field of stellar astrophysics as they provide a myriad of information. Eclipses in their light curves (LCs) are affected by the stellar radii, the orbital inclination, the luminosities of the stars, and several other factors. Modelling these LCs can thus be seen as an inverse problem of finding the set of properties for these stars that reproduces the observed phenomena.

However, LCs give us little to no information about the stellar masses. The mass of a star is arguably the most critical factor affecting its evolution. Double-lined EBs for which both components are visible in the observed spectra, provide us with the radial velocity (RV) semi-amplitudes, and hence the mass ratio and orbital elements, like the eccentricity and orbital separation between the stars. Further detailed analyses of the spectra can also anchor T_eff, log g, metallicity, and individual chemical abundances.

Detached EBs (DEBs) are such systems with well-separated stars, which have not undergone any interaction in the form of mass transfer. The stars in such systems are thus assumed to have evolved independently. Such systems are very important to paint an accurate picture of the age and evolution of stars by comparison with the theoretical models. They are also the gold standard when it comes to determining absolute stellar parameters as the accuracy of the derived parameters can be better than 1% (Gaulme & Guzik 2019; Kahraman Aliçavuş et al. 2023b).

While EBs provide us with model-independent fundamental and atmospheric stellar parameters, our understanding of the structure and processes of the stellar interiors mostly remain model-dependent. However, the phenomenon of stellar pulsations has an intrinsic driving mechanism. Accurate modelling of observed pulsations thus provides us with a tool to probe the internal structure of stars. Pulsating stars in EBs come closest to offering both in a single system. A significant fraction of pulsating stars are being discovered in EBs as a consequence of the availability of high-precision photometry from space-based telescopes like Kepler (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015).

δ Scuti stars are A to F-type main-sequence (MS), pre-MS, or early post-MS stars, typically with masses ranging from 1.5 to 2.5 M_⊙, although in certain cases this lower limit can extend down to 1.3 M_⊙ (McNamara 2011). They generally exhibit pressure mode oscillations with periods between 18 min and 8 h and amplitudes below $0\stackrel{\text{m}}{.}1$ in the V-band (Aerts et al. 2010). Some of these systems have been studied before to understand the effect of binarity on stellar pulsations. Some exhibit tidally induced and tidally tilted pulsations and help us gauge the effect of tidal interaction (Fuller et al. 2020). However, the list of such systems with comprehensive analysis remains short.

In this study, we present a detailed analysis of five EBs with candidates for δ Scuti-type pulsators as one of the components. We explain the process of sample selection and obtaining the data. Analyses for this study can be broadly divided into the following sections: photometric analysis, RV and spectroscopic analysis, pulsation analysis, and evolutionary modelling. We then discuss the results and present our conclusions and summary.

2. Target selection and observations

The targets for this project were selected in various ways. First, by using data obtained from the catalogue of detached EBs generated by the CRÉME survey (Ratajczak et al. 2013; Hełminiak et al. 2014, 2015, 2019b, 2021, 2022). The survey meticulously monitored a vast sample of more than 300 EBs and acquired multi-epoch spectra to enhance our understanding of these systems. Using a pre-cut on the projected minimum masses, Msin³(i), acquired from RV modelling, systems were selected in the mass range corresponding to that of δ Scuti-type pulsators. The LCs of these EBs were then visually inspected to look for pulsations. Three such systems are presented here: TICs 81702112 (HD 97329), 189784898 (V Cir), and 165459595 (V1109 Cen).

A similar approach was applied to targets from the catalogue DEBCat (Southworth 2015), with the focus on targets that only had ground-based photometric data published. In this work, we analyse one such system – TIC 308953703 (HIP 7666). Finally, one of us (MP) recently reported eclipses and pulsations discovered in the LC of TIC 386622782, a previously known spectroscopic double-lined (SB2) binary HR 2214. After this announcement, the target was included in the CRÉME sample, and scheduled for additional spectroscopic observations.

2.1. Photometry

High-precision TESS photometry, in the form of 2-min cadence LCs, is available for all the targets in the sample¹. These LCs are available through the Mikulski Archive for Space Telescopes (MAST) and were retrieved using the LIGHTKURVE package (Lightkurve Collaboration 2018). This short-cadence photometry ensures a good resolution of the pulsation frequencies, especially beneficial for short-period pulsators like δ Scuti-type variables, and ensures a reliable model to estimate stellar masses and radii. We chose the simple aperture photometry (SAP) fluxes for all the targets to avoid distortions in the LCs introduced through Pre-search Data Conditioning (PDC)-SAP fluxes. The data was converted to differential magnitude and shifted to a baseline of 10 to facilitate the input file requirements of the LC modelling code. The choice of this baseline was arbitrary and does not affect the physical parameters in the solution. The LCs of all five targets used in this study are shown in Fig. 1

Fig. 1. Short (2-min) cadence LCs of the targets, each for a single TESS sector.

2.2. High resolution spectra

Most of the spectra used for analysis were obtained using the CHIRON spectrograph (Tokovinin et al. 2013), attached to the 1.5 m SMARTS telescope installed at Cerro Tololo Inter-American Observatory (CTIO). Aiming for higher efficiency, the spectrograph was used in the fiber mode, providing a resolution of ∼28 000. Extracted and wavelength-calibrated spectra were obtained with the pipeline developed at Yale University (Tokovinin et al. 2013) and provided to the user. However, barycentric velocity corrections were done in-house using the bcvcor procedure within IRAF (Tody 1986). Forty echelle orders were combined and continuum corrected to be used for RV and spectral analysis.

CHIRON was used to observe TICs 81702112, 189784898, 165459595, and 386622782. For the latter, our new data were complemented with the archival measurements from Abt & Levy (1985). We would like to point out that our observing programme on CHIRON was stopped abruptly when the Association of Universities for Research in Astronomy announced the shutdown of the SMARTS telescopes in September 2023. We would have gathered and presented more data, especially for TIC 16545959, if not for this unfortunate event.

For TIC 189784898, we also obtained six spectra taken with the FEROS spectrograph at the MPG-2.2 m telescope (Kaufer et al. 1999), providing a higher resolution of 48 000. These data were taken in 2004 during an engineering run, and are available from the ESO Archive. No new spectroscopic data was obtained for the target TIC 308953703. The RV measurements obtained by Feng et al. (2021) were adopted while carrying out the RV+LC analysis. A short summary of the spectroscopic observations is provided in Table 1.

3. Spectroscopy

3.1. Radial velocities

To calculate the RV values, we used the two-dimensional cross-correlation technique implemented in the TODCOR program (Zucker & Mazeh 1994). Synthetic spectra used as templates were calculated using the ATLAS9 model atmosphere (Kurucz 1979). The values of T_eff, metallicity, and vsin(i) were taken as 6800 K, 0.0, and 30 km s⁻¹, respectively. Errors were calculated using a bootstrap approach (Hełminiak et al. 2012).

Orbital solutions for the extracted RVs were calculated using V2FIT (Konacki et al. 2010). The routine estimates the orbital parameters of a double-Keplerian orbit fitted to the data, and is minimised using a Levenberg-Marquardt scheme. In addition to this, simultaneous LC-RV modelling was performed for four of the five targets. The exception was TIC 81702112, due to high root mean square (rms) values of the secondary RV fluctuations leading to erroneous outcomes.

We fit for the orbital period, P_orb, zero-phase, T_p, systemic velocity, γ, velocity semi-amplitudes, K_1, 2, eccentricity, e and longitude at periastron passage, ω. The orbital fits are displayed in Fig. 2.

Fig. 2. Orbital fits obtained for the RV curves for the targets. Here, primary refers to the more massive star and secondary is the less massive companion. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595.

3.2. Broadening functions

In the context of LC and spectral analyses, similar to any multi-parameter problem, it is advisable to obtain and fix as many parameters from various methods as possible to mitigate degeneracy. In the context of spectral disentangling and fitting, we resorted to broadening functions (BFs) to mainly obtain two necessary parameters: the rotational velocities of the stars and their light contributions.

The BF depicts the spectral profiles in velocity space, containing the signatures of RV shifts across various lines, along with intrinsic stellar phenomena such as rotational broadening, spots, and pulsations (Rucinski 1999). The implementation of this method is described in detail by Moharana et al. (2023). The fit was obtained for all the spectra for which the peaks of both components were well separated, avoiding in- or near-eclipse epochs. The final parameter values were calculated by taking the average of the ones from individual fits. The line profiles are extremely noisy in the case of TIC 165459595, preventing us from performing BF analysis for this target. The obtained parameters are stated in Table 2.

3.3. Spectral disentangling

A lot can be inferred about a star by studying its spectra. Estimates of temperature and surface gravity obtained from spectral analysis can provide useful checks on the parameters, such as flux ratios and relative radii, obtained from the LC modelling. They also provide an estimate of the global metallicity, which is an important parameter when constraining age and evolution. However, spectra obtained for binaries need to be disentangled into separate components before we can proceed with the modelling.

One way to do this involves obtaining spectra at different epochs in the orbit. Due to the different Doppler shifts, the spectrum corresponding to each star will have a different shift in the observed composite spectrum. By using a number of such spectra and the associated RVs, it is possible to reconstruct an individual component spectrum independent of any template.

One implementation of this technique to disentangle stellar spectra is shift-and-add (Shenar et al. 2020; Quintero et al. 2020; González & Levato 2006). It is an iterative process whereby the component spectra are shifted according to the corresponding RVs and added at each iteration. This process uses a set of N spectra, each at a different epoch in the orbit. The algorithm begins by correcting the spectra of component A for the RVs and adds them up. This serves as a first approximation of the component A spectra. This co-added spectra is then subtracted from all the input spectra to obtain the first approximation of component B. This process is repeated for a large number of iterations in order to refine and improve the accuracy of the separation process. This process also improves the signal-to-noise (S/N) by the square root of N. Figure 3 illustrates the spectral disentangling process for TIC 81702112, showing one of the input spectra alongside the resulting component spectra obtained after applying the shift-and-add algorithm.

Fig. 3. Input and output spectra utilised and obtained during the spectral disentangling for the target TIC 81702112. Black: one of the observed spectra that was used as an input for the disentangling algorithm. This is a composite of contributions from both the stars. Blue: disentangled spectrum of the primary component. Red: disentangled spectrum of the secondary component.

To implement this algorithm, we used the code DISENTANGLING_SHIFT_AND_ADD (Shenar et al. 2020, 2022). We executed the procedure for a wavelength range of 500–580 nm to avoid broad-wing features while ensuring that sufficient narrow lines with decent S/N were present. We used the RV semi-amplitudes obtained using V2FIT. The flux ratio calculated using BFs was provided as an input, which was used by the routine to scale the disentangled components. For TIC 165459595, we used the flux ratio obtained from the JKTEBOP solution. In case of TIC 189784898, we only used the six FEROS (R = 48 000) spectra obtained for this target. Due to the high S/N of the input spectra and higher spectral resolution, the individual decomposed spectra were of better quality than for the other cases. The procedure was run for 30 000 iterations in each case.

The initial detrended spectrum has trends in the continuum as a result of the imperfections in the continuum normalisation of composite spectra, line profile variations etc. (Hensberge et al. 2008). These trends or biases are additive and can be removed by fitting the continuum and then subtracting it from the disentangled spectra. This gave us the final spectra that we used for spectral analysis.

3.4. Spectral analysis

We used the ISPEC framework (Blanco-Cuaresma et al. 2014a; Blanco-Cuaresma 2019) for the purpose of atmospheric parameter determination using the disentangled spectra. Further continuum normalisation within ISPEC for example, using third-order splines was avoided to ensure that line depths remained least affected. Shifts in the RVs resulting from the disentangling process were removed before analysis.

Within the synthetic spectral-fitting technique, ISPEC provides users with the ability to generate synthetic spectra on the go using several choices of radiative transfer codes. We used the code SPECTRUM (Gray & Corbally 1994) to generate synthetic spectra in the wavelength region of our interest, owing to its speed and accuracy. We used the line lists from the Gaia-ESO Survey (GES; Heiter et al. 2015), covering the wavelength range of 420-920 nm. ISPEC uses χ² minimisation to choose the best-fitting spectrum from the ones synthesised.

During the fitting process, we fixed the spectral resolution to the known value, and the values of surface gravity (log(g)) from the LC solution to minimise the degeneracy in fitted values. We performed the fitting procedure in two steps. Step one used the GES line-list² well suited to determine abundances in order to determine the global metallicity. In a trial run, we kept the vsin(i) free alongside other parameters to check if it roughly matched the value obtained from BF analysis. We then fixed it and estimated the metallicity. In this step, we excluded the line regions with broad-wing features (e.g. the Mg triplet), to make sure that fitted parameters were not affected by incorrect continuum normalisation. In the second step, we fixed the metallicity to the value obtained from step one, and using the line-list best suited for determining atmospheric parameters we fitted for T_eff, vsin(i) and v_mic. Macroturbulence velocity (v_mac) was determined using an empirical relation as in Blanco-Cuaresma et al. (2014b), which is intended to be used for FGK-type stars (∼4000–7000 K). Therefore, in the case of TIC 386622782 and TIC 165459595, we set it to zero to avoid degeneracy with vsin(i). The final parameter estimates obtained from this analysis are stated in Table 3, and the best-fit spectra are shown in Fig. 4.

4. Light Curve modelling

An expeditious look at the TESS LCs of our sample shows that they consist of well-separated stars. For the process of modelling the LCs, we needed something that was efficient but also robust and well-tested. We used version 40 of the JKTEBOP code, (Southworth et al. 2004a,b; Southworth 2013) which itself is based on the EBOP program (Popper & Etzel 1981; Etzel 1981). Using the provided parameters, it computes the LCs through integration of concentric circles that approximate the surface of each star. It accounts for the oblateness of the stars by modelling them as bi-axial ellipsoids, and uses their projections to calculate proximity effects, which is a sufficient approximation for oblateness levels below ∼4%.

JKTEBOP adopts the relative radii r_A and r_B, stellar radii relative to the orbital separation, as parameters governing the shape of the eclipses. This is implemented using two input values, r_A/r_B and r_A+r_B. We set both as free parameters. Other parameters set free during the fitting process are the inclination (i), surface brightness ratio (J), eccentricity (e), and argument of periastron (ω). This code also provides an option to fit for a third light (l₃). In our test runs, we found the values of l₃ to be either very close to zero or nonphysical in nature. Scatter arising from the pulsations is an added reason why we simply fixed the l₃ parameter to zero (in all cases except for TIC 165459595 due to its proximity to two other stars). The coefficients of gravitational and limb darkening were set values sourced from Claret (2017). This version of the code also allows the user to fit polynomials and (up to 9) sines to account for additional trends in the data affecting the binary LC. We used polynomial functions to account for some long-term trends and also sine functions in some cases to model possible signatures from a stellar spot.

For cases in which the RVs are available, we used them within the code to perform a simultaneous LC-RV fit. This provides absolute stellar masses and radii as outputs instead of relative radii. For cases in which RV data has high uncertainties, we used the K₁, K₂ values obtained using V2FIT, and JKTABSDIM (Southworth et al. 2005) to calculate the absolute parameters of the system. JKTABSDIM was also used to obtain distance estimates for all the targets.

JKTEBOP provides a way to obtain reliable error estimates for the fitted parameters by implementing an in-built Monte Carlo (MC) algorithm. We exploited this feature and performed 25 000 MC runs for each target to obtain the parameter uncertainties. The corner plots for the key parameters are made available in Figure A.2. The correlations observed in these corner plots are well known, and are often showed in other studies of EBs (Panchal et al. 2023). The differences between each corner plot are likely due to the different character of the LC: some show a total eclipse, while others are only grazing, and sometimes eclipses are nearly identical and at other times they are very different, as in the case of TIC 386622782.

5. Age and evolution

Stellar mass is arguably the most crucial parameter determining the evolution of a star. Therefore, stellar mass and the corresponding radii obtained through LC modelling are good parameters to compare against the stellar isochrones. However, the mass–radius plane is degenerate to the combination of age–metallicity values, requiring additional constraints. This problem can be relieved by checking if the star is located at the same evolutionary stage in the mass–T_eff (Hełminiak et al. 2019a).

In order to perform a simultaneous minimisation, we implemented a multi-parameter fitting procedure based on the Mahalanobis distance (Mahalanobis 1936). Using the EZMIST³ package, we created a database of isochrones from MESA Isochrones and Stellar Tracks (MIST) (Dotter 2016; Choi et al. 2016), from ages of 1 Myr to 10 Gyr in steps of 0.1 on a logarithmic scale. The range of metallicities from −0.25 to 0.25 was covered, and the initial abundance value was Z = 0.0142.

For the fitting procedure, we created a covariance matrix for the combination of masses, radii, and T_eff, independently for each star in the binary. We then calculated the Mahalanobis distance for this point in phase space to every row of the isochrone database. The minimum Mahalanobis distance, corresponding to the minimum χ², was then used to visualise the isochrones around the corresponding age.

Since the two stars in the binary are being fit independently within this procedure, it is possible to obtain different ages for both stars. The assumption of co-evolution is then to be applied by the user to choose from the best-fit ages around the obtained results for the two stars.

We tested this procedure on two systems, KIC 7821010 and KIC 9641031, from the Hełminiak et al. (2019a) study, to validate its effectiveness. The age estimates for these targets from the study were reproduced within 2-σ bounds around the parameters. The procedure is especially useful due to its flexibility with respect to the number of stars and the number of parameters that can be used for each of the stars. Within the scope of this paper, we did not compute any statistical errors within the fitting routine. Instead, we derived isochrones corresponding to the parameter uncertainty limits and used these to define the uncertainty range for the age estimates.

6. Results and discussions

In this section, we present the results from our analyses for the four new candidate δ Scuti pulsators in EBs, and also revised parameters for TIC 308953703 using TESS photometry.

6.1. TIC 81702112 = HD 97329

TIC 81702112 was classified as a variable star in the All Sky Automated Survey (ASAS) catalogue of variable stars (Pojmanski 2002) as ASAS J111134−4956.2. A preliminary LC+RV solution, based on ASAS photometry and RV measurements of a precision lower than ours, has been presented in Kozłowski et al. (2016). The system was revisited by the TESS, obtaining 2-min cadence photometry in sectors 10, 63, and 64. A brief look at the LC reveals that at least one component of this EB is pulsating.

An orbital solution for this target was obtained using ten CHIRON spectra. We initialised the JKTEBOP model using q, e, and ω obtained using the V2FIT, attempting a simultaneous LC-RV fit. However, after several unsuccessful attempts, we realised that the RVs of the secondary star with large rms may be responsible for the erroneous fitting outcomes. This large rms may be due to the high-amplitude pulsations, hinting that the secondary star may be the pulsator. We then used RVs only from the primary component and a good fit was obtained. We verified that the parameters for this single-component RV fit were consistent with the ones from V2FIT.

Due to the high-amplitude pulsations affecting the eclipse profiles, and hence the estimation of stellar radii, it was necessary to account for them during the modelling. We extracted the dominant pulsation frequencies using PERIOD04 (Lenz & Breger 2005), and employed a total of three sine waves using these fixed frequencies. The starting epoch and amplitude of these sine waves were kept variable along with (J), (i), (r_A+r_B), (r_A/r_B), the time of superconjunction (t_o), and the orbital period (P_orb). The best-fit model obtained after following this treatment is shown in Fig. 5.

Fig. 4. Model spectra (in red) overplotted on the disentangled spectra (in blue) of the primary and secondary components.

Fig. 5. 120-s cadence TESS LC of TIC 81702112 (filled circles), overplotted with the best-fit JKTEBOP model. The variable amplitude of pulsations was fitted using a combination of sines.

Using the best-fit parameters from the final model with pulsations, we made another model without the addition of sinusoidal functions. This “clean” eclipsing model was then removed from the LC, and the residuals were used for pulsation analysis. This model is shown in Fig. 6a.

Fig. 6. JKTEBOP binary models without using sines to account for pulsations. On the left panel is the primary eclipse and on the right is the secondary eclipse, with residuals plotted in the panel below. The error bars represent 3-σ regions around the mean parameter values. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.

Through multiple attempts made for the error estimation via MC runs, it was clear that reliable parameters were obtained when sines were simultaneously fitted to account for higher-amplitude pulsations, highlighting the utility of this feature in JKTEBOP. The errors resulting from the MC runs translate to 1.8% and 1.4% errors on the absolute radii of the primary and secondary, respectively. The orbit is eccentric, which may have a role to play in the variability of the amplitude of the pulsation frequency.

The spectroscopic analysis in ISPEC confirms that the stars have almost identical [M/H], very similar T_eff, and vsin(i). Using the parameter estimates from EB modelling and ISPEC, we determine the age of this system to be ∼1.125 Gyr, as is displayed in Fig. 7a. Both components seem to be at the MS stage; however, the more massive one is at the very end of this phase. The estimated distance lands around 303 ± 7 pc, which is in excellent agreement with the GDR3 value of 300 ± 2 pc. This is supporting evidence towards validation of the obtained radii and T_eff. An E(B − V) of 0.005 was assumed to get a consistent distance between different photometric bands within JKTABSDIM.

Fig. 7. Results from the isochrone fitting routine, indicating the positions of the stars in the mass-radius (top) and mass–(Teff) (bottom) diagrams. The black line indicates the best age estimate obtained for the set of input stellar parameters. The filled red and blue circles indicate the more and less massive components, respectively. In the case of TIC 165459595, the error bars represent the 2-σ region instead of 3-σ. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.

The synchronisation velocities given by JKTABSDIM, for the less massive and more massive components, are ∼17 and 21 km/s, respectively. These values are valid for circular orbits. In the case of eccentric orbits, a pseudo-synchronisation is achieved when the ratio of orbital period and rotation period equals

$$\begin{array}{c}\hfill \frac{P_{\mathrm{orb}}}{P_{\mathrm{rot},\mathrm{ps}}}=\frac{1+7.5e^2+5.625e^4+0.3125e^6}{{(1-e^2)}^{3/2}(1+3e^2+0.375e^4)}.\end{array}$$(1)

Using Eq. (1) (Hut 1981; Mazeh 2008), we derived the pseudo-synchronisation velocities (v_rot, ps) for the less massive and more massive star to be ∼21 and 26 km/s, respectively. The values are in 2-σ agreement with the vsin(i) values for these components obtained from the spectral analysis.

6.2. TIC 189784898 = V Cir

The classification of TIC 189784898 as a variable star dates back to 1932 (Gaposchkin 1932). This variability was later accredited to eclipses (Strohmeier 1967; Alfonso-Garzón et al. 2012). More recently, Zasche (2010) presented LC analysis of data obtained with the INTEGRAL-OMC instrument. The target was revisited by the TESS and 2-min cadence photometry was obtained in sector 38.

We initiated the simultaneous LC-RV modelling using the orbital parameters from V2FIT. The LC demanded a single polynomial together with a total of four sine functions alongside the binary model to account for trends likely arising due to stellar spot(s). These were fixed to P_orb, P_orb/2, P_orb/3, and P_orb/4, while keeping the starting epoch and amplitude free. Initially, we attempted to fit two sines per component to check their contribution towards the final fit; however, we found that spots on the primary component were preferred in order to obtain the best fit. Due to the lower number of RVs and variability near eclipses due to spots, we chose to fix e and ω to zero. The best-fit obtained following this process is shown in Fig. 6b.

The ISPEC analysis reveals a hotter primary that is rotating more slowly than the larger secondary component. The metallicity estimates of the two stars match. It is slightly lower than the solar value; however, the errors are large enough to prevent us from saying this with absolute certainty. The high S/N disentangled spectra also allowed us to check the spectroscopic log(g), which agreed with the photometric log(g) within errors for both the components.

The distance calculated using JKTABSDIM, 575 ± 23 pc, agrees with the GDR3 distance of 549 ± 5 pc. To achieve a consistent distance between different bands within JKTABSDIM, the E(B − V) of 0.2 needed to be used. The age for this system is estimated to be ∼1.75 Gyr. As is shown in Fig. 7b, the more massive component has already evolved into a sub-giant, and the less massive star is at the end of its MS.

6.3. TIC 386622782 = HR 2214

The discovery and classification of TIC 386622782 as a “double star” dates back to 1950 (Muller 1950; Wilson 1950). The system is a visual binary, with one of the components being an SB2. This object has been observed over multiple campaigns (Cannon & Pickering 1993; Royer et al. 2002; Pourbaix et al. 2004; Adamczak & Lambert 2014) in attempt to constrain the spectral type and the rotational and radial velocities. It is described across the literature as a chemically peculiar star (Slettebak & Nassau 1959; Abt & Levy 1985; Renson & Manfroid 2009; Ghazaryan et al. 2018). The eclipsing (and pulsating) nature of this system was first noted during an inspection of the TESS data by M. Pyatnytskyy in January 2022, and reported to the International Variable Star Index (PMAK V115⁴).

A first look at its TESS LC shows primary eclipses at a temporal separation of 23.80841 days, and a complex pulsation signal containing multiple frequencies. A careful inspection of the phase-folded LC reveals a tiny secondary eclipse at a phase of ∼0.49. Modelling the tiny secondary eclipse in JKTEBOP was challenging because the amplitude of pulsations is of the same order as that of the eclipse depth. Fixing the value of e and ω obtained from V2FIT, we performed a simultaneous LC-RV fit. This was necessary to make sure that the model does not “miss” the tiny secondary eclipse. We primarily aimed to constrain the J, i, r_A, and r_B. A total of six polynomials were used to offset the trends affecting different sections of the LC. The best-fit model, as is seen in Fig. 6c, is able to reproduce the tiny secondary eclipse. The geometry of this eccentric system gives rise to such a small secondary eclipse in spite of this star contributing more than 80% of the total flux in the observed band. In this case, the availability of a large number of RV measurements played a crucial role in obtaining a reliable binary model. The distribution of the parameters during MC runs is displayed in a corner plot in Fig. A.1.

The metallicity values for the two components are consistent with each other. Using the obtained T_eff estimates, the distance to this system was found to be 103 ± 3 pc. This value lies just within the 3-σ range of the last updated value from the Hipparcos catalogue (van Leeuwen 2007). We found that an isochrone of age ∼0.8 Gyr reproduced the two stars in fair agreement with the parameters obtained from our analysis. However, the temperature of the less massive component is nearly 4-σ away from the corresponding value on the isochrone. This isochrone fit is shown in Fig. 7c.

6.4. TIC 165459595 = V1109 Cen

During our sample selection from the CRÉME database, TIC 165459595 was included as one of the components was near the 2.5 M_⊙ upper limit. While collecting further data on the sample, it was noted that TIC 165459595 has recently been flagged as one of the many EBs observed by TESS that exhibit pulsations (Shi et al. 2022; Chen et al. 2022). A detailed exploration of its physical as well as oscillatory nature is therefore timely.

An orbital solution for this target was obtained using six CHIRON spectra. The unsatisfactory quantity of spectra and high rotational velocities of both components hampered the calculation of precise K₁, K₂, resulting in much larger errors on masses compared to the other targets in the sample. We initialised the JKTEBOP model using q, e, and ω obtained from V2FIT, attempting a simultaneous LC-RV fit using the 2-min cadence TESS photometry from sector 37. The MC errors of ∼0.3% obtained on the relative radii translated to ∼1.4% errors on the absolute radii of the stars. The best-fit model is shown in Fig. 6d

In spite of the noisy nature of composite spectra and the low number of observations, we tried disentangling the components and modelling them within ISPEC. In this case we fixed all of the parameters, except the T_eff and v_mic. After a few unsuccessful attempts, we chose not to fit for the [M/H]. The two-fold reasons for erroneous results could be high rotation velocity and relatively low S/N of the decomposed individual spectra. To determine the T_eff, the [M/H] was set to the GDR3 estimate of −0.44. The bottom panel of Figure 4 shows that the model obtained using ispec (red) does not fit very well, particularly in the 535–538 nm section. This could be due to poor S/N or incorrect continuum normalisation or both.

The process of determining the age for this system was challenging due to the unavailability of reliable spectroscopic temperature constraints. We attempted to constrain it using the flux ratio obtained from the EB modelling. We did this using ISOFITTER, explained in further detail in Moharana et al. (2023), which searches through a grid of MIST isochrones and finds the best fit for the given set of parameters. Keeping the l₃ parameter free results in a L_B/L_A value of ∼0.9, meaning the less massive primary star is brighter. This corresponds to an age of ∼2.37 Myr. However, when l₃ was kept fixed to zero, we obtained an L_B/L_A that is slightly greater than 1. The difference between the residuals obtained after subtracting these two models from the LC is negligible. This prefers the slightly older age of ∼100 Myr. The errors in the mass estimates for the stars are quite large, making it very difficult to better constrain the age. As is shown in Fig. 7d, the best-fit isochrone represents stars with temperatures around 13 000 K.

The calculation of distance for the other targets, within JKTABSDIM, used the Kervella et al. (2004) relations. However, these are valid for temperatures up to 10 000 K. The distance calculated using the relations of Girardi et al. (2002) is ∼400 pc and does not agree with the GDR3 distance of 310 ± 2. Also, consistency between the different photometric bands is obtained only if the E(B − V) is set to ∼0.45.

6.5. TIC 308953703 = HIP 7666

TIC 308953703 has been revisited multiple times over the past 40 years (Oja 1985; Cannon & Pickering 1993; Liakos & Niarchos 2017) and is classified as an EB with one of the components exhibiting pulsations (Escolà-Sirisi et al. 2005). Recently, Feng et al. (2021) performed an extensive multi-site photometry and spectroscopic study to derive the absolute parameters of the binary, and study the pulsations in further detail. Three significant pulsation frequencies were identified and mode identification was attempted. However, the authors suggested that photometry from space missions is needed to improve the reliability of these mode identification results.

This target now has a 2-min cadence TESS (sector 58) photometry. We truncated the LC to four orbits to optimise computational time, while obtaining a good resolution between pulsation frequencies. We initialised our LC model together with the RV data obtained by Feng et al. (2021). During the initialisation of our model, we identified a typographical error in one of the RV timestamps provided. Specifically, the value 329.75669 was incorrectly recorded as 331.75669. The RV data has insufficient phase coverage, which adds the possibility of a small but measurable value of e. To check this, we set e as a free parameter, at a fixed omega of 90° (expected from the apparent phase of the secondary minimum). Due to the lower amplitude of pulsation frequencies, it was possible to obtain a good fit without implementing additional sines during the JKTEBOP runs, except one for the out-of-eclipse ellipsoidal variability.

A good fit is obtained at e = 0.02697, and it is shown in Fig. 6e. This resulted in slightly different values of masses and radii compared to the Feng et al. (2021) results. These updated parameters were used to determine the age of the system. Feng et al. (2021) discussed the inconsistency between the parameter estimates for the cooler secondary star and the one from the MESA evolutionary tracks. As is shown in Fig. 7c, the updated masses and radii for the primary and secondary can now be represented by a single isochrone of age ∼1.75 Gyr, where both the stars are on the MS.

We used the temperature estimates together with the updated parameters to estimate the distance to this system. The value provided by JKTABSDIM for the J-band – 292 ± 7 pc – is consistent within 3-σ of the GDR3 value of 313 ± 1.3.

7. Pulsation analysis

To determine the pulsation frequencies for each system, we used the residuals left after subtracting the binary model from the LC. This ensures that effects arising due to ellipsoidal deformation and stellar spots are also removed to a good extent. We used v.1.2.9 of the PERIOD04, based on Lenz & Breger (2005), to perform the time series analysis. A typical pre-whitening procedure was followed to extract all the frequencies. The S/N ratio for detection, after pre-whitening for the previous frequencies, was kept at 5. All the significant frequencies were then simultaneously fitted using sinusoids in order to obtain their amplitudes and phases. A list of all the extracted frequencies, with their analytical errors (Montgomery & O’Donoghue 1999), for each of the systems is provided in Appendix C.1. Uncertainties on amplitude and phase were calculated using an MC simulation within PERIOD04, keeping the frequency fixed to uncouple frequency and phase uncertainties.

7.1. TIC 81702112

The strongest pulsation frequency has an amplitude of ∼10 mmag (see Fig. 11a), which is variable over the orbital phase. Such variability in the amplitude of pulsation frequency has been discussed in a few recently discovered systems (Handler et al. 2020; Fuller et al. 2020). The extent of this amplitude modulation over the orbital phase can be seen in Fig. 8.

Fig. 8. Amplitude of frequency with the highest power plotted over the orbital phase.

Suffering from the lack of total eclipses, it is difficult to state with confidence which of the two stars is pulsating. However, to get an idea, we followed the technique of analyzing the Fourier spectra in the primary and secondary eclipses. We chopped the residuals from the binary modelling to obtain these eclipse snippets, and overplotted their periodograms on top of each other as is shown in Figure 9. The amplitude of pulsations in the primary eclipse is slightly larger compared to the ones in the secondary eclipse. This is not so evident when visually inspecting the eclipses because of the difference in depths. The difference in pulsation amplitudes during the eclipses, though not large, hints towards the secondary being the pulsator.

Fig. 9. Frequency power spectra for primary (blue) and secondary eclipses (red) of TIC 81702112. Overplotted in black is the frequency spectra of the out-of-eclipse portions of the LC.

7.2. TIC 189784898

The target was selected based on the stellar masses of the two stars. Although the LC did not show clear evidence of pulsations, initial checks did hint towards a couple of very low-amplitude pulsation signals. The binary subtracted residuals, however, show two very faint peaks (S/N ∼ 3; see Fig. 11b) that are extremely close to the noise floor. These peaks fall below our set S/N limit. We therefore mark this as a false positive.

This target has 2-min cadence data only from a single TESS sector so far. More photometry will help in boosting the S/N, which is necessary to test the nature of these low-amplitude frequency peaks.

7.3. TIC 386622782

The large number of independent and combination frequencies, mentioned in Table C.2, result in the complex pulsations visible in the LC. Considering the large flux contribution from the secondary, it is most probable that it is the pulsating component. However, it cannot be ruled out that the other star may also be pulsating, although the relative contribution of those pulsations to the visible frequency spectra would be minimal.

7.4. TIC 165459595

The frequencies for this target fall in the range of δ Scuti-type oscillations. Stars with an initial mass greater than 1.5 times M_⊙ can cross the instability strip before entering the MS (Handler 2009). Such oscillations can therefore be explained despite of the young age of the system.

The T_eff and log(g) for δ Scuti-type stars typically have values between 6300 and 8500 K and 3.2–4.3, respectively, corresponding to the A-F type variables. However, for this target the T_eff estimates from the isochrone fits and the log(g) from the LC modelling would suggest a B-type classification. One such candidate with higher T_eff, AL Scl, was recently found by Kahraman Aliçavuş et al. (2023a). Better constraints on the atmospheric parameters are necessary to explain the pulsations observed in this case.

7.5. TIC 308953703

The extracted independent frequencies, f₁ and f₃, and their amplitudes are close to the values obtained by Feng et al. (2021). Since they have already performed a detailed asteroseismic study for this target, we have not attempted any further analysis.

Using the P_orb values from LC modelling and the period corresponding to the highest-amplitude frequency, we located our targets on the log P_orb–log P_pul diagram (see Fig. 10). We overplotted the P_orb − P_pul empirical relation obtained by Liakos (2020), and used it check the expected P_pul for TIC 189784898.

Fig. 10. Location of the dominant pulsating frequency for the studied systems relative to the Porb − Ppul empirical relation (dotted line) derived by Liakos (2020), for such systems with Porb < 13 d.

Fig. 11. Periodograms for residuals obtained after subtracting JKTEBOP models from the corresponding LCs. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.

8. Conclusions

We have presented a comprehensive analysis for four candidate DEBs with at least one component exhibiting δ Scuti-type pulsations. We have confirmed the presence of pulsations in three systems, and obtained a set of precise parameters for both stars. We have also presented updated parameters for TIC 308953703 based on the LC modelling of the short-cadence TESS photometry.

The variability of the pulsation frequency in the case of TIC 81702112 makes this target an extremely interesting case study. However, detailed modelling to explain the mechanism responsible for this is beyond the scope of this paper. It is necessary to use a combination of techniques, such as multi-band photometry for mode identification, detailed binary evolution models, and asteroseismic modelling, to obtain a better understanding of this unique target.

We are able to put good constraints on the physical properties and evolutionary stage of all the targets in the sample except for TIC 165459595. The main limitation in this case is the unavailability of high-quality spectra. However, the obtained masses, radii, and flux ratio were enough to estimate the stellar ages in the range of a few million years. Our understanding of δ Scuti-type pulsations in pre-MS stars is still in the nascent stages. Observational studies of such targets are crucial to aid the ongoing theoretical studies (Steindl et al. 2021; Murphy et al. 2023) of such stars.

Building upon the detailed study by Feng et al. (2021), we have been able to update the physical parameters for TIC 308953703 and prove that the evolutionary stages of both stars can in fact be explained within a co-evolution scenario. This study attempts to further enrich the sample of well-studied δ Scuti-type pulsators in EBs. Such a sample is crucial to understand the extent of the effect of binarity on stellar pulsations. With high-quality space photometry, the detection of such systems is rising rapidly. This presents a great opportunity to use these unique and powerful tools to study the structure and evolution of such stars in great detail.

Acknowledgments

T.P. would like to thank A. Miszuda for the fruitful discussions and suggestions. T.P. acknowledges the support provided by the Polish National Science Centre (NCN) through grant no. 2021/43/B/ST9/02972 and 2017/27/B/ST9/02727. K.G.H. acknowledges the support provided by the NCN with the grant 2023/49/B/ST9/01671. A.M. acknowledges the support provided by the NCN with the grant 2021/41/N/ST9/02746. This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018). This research uses data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by NASA’s Science Mission directorate.

References

Appendix A: JKTEBOP MC runs

Fig. A.1. Corner plot for the distribution of parameter values from the MC runs within JKTEBOP. For i we use 95% of the sampling to avoid values where the secondary eclipse is left undetected.

Fig. A.2. Corner plots for the other targets corresponding to their best-fit JKTEBOP models.

Appendix B: Additional Table

Appendix C: Pulsation Frequencies

All Tables

All Figures

	Fig. 1. Short (2-min) cadence LCs of the targets, each for a single TESS sector.
In the text

	Fig. 2. Orbital fits obtained for the RV curves for the targets. Here, primary refers to the more massive star and secondary is the less massive companion. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595.
In the text

	Fig. 3. Input and output spectra utilised and obtained during the spectral disentangling for the target TIC 81702112. Black: one of the observed spectra that was used as an input for the disentangling algorithm. This is a composite of contributions from both the stars. Blue: disentangled spectrum of the primary component. Red: disentangled spectrum of the secondary component.
In the text

	Fig. 4. Model spectra (in red) overplotted on the disentangled spectra (in blue) of the primary and secondary components.
In the text

	Fig. 5. 120-s cadence TESS LC of TIC 81702112 (filled circles), overplotted with the best-fit JKTEBOP model. The variable amplitude of pulsations was fitted using a combination of sines.
In the text

	Fig. 6. JKTEBOP binary models without using sines to account for pulsations. On the left panel is the primary eclipse and on the right is the secondary eclipse, with residuals plotted in the panel below. The error bars represent 3-σ regions around the mean parameter values. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.
In the text

Fig. 7. Results from the isochrone fitting routine, indicating the positions of the stars in the mass-radius (top) and mass–(Teff) (bottom) diagrams. The black line indicates the best age estimate obtained for the set of input stellar parameters. The filled red and blue circles indicate the more and less massive components, respectively. In the case of TIC 165459595, the error bars represent the 2-σ region instead of 3-σ. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.

In the text

	Fig. 8. Amplitude of frequency with the highest power plotted over the orbital phase.
In the text

	Fig. 9. Frequency power spectra for primary (blue) and secondary eclipses (red) of TIC 81702112. Overplotted in black is the frequency spectra of the out-of-eclipse portions of the LC.
In the text

	Fig. 10. Location of the dominant pulsating frequency for the studied systems relative to the Porb − Ppul empirical relation (dotted line) derived by Liakos (2020), for such systems with Porb < 13 d.
In the text

	Fig. 11. Periodograms for residuals obtained after subtracting JKTEBOP models from the corresponding LCs. (a) TIC 81702112. (b) TIC 189784898. (c) TIC 386622782. (d) TIC 165459595. (e) TIC 308953703.
In the text

	Fig. A.1. Corner plot for the distribution of parameter values from the MC runs within JKTEBOP. For i we use 95% of the sampling to avoid values where the secondary eclipse is left undetected.
In the text

	Fig. A.2. Corner plots for the other targets corresponding to their best-fit JKTEBOP models.
In the text