© V. Vaulato et al. 2022

1 Introduction

Since the pioneering detections of 51 Peg b, the first exoplanet orbiting a solar-type star (Mayor & Queloz 1995), and HD209458b, the first transiting hot Jupiter (Henry et al. 2000; Charbonneau et al. 2000), contemporary exoplanetary research has mostly focused on late-type (FGKM) main-sequence hosts. This is not just because of the interesting analogies with our own Solar System, but also because these stars are particularly suited for the most fruitful detection techniques in terms of discoveries: radial velocities (RV) and transits. These two techniques are hampered by the presence of stellar activity and pulsations, making the detection and characterization of planets around very young, evolved, or early-type stars extremely challenging. Along these lines, the somewhat standard road to discovery has become the photometric detection of transit-like signals, followed by an RV monitoring to confirm the planetary nature of the candidate and to measure its mass.

The vast majority of transiting planets known today have been discovered by dedicated space-based telescopes: starting with the pioneering work of the CoRoT satellite (2006–2013; built by an international consortium led by CNES; Auvergne et al. 2009), which was followed by the Kepler mission by NASA (2009–2014; Borucki et al. 2010). Kepler was stopped because of a technical failure and was then restored for a second-phase observing program called K2 (2014–2018; Howell et al. 2014). Finally, the advent of the NASA Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), launched in 2018 and currently operating, doubled the number of known candidate exoplanets by starting the first all-sky space-based transit survey. TESS is optimized to search for Earth- to Neptune-sized planets transiting bright and nearby main-sequence stars with a particular focus on cool, red dwarfs (K and M) because their smaller radii and lower masses make the transit and RV signals much stronger. TESS objects of interest are three magnitudes brighter than those of Kepler on average (Barclay et al. 2018), enabling an easier and more effective spectroscopic follow-up analysis. Most importantly in this context, and unlike its predecessors, TESS has the capability to download photometric data not only for a limited sample of preselected targets, but also for virtually any other bright star in the sky, although at a longer cadence with respect to the core sample. This is made possible by the availability of full-frame images (FFIs), originally downloaded every 30 min, and since sector 27, at a 10-min cadence. FFI photometry opened a vast array of opportunities in many astrophysical fields, including extending the search for planets to a wider range in the stellar parameter space (Montalto et al. 2020; Nardiello et al. 2019).

Several alternative detection techniques have been proposed and implemented to search for exoplanets around stars other than solar-type or M dwarfs, or in different Galactic environments. A-type stars, in particular, are one of the least explored class of planet-hosting stars. Of more than 5000 confirmed planets with known mass and/or radius listed in the NASA Exoplanet Archive, only 22 are hosted by 21 main-sequence stars earlier than the F0V spectral type (T_eff ≥ 7220 K; Pecaut & Mamajek 2013). With only four exceptions coming from direct imaging surveys, they were all detected through the transit technique (3 by HAT-Net, 3 by WASP, 5 by KELT, 2 by MASCARA, 3 by Kepler, 2 by TESS) and were confirmed by means of spectral tomography (Collier Cameron et al. 2010) and/or statistical validation techniques because classical RV measurements are usually not effective on such stars, whose spectral lines are few and Doppler broadened by their very high rotational velocities (v sin i ~ 100 km s⁻¹; Rainer et al. 2021).

Only three of the hosts described above are confirmed or suspected δ Scuti variables, that is, population I pulsators located in the classical Cepheid instability strip (Bowman & Kurtz 2018; Guzik 2021 and references therein): WASP-33b (Herrero et al. 2011), HAT-P-57b (Hartman et al. 2015), and WASP-167b = KELT-13b (Temple et al. 2017). Because ~50% of the stars in this region of the parameter space are expected to show δ Scuti pulsations to some extent (Bowman 2017), some astrophysical reason must explain the lack of planets around δ Scuti stars, and/or a strong observational selection effect is at play. The latter explanation is likely true, as the stellar pulsation pattern is known to limit the efficiency of transit-search algorithms even after they are filtered in the frequency domain by sophisticated techniques (Ahlers et al. 2019; Hey et al. 2021).

In this context, the pulsation timing technique (PT; Hermes 2018) is a completely independent method for detecting massive companions around stars showing coherent and well-detectable pulsations, such as δ Scuti stars, pulsating white dwarfs, or hot subdwarfs. The underlying physical principle is the so-called light-travel-time effect (LTTE) that was first discussed in detail by Irwin (1952): Because the speed of light c is finite, the motion of the target along the line of sight will translate into a phase shift of any periodic astrophysical signal originating from the star, including pulsations. The arrival times of the signals can then be compared with those predicted by a linear ephemeris (implying a constant periodicity dP/dt = 0) by computing the so-called O–C diagram (observed-calculated; Sterken 2005). Any deviation from a strict periodicity dP/dt ≠ 0 (assuming that the pulsation themselves are intrinsically stable) would then reveal itself as a nonzero second derivative in the O–C diagram (see Sect. 3 for details). The LTTE signal from a single, massive perturber on a Keplerian orbit is formally analogous to an RV curve having half of the eccentricity e and the argument of the periapsis ω decreased by π/2 (Irwin 1952); if the full shape of the LTTE signal is retrieved, then the orbital period P, the minimum mass M sin i, and other orbital elements of the perturber can be reliably measured. The PT approach shares some similarities with the astrometric method (Black & Scargle 1982) in that it measures the displacement of the target star with respect to an inertial rest frame, but projected along the line of sight rather than on the sky plane. Notably, also the underlying selection effects are similar because both methods are more sensitive to massive bodies at large star–planet separations. PT therefore probes a region of the parameter space that is complementary to the space that is investigated by transits and RVs, which are more sensitive to close-in perturbers (Sozzetti 2005).

The PT technique requires photometric data spanning a temporal baseline that is long enough to sample the orbital period of the perturber and a signal-to-noise ratio (S/N) high enough to constrain the phase shift of the oscillations to a level comparable to the expected LTTE amplitude¹. The main pulsation modes of the star have to be individually identified and measured on the Fourier spectrum of the light curve, implying that continuous time series spanning at least ~10 days are mandatory to achieve the needed resolution in the frequency space (~0.1 cycles per day) when dealing with the harmonic content of typical δ Scuti stars (Murphy et al. 2014). This is the main reason why sparse sampling imposed by the day-night cycle is the main limiting factor of ground-based observations. Targets with an extremely simple harmonic content, such as pseudo-sinusoidal pulsators, represent an exception, and a few pioneering results were published by analyzing ground-based data (Paparo et al. 1988; Barlow et al. 2011b). Interestingly, the latter claim was subsequently confirmed by an independent RV follow-up, demonstrating the reliability of this technique (Barlow et al. 2011a).

The full power of the PT technique, however, revealed itself when the nearly uninterrupted 4 yr baseline of Kepler photometry was exploited. Some works (Shibahashi & Kurtz 2012; Murphy et al. 2014, 2018, 2020; Shibahashi et al. 2015; Murphy & Shibahashi 2015), although they focused on the detection and characterization of stellar-mass companions, demonstrated that under favorable assumptions, even LTTE signals from planetary-mass companions (<13 M_jup) are in principle detectable through space-based photometry (Murphy et al. 2014). Murphy et al. (2016) presented the PT detection of a massive planetary candidate for the first time (11.8 ± 0.7 M_jup, P = 840 ± 20 days). This candidate orbits a metal-poor δ Scuti star (KIC 7917485). Together with the discovery by Silvotti et al. (2007) of a giant planet around the hot subdwarf B star V391 Peg, these are the only PT discoveries of planetary-mass objects published so far. It should be noted, however, that a follow-up paper on V391 Peg with new data (Silvotti et al. 2018) failed to fully reproduce the previous claim, probably due to nonlinear interactions between pulsation modes (see also Bowman et al. 2021). This acts as a caution.

While several studies of δ Scuti stars have been published that exploited TESS short-cadence data (Antoci et al. 2019; Hasanzadeh et al. 2021; Southworth 2021 for a review), none of them applied the PT technique. The goal of this paper is to focus on the PT analysis of a small sample of particularly suited targets observed by TESS at short cadence as a pilot study for the exploitation of TESS (and, later on, PLATO) data on a larger scale. In Sect. 2, we describe how our targets were identified and how their relevant stellar parameters were collected. In Sect. 3, we summarize all the steps of our data analysis, starting from the filtering of raw light curves through the harmonic analysis to measure the phase shifts as a function of time. We determine the best-fit orbital solutions through an LTTE model in Sect. 4 and discuss the results in Sect. 5.

2 Target selection and characterization

The vast majority of targets observed by TESS at short and ultra-short photometric cadence (120 and 20 s, respectively) are FGKM dwarfs selected for the core transit-search survey, drawn from the so-called candidate target list (CTL; Fausnaugh et al. 2021). However, stars of different spectral types, including pulsating variables, are ingested mostly from GI/DDT proposals or from the TESS Asteroseismic Science Operation Center (TASOC²), which provides the TESS Asteroseismic Science Consortium (TASC) with an asteroseismological database of the mission (Schofield et al. 2019).

In order to select our targets, we chose the catalog compiled by Chang et al. (2013) as a starting point. This catalog lists a sample of 1578 δ Scuti stars that were identified at high confidence level by several previous surveys. We cross-matched this catalog with the TESS CTL tables for sectors from 1 to 42 included to obtain only δ Scuti stars that were observed by TESS at 2-min cadence. As an additional constraint, we further restricted our sample by (1) requiring a minimum data coverage of at least seven TESS sectors, not necessarily contiguous, to build O–C diagrams with a number of points and a temporal baseline long enough to detect LTTE signals at an approximately 1-yr timescale (see Sect. 4); (2) selecting targets brighter than T = 12 in the TESS photometric system, to avoid being limited by photon and background noise; and (3) excluding all the stars that have been identified as binaries in the Chang et al. (2013) catalog. These additional constraints are justified by the nature of this study, which is not focused on a complete sample, but is rather intended as a pilot study to investigate the limiting factors of TESS (especially due the systematic errors) on a very small sample of the most favorable targets.

The final output of this selection process is a short list of 12 targets, which were then individually examined both through a literature search and by inspecting their TESS light curves. In particular, we confirmed that all of them are actually δ Scuti pulsators and carried out a preliminary harmonic analysis on them by computing the generalized Lomb–Scargle periodogram (GLS; Zechmeister & Kürster 2018) of a single sector with the same algorithms as described in Sect. 3. While most targets are classical δ Scuti stars showing small-amplitude (~0.01–0.03 mag) pulsations and a complex pattern of different modes in the frequency domain, two of them stand out as high-amplitude pulsators (≥0.1 mag) with a particularly clean periodogram, in which most of the signal is due to a single, well-defined pulsation mode and its harmonics (Figs. 1 and 2): hereafter we refer to them as Chang 134 (from its ID number in the Chang et al. 2013 catalog) and V393 Carinae. This combination of a large pulsation amplitude and a very simple and coherent frequency pattern makes these targets well suited for a precise timing analysis, especially for our pilot study. For this reason, we focused this analysis on them and set aside the other targets (five of which are multi-mode pulsators, but with modes that are easy to identify) for the next paper of the series. A short review of the available scientific literature is given in the following subsections, and the basic astrophysical parameters of the two targets are reported in Table 1. It is worth mentioning that both our targets could fit the high-amplitude δ Scuti subclass (HADS; Breger 2000; Antoci et al. 2019), which is characterized by large-amplitude pulsations, low rotational velocities, and time-domain spectra that are dominated by radial modes. We did not attempt to perform a detailed identification of their pulsation modes, however, because our work is focused on a dynamical technique (where the pulsation pattern is exploited just as a coherent astrophysical clock) rather than on stellar physics.

Fig. 1 TESS photometry of Chang 134 = TYC 9158-919-1. Upper panel: one-day section of the light curve from TESS sector 1. Lower panel: GLS periodogram of the whole light curve including nine TESS sectors, stitched together and filtered as explained in Sect. 3.

2.1 Chang 134 = TYC 9158-919-1

Chang 134, also known as TYC 9158-919-1 or TIC 431589510, was first identified as a variable star by ASAS (Pojmanski 2002) and classified as a generic δ Scuti star (DSCT variability class). No targeted follow-up study on it has been published since 2002. An asteroseismological measurement of log(g) and T_eff was included in the catalog by Barceló Forteza et al. (2020) through the v_max seismic index extracted from the TESS light curves in an automated fashion. All the astrophysical parameters we collected are reported in the second and third columns of Table 1.

The PT technique, just like astrometry or RVs, cannot directly infer the mass of perturbing body, but rather its ratio m_p/M_* with the stellar mass. The latter is needed at the analysis stage (Sect. 4) to properly interpret our results for the two targets. For Chang 134, an estimate of stellar mass (M_* = 1.40 ± 0.16 M_⊙) and age was first published in the catalog by Mints & Hekker (2017) using UniDAM models, but without taking advantage of spectroscopic data or Gaia parallaxes, which were not available at that time. Later, Chang 134 was included in the large-scale analysis by Queiroz et al. (2020), combining high-resolution spectra from the APOGEE-2 survey DR16 with broadband photometric data taken from different sources and Gaia DR2 distances. Stellar parameters are derived as the posterior distribution returned by the Bayesian isochrone-fitting code StarHorse (Queiroz et al. 2018). The stellar mass reported for Chang 134 is M_* = 1.38 ± 0.03 M_⊙, where the central value is the 50th percentile of the StarHorse posterior distribution, while the symmetrized uncertainty is calculated as the half-difference between the 84th and the 16th percentiles of the same distribution.

Fig. 2 TESS photometry of V393 Carinae = HD 66260. Upper panel: one-day section of the light curve from TESS sector 1. Lower panel: GLS periodogram of the whole light curve including 14 TESS sectors, stitched together and filtered as explained in Sect. 3.

2.2 V393 Carinae = HD 66260

Unlike Chang 134, the much brighter V393 Car (also known as HD 66260 or TIC 364399376) has been the subject of several targeted studies. It was first discovered as a variable star and classified as a δ Scuti radial pulsator by Helt (1984). Balona & Evers (1999) found in a more detailed analysis that the dominant mode is rather nonradial; the debate of this point is still ongoing. A first claim of change in the period of the main pulsation mode appeared in Garcia et al. (2001), who measured a 7.18-min discrepancy with respect to the ephemeris by Helt (1984). Interestingly, further ground-based photometric follow-up by Axelsen (2014) retrieved the very same period of Helt (1984) within 0.33 s, and found no conclusive evidence of additional overtones. The present work therefore also represents an opportunity to confirm or disprove the claimed period change, and to determine whether it is caused by a secular or an oscillating term.

A measurement of v_max, log(g), T_eff was listed by for Chang 134 by Barceló Forteza et al. (2020). Unfortunately, V393 Car is missing from the Queiroz et al. (2020) catalog, but is listed in the Anders et al. (2022) catalog, which applied StarHorse on the Gaia EDR3 photometric and astrometric measurements, yielding M_* = 2.02 ± 0.10 M_⊙ for this target. As accurate stellar parameters of V393 Car derived from spectroscopy are lacking in the literature, we attempted to obtain an independent estimate of the stellar mass as a crosscheck by applying the empirical relations derived by Moya et al. (2018) as a function of other stellar parameters. In particular, we exploited the relation $${\log }_{10}\left(M_{*}\right)=-a+b{T}_{\text{eff}}+c{\log }_{10}\left(L_{*}\right),$$(1)

where a = −0.119 ± 0.003, b = 2.14 × 10⁻⁵ ± 5 × 10⁻⁷ and c = 0.1837 ± 0.0011 are the coefficients calculated by Moya et al. (2018). This relation is valid within the range of temperatures 4780 ≤ T_eff ≤ 10990 K. Adopting T_eff from the seismic value given by Barceló Forteza et al. (2020) and L_* from the Gaia DR2 catalog³, and propagating the errors, we computed a stellar mass M_* = 1.98 ± 0.17 M_⊙ for V393 Car, which is perfectly consistent with the previous value. We therefore adopt the M_* = 2.02 ± 0.10 M_⊙ estimate from Anders et al. (2022) in the subsequent analysis.

We are aware that a detailed asteroseismic analysis performed on all the available TESS light curves would likely yield much more accurate fundamental parameters for both our targets (Hasanzadeh et al. 2021). However, since an error of about 5% on stellar mass will not be the dominant source of uncertainty in our final parameters (as we discuss in Sect. 4), an analysis like this is beyond the scope of our work.

2.3 TESS light curves

Chang 134 and V393 Car were both observed by TESS in 2-min cadence mode. At the time when our analysis started, Chang 134 was observed in 9 nonconsecutive sectors (1-2-3-6-13-27-28-29-36) between July 2018 and March 2021. V393 Car was observed in 14 nonconsecutive sectors (1-4-7-8-10-11-27-28-31-34-35-36-37-38) between July 2018 and April 2021. All the photometric data analyzed in this work were processed and extracted from the raw data by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) and are publicly available on the Mikulski Archive for Space Telescopes (MAST⁴). Specifically, we built the light curves for our analysis by extracting the stellar flux from the pre-search data conditioning simple aperture photometry column (PDCSAP; Smith et al. 2012; Stumpe et al. 2012) because it results in cleaner time series since systematic long-term trends are removed.

3 Harmonic analysis

Our approach to data analysis is based on the fit of sums of harmonic functions to preselected segments of our light curves in the form $${\displaystyle \sum _{i}f\left(A_i,P_i,{\varphi }_i\right)}={\displaystyle \sum _i{A}_i\times \cos \left(\frac{2\pi t}{P_i}+2\pi {\varphi }_i\right),}$$(2)

where A_i is the amplitude, P_i the period, and ϕ_i the phase of each component normalized between 0 and 1. The phase of the signal at the dominant pulsation mode ϕ₀ is then retrieved for each segment, giving the so-called phase shift ϕ₀(t) as a function of time. From this quantity, the usual O–C diagram in time units can be computed in a straightforward way by just multiplying the phase shift by the pulsation period P₀.

We developed an independent pipeline to perform all the needed tasks from the raw light curves down to the final O–C diagram. Its flow chart (described in the following sections) is conceptually similar to the phase modulation (PM) method described by Murphy et al. (2014), with some substantial differences, in particular, in how the actual fit is performed, that is, through a Markov chain Monte Carlo (MCMC) analysis in our case rather than with a simple least-squares algorithm. Most of our processing steps were implemented by scripting and modifying routines from the VARTOOLS code version 1.39⁵ released by Hartman & Bakos (2016), to which we refer for a more in-depth description of the individual algorithms.

3.1 Light curve preconditioning

We started our process by downloading TESS SPOC light curves of our targets from the MAST archive for each observing sector (23 combined sectors). Each data point was sampled every 2 min as a default for short-cadence observations.

As a first step, we discarded all data points flagged as bad, that is, with a quality factor q ≠ 0 (QUALITY column). We then extracted the time column TIME from each FITS file and converted it into the BJD_TDB standard (Eastman et al. 2010). Finally, we extracted the PDCSAP_FLUX column and its associated error PDCSAP_FLUX_ERR and converted them into the magnitude system. The choice of PDCSAP over SAP is justified by our need that systematic errors, especially those manifesting themselves as long-term trends, are corrected for as much as possible to avoid unnecessary noise on our periodograms. As our analysis is focused on retrieving the phase of the pulsation signal, we are not concerned by any small perturbation in the amplitudes that might be introduced by the PDCSAP processing (Cui et al. 2019). In order to filter out the most obvious outliers that survived the q = 0 selection, we carried out an iterative clipping at 15 σ with respect to the mean. This specific threshold was empirically set by confirming that the shape of the main pulsation mode was left unchanged.

The final step of our preconditioning recipe was to split the light curves from each sector into smaller chunks that were individually analyzed later to become single points in the O–C diagram. Working on continuous segments, about 10 days has been shown to be a reasonable compromise between the need of getting 1) a frequency resolution that is high enough to reliably measure the phases of the individual pulsation modes and 2) a time resolution on the O–C diagram as low as possible (Murphy et al. 2014). Because for most TESS sectors the only significant interruption is the one-day central gap to allow the data downlink toward Earth at each perigee of the spacecraft, we chose that gap as a natural boundary and split each sector accordingly into two mostly continuous segments of about 14 days each, to which we refer hereafter as “orbits.” We are therefore left with 18 orbits for Chang 134 and 28 orbits for V393 Car. Each orbit was individually normalized to unit flux (and magnitude zero) before we continued to the next stages.

3.2 Frequency identification

In order to identify the most significant pulsation modes of our targets we made use of the Generalized Lomb-Scargle periodogram (GLS, Zechmeister & Kürster 2018), an improved version of the Lomb-Scargle periodogram (LS, VanderPlas 2018) which is less sensitive to aliasing and returns more accurate frequencies. Because we searched for LTTE effects, all the frequencies might in principle be slightly varying as a function of time due to the phase shift itself. At this stage, we therefore need to calculate the average value over the full span of our observations for each frequency. This can be accomplished by stitching all the normalized orbits into a single light curve and running GLS on it. We limited our search from 1 to 50 cycles per day, thus including the typical frequency range of δ Scuti modes (Balona & Dziembowski 2011) and being well within the Nyquist limit for our sampling rate (~360 cycles per day). The output are two periodograms for Chang 134 and V393 Car that are visually similar to those shown in Figs. 1 and 2, but with a much higher S/N. In order to improve the estimate of the pulsation frequencies, the light curve is usually fitted by a multisinusoidal function using the periodogram peaks as a starting point (Silvotti et al. 2018). In our case, however, this additional step is not essential given the extremely high S/N of our data.

By identifying the strongest frequencies v_i (and corresponding periods P_i) as the most prominent peaks in the periodogram, we found one dominant frequency for Chang 124 that was followed by its multiple integer harmonics, sorted hereafter by decreasing GLS power, $$\begin{array}{cc}\{\begin{array}{c}{P}_0=0.12942447\text{\hspace{0.17em}days}\\ P_1=0.06471223\text{\hspace{0.17em}days}\\ P_2=0.04314149\text{\hspace{0.17em}days}\end{array}& \{\begin{array}{c}{v}_0=7.72651415{\text{day}}^{-1}\\ v_1=15.4530283{\text{day}}^{-1}\\ v_2=23.1795424{\text{day}}^{-1}\end{array},\end{array}$$(3)

and in the same way for V393 Car, $$\begin{array}{cc}\{\begin{array}{c}{P}_0=0.14129519\text{\hspace{0.17em}days}\\ P_1=0.07064759\text{\hspace{0.17em}days}\\ P_2=0.04709840\text{\hspace{0.17em}days}\end{array}& \{\begin{array}{c}{v}_0=7.07738170{\text{day}}^{-1}\\ v_1=14.1547634{\text{day}}^{-1}\\ v_2=21.2321451{\text{day}}^{-1}\end{array}.\end{array}$$(4)

A complete list of all the significant frequencies detected for Chang 134 and V393 Car can be found in Tables B.1 and B.2, respectively.

For both stars, we obtained v₁ = 2v₀ and v₂ = 3v₀, that is, the three most prominent frequencies are the second and third harmonics of the main pulsation mode. Chang 134 does not show any significant peak outside the harmonic series of v₀. The GLS power of the strongest peak for V393 Car outside the harmonics of v₀ (at ~ 12.58288 day⁻¹) is three orders of magnitude smaller than the main mode. While irrelevant for our analysis due to its weakness, we note that this minor peak, detected at very high confidence in the TESS data, confirms the presence of an additional pulsation mode at 12.58 or 13.58 cycles per day claimed by Helt (1984) that was never confirmed by subsequent works (Garcia et al. 2001; Axelsen 2014). The very clean window function of the TESS data enables us to state that 13.58 was just a one-day alias of the real peak at 12.58 cycles per day.

3.3 Harmonic fit

After the average frequencies v_i = 2π/P_i of the most powerful pulsation modes were identified, we proceeded to fit a harmonic series in the form of Eq. (2) to our light curves. As several TESS orbits are affected by long-term systematic errors even after the PDCSAP treatment (on a timescale of typically days to weeks), we added a polynomial term as a function of time to our model to mitigate their impact. After some trial and error, we concluded that a polynomial of order 2 (i.e., a quadratic baseline) is enough to provide us with an accurate fit. The full model M we fit is therefore $$M=b_0+b_1\left(t-t_{\text{m}}\right)+b_2{\left(t-t_{\text{m}}\right)}^2+{\displaystyle \sum _{i=0}^{N-1}{f}_i\left(A_i,P_i,{\varphi }_i\right),}$$(5)

where t_m is the median time of each orbit, set to minimize the correlations between the fit b_i parameters. The P_i periods must be fixed to their average values determined in the previous section in order to obtain meaningful phase shifts; for N harmonic components included, the number of free parameters is therefore 3 + 2N. While N might in principle be arbitrarily large, for instance, when all the components that are even barely detectable in the frequency spectrum are included, almost all the information on the phase shift is contained in the three most significant peaks, which sum up to ≫99% of the total power. We therefore adopt N = 3 for our model by including the Pi constants from Eqs. (3) and (4), implying nine free parameters for our fit. We verified this assumption by adding further harmonic terms as a test, which resulted in statistically indistinguishable O–C diagrams at the expense of a much more intensive computation.

We fit the model in Eq. (5) through a differential evolution Monte Carlo Markov chain algorithm (DE-MC; Ter Braak 2006), in which multiple MC chains are run in parallel and learn from each other how to converge to a global minimum in the parameter space, rather than running independently as in the classical approach. With respect to the least-squares techniques used by most previous works to fit harmonic series to light curves, DE-MC allowed us to obtain much more reliable error estimates on the final best-fit parameters because the posterior distribution takes any correlation between them into account. In order to increase the efficiency, we first fit only the oscillating term ∑_if_i of our model to each orbit with a Levenberg–Marquardt method (Moré 1978), to obtain reasonably good starting points for A_i and ϕ_i. Then we ran the DE-MC code, again orbit by orbit, to fit the full model with uninformative priors centered on the A_i, ϕ_i values found above, while the b_i parameters were initially set to zero. After 1 000 000 steps, the first 10% of the chain was discarded as burn-in phase and the posterior distributions of the nine fit parameters were examined to ensure that convergence was reached.

At this stage, when we plotted the residuals of the light curves after the best-fit model had been subtracted, we realized that a non-negligible number of outliers had escaped the previous filtering steps. These data points are mostly located close to momentum dumps of the spacecraft or at the beginning of each orbit when TESS cameras are not yet thermally stabilized. To deal with them, we discarded all the outliers at more than 5σ and fed the clipped light curves for a second run of the DE-MC fit, with the same configuration. The final best-fit parameters, especially the phases ϕ_i to which we are most interested in, are statistically consistent with those delivered from the first run, but usually with slightly smaller uncertainties.

Finally, the phase shifts ϕ_i, by our definition within the [0,1] interval, were converted into O–C delays through a P_i factor. In principle, if the phase shifts were associated with independent pulsation modes and due solely to an LTTE, all the ϕ_i computed for the same orbit should result in the same O–C value, within the measurement errors: ϕ₀P₀ = ϕ₁P₁ = ϕ₂P₂, and computing a weighted average 〈ϕ_iP_i〉 would be the rigorous way to obtain the overall phase shift. This would also be a great opportunity to independently verify whether the LTTE model is the best explanation because different modes should yield the same O–C (Hermes 2018). In our case, however, all the ϕ_i come from a single pulsation mode, and the harmonics do not carry additional information with respect to the fundamental frequency v₀, P₀, where most of the spectral power lies: σ(ϕ₀) ≪ σ(ϕ₁) ≪ σ(ϕ₂) and thus ϕ₀P₀ ≃ 〈ϕ_iP_i〉. For this reason, after confirming that the assumption above is true, we defined ϕ₀P₀ as our O–C. We associated the median epoch of the orbit t_m as a time stamp for each value as introduced above: (O–C)(t_m) = ϕ₀ × P₀. The resulting O–C diagrams for Chang 134 and V393 Car are plotted in Figs. 3 and 4. The formal error bars are 1.5 s on average (range: 1.1–2.6 s) and 0.16 s (range: 0.14–0.19 s), respectively.

Fig. 3 Orbital solution for Chang 134. Upper panel: O–C diagram for the main pulsation mode of Chang 134 in seconds. The best-fit LTTE solution (including a quadratic trend) is overplotted as a continuous cyan line, and the quadratic baseline is plotted as a dashed orange line (see Sect. 3 for details). Lower panel: residuals from the best-fit model in seconds.

Fig. 4 Quadratic orbital solution for V393 Car. Upper panel: O–C diagram for the main pulsation mode of V393 Car in seconds. The best-fit LTTE solution with a quadratic baseline (O–C)quad (Eq. (7)) is overplotted as a continuous cyan line (see Sect. 3 for details). The quadratic baseline a0 + a1t + a2t2 is plotted separately as a dashed orange line. The nominal error bars (~0.16 s on average) are smaller than the point size. Lower panel: residuals from the best-fit model in seconds.

4 LTTE modeling

Even a quick look at both our O–C diagrams (Figs. 3 and 4) reveals that neither can be reasonably fit by a straight line, that is, the resulting reduced ${\chi }_r^2$ is ≫ 1. The presence of a nonzero second derivative implies that the phase of the signal is evolving with time, and we can interpret this modulation by assuming that it is due to an LTTE from unseen companions. Following Kepler et al. (1991), we parameterize our model as $$\text{O–C}=\left(t_0+\text{Δ}P·t\right)+\frac{1}{2}\frac{\dot{P}}{P}{t}^2+\alpha \text{cos}\left(\frac{2\pi t}{P_{\text{orb}}}+\phi \right),$$(6)

where t₀ is an arbitrary reference time, ∆P is the difference between the actual average pulsation period P and that estimated in Eqs. (3) and (4), $\dot{P}=\text{d}P\text{/d}t$ accounts for a linear variation of P as a function of time, and a is the amplitude of an LTTE signal from a perturber on a circular⁶ orbit having an orbital period P_orb and with a phase φ.

From a physical point of view, the first term of Eq. (6) does not carry information because it represents just a change in the slope of the O–C and does not imply any variation of the pulsation period. The second term, in the LTTE framework, represents a constant acceleration with respect to the barycenter of the system, possibly caused by a massive, perturbing body whose P_orb is much longer than our observing baseline (≈1000 days), such as stellar companions on very wide orbits. We should be aware here that several physical mechanisms other than LTTE can result in a P ≠ 0 term, including nonlinear interactions between different pulsation modes and stellar evolution effects. We return on this point in Sect. 4.2.

Last, the third, oscillating term of Eq. (6) is the LTTE modulation we are most interested in. As the first two terms play the role of a quadratic baseline in our fit, we make it more explicit by rewriting the model with a₀ = t₀, a₁ = ∆P and $a_2=0.5\left(\dot{P}\text{/}P\right)$, $${\left(\text{O–C}\right)}_{\text{quad}}=a_0+a_1\left(t-t_{\text{m}}\right)+a_2{\left(t-t_{\text{m}}\right)}^2+\alpha \cos \left(\frac{2\pi t}{P_{\text{orb}}}+\phi \right),$$(7)

where the median epoch of each orbit t_m has been subtracted from the time variable in order to minimize the correlations between a₀, a₁, and a₂.

In order to confirm the uniqueness of our best-fit solution and the significance of the quadratic coefficient a₂, we also performed a second fit to each data set by fixing a₂ = 0 in our model, that is, by imposing a linear baseline, $${\left(\text{O–C}\right)}_{\text{lin}}=a_0+a_1\left(t-t_{\text{m}}\right)+\alpha \cos \left(\frac{2\pi t}{P_{\text{orb}}}+\phi \right).$$(8)

According to the Irwin (1952) model, for a circular orbit (eccentricity e = 0) and assuming that the perturbing body is much less massive than the star (m_p ≪ M_*), the amplitude of the LTTE signal is $$\alpha =\frac{a\sin \left(i\right)}{c}\frac{m_{\text{p}}}{M_{*}},$$(9)

where c is the speed of light, a is the orbital semimajor axis of the perturber and i is the inclination of its orbital plane with respect to the sky plane. In other words, when a and P_orb are measured from the O–C data, the stellar mass is known (from Table 1) and a is derived from Kepler’s laws, we can constrain the minimum mass of the perturber as m_p sin i, with the same degeneracy on i as for the RV and astrometric techniques, $$m\sin \left(i\right)=a\times c{\left(\frac{M_{*}}{P_{\text{orb}}}\right)}^{2/3}{\left(\frac{G}{4{\pi }^2}\right)}^{-1/3}.$$(10)

We started the fitting process by computing a GLS periodogram on our O–C diagrams to search for periodic modulations. For both targets, a single, well-defined peak stands out, at about 82 days and 1000 days for Chang 134 and V393 Car, respectively. We then fit our full LTTE model with both a quadratic and linear baseline (parameterized as in Eqs. (7) and (8), respectively) on our O–C diagrams through DE-MC, setting uninformative priors on all the five or six free parameters involved: a₀, a₁, a₂, α, P_orb, and φ in the (O–C)_quad model, and α₀, α₁, α, P_orb, and φ in the (O–C)_lin model, and centering the boundaries of P_orb on the values previously found by GLS to speed up the convergence. As done at the harmonic analysis stage, after 1000 000 steps, the first 10% of the chain was then discarded as burn-in phase, and we extracted the best-fit values and the associated errors of our parameters from their posterior distributions. All these values are reported in Table 2 for the four independent fits. The mass of the perturbing body m_p sin i was calculated as well as a derived parameter by propagating the relevant errors in Eq. (10).

4.1 Orbital solution for Chang 134

The best-fit (O–C)_quad model for Chang 134 is overplotted on our O–C in the upper panel of Fig. 3 as a solid cyan line, with the quadratic baseline a₀ + a₁t + a₂t² plotted as a dashed orange line. The resulting χ² of the residuals (lower panel of Fig. 3) is 56.7, with the reduced χ² being ${\chi }_{\text{red}}^2={\chi }^2\text{/DOF}\simeq 4.36$, where DOF = 18 − 6 = 12 are the degrees of freedom of our fit. It is worth noting that the best-fit value for a₂ is statistically consistent with zero (Table 2, second column), that is, we do not see any evidence of long-term LTTE effects because the derived $\dot{P}$ is also consistent with zero. This is confirmed by the fact that the best-fit linear model (O – C)_lin (Table 2, first column) yielded the same results within the error bars for all the fit parameters, with a virtually indistinguishable χ² (56.4 vs. 56.7). In other words, because the (O–C)_quad model is not favored by our data, we hereafter adopt the results from the (O–C)_lin fit.

The best-fit returned the amplitude of the LTTE term to be α = 6.11 ± 0.65 s (according to Eq. (6)), and after propagating the errors through Eq. (9), we obtain a minimum mass of m_p sin i = 43.1 ± 4.7 M_jup for the hypothetical companion. This is within the brown dwarf regime (13–80 M_jup; Grieves et al. 2021).

4.2 Orbital solution for V393 Car

Unlike the case of Chang 134, the results from V393 Car fits are more difficult to interpret because the two models (O–C)_quad and (O–C)_lin, while having a similar shape, led us to different orbital solutions. The two best-fit models are overplotted on our O–C in the upper panel of Figs. 4 and 5 as a solid cyan line, respectively, with the linear and quadratic baseline a₀ + a₁t + a₂t² plotted as a dashed orange line.

While it is clear that at least an oscillating term is needed to fit the O–C data, the inclusion of a quadratic term yields different parameters for the LTTE sinusoid: P ≃ 1071 days with amplitude α ≃ 107 s for the (O–C)_lin fit versus P ≃ 723 days with amplitude α ≃ 48 s for the (O–C)_quad fit. The discrepancy of a factor of ~2 in the α parameter translates into a similar ratio for the derived hypothetical companion masses: m_p sin i = 175.7 ± 9.2 M_jup, or 0.17 ± 0.01 M_⊙ (linear model) and m_p sin i = 102.0 ± 5.3 M_jup, or 0.10 ± 0.01 M_⊙ (quadratic model). We emphasize that for both scenarios, the minimum mass would be consistent with a cool dwarf of spectral type M5V and M6V, respectively (Pecaut & Mamajek 2013).

A technique that is commonly employed for model selection problems like this is the evaluation of metrics such as the Akaike information criterion (AIC; Akaike 1974) or the Bayesian information criterion (BIC; Schwarz 1978). In the case of Gaussian distributions, the latter can be written as BIC = χ² + k ln n, where k is the number of free parameters and n is the number of data points. If applied literally to our problem, the BIC would overwhelmingly favor our (O–C)_quad scenario (∆BIC ≫ 10). Unfortunately, this result is distorted by the fact that the χ² of the residuals for both fits (lower panel of Fig. 4 and 5) is exceptionally high due to the extremely low formal error bars with respect to the systematic errors at play (on which we comment in Sect. 5): for the (O–C)_quad fit, for instance, we obtained χ² = 8627, with the reduced χ² being ${\chi }_{\text{red}}^2={\chi }^2\text{/DOF}\simeq 392$, where DOF = 28 − 6 = 22 are the degrees of freedom of our fit. This fact, and the very sparse phase sampling of our signal, prevents us from reaching a firm conclusion on the properties of this candidate companion, which will require further data to be confirmed and better constrained. New TESS observations of V393 Car are planned for Cycle 5, during four consecutive sectors (61–64; from January to May 2023), and again in Sector 68 (August 2023). Since the predictions from the (O–C)_quad and (O–C)_lin models will diverge by ~12 h by mid-2023, we expect that the extension of our analysis to the new TESS light curves will definitely be conclusive about the correct scenario.

As a side note, we could wonder whether the quadratic term we found when we assumed a quadratic baseline (a₂ = 650 ± 4 ppm s⁻¹ as fitted on our O–C plane, translating into $\dot{P}\text{/}P\simeq 5.5\times {10}^{-6}$ yr⁻¹) might be consistent with other non-LTTE mechanisms, as previously mentioned. Significant nonlinear interactions between pulsating modes (Silvotti et al. 2018; Bowman et al. 2021) are not at play here because the spectral power outside the main pulsation mode and its harmonics is essentially negligible. On the other hand, the typical $\dot{P}\text{/}P$ expected from evolutionary effects range from 10⁻⁹ to 10⁻⁷ yr⁻¹ (Breger & Pamyatnykh 1998; Xue et al. 2022), that is, it is much lower than our fitted value.

The same reasoning can be applied to Chang 134, for which we measure a nonsignificant quadratic term a₂ = −2.2 ± 7.5 ppm s⁻¹ corresponding to $\dot{P}\text{/}P\simeq \left(-1.8\pm 6.6\right)\times {10}^{-8}$ yr⁻¹.

Unlike V393 Car, this upper limit is still within the predictions for evolutionary effects. For instance, if we adopt $\dot{P}\text{/}P\simeq 3\times {10}^{-9}$ yr⁻¹ (as measured on AE UMa by Xue et al. 2022), we obtain a total contribution on the O–C of ~0.35 s over the full baseline of our observations (~1000 days), which is well below our current measurement error.

Fig. 5 Linear orbital solution for V393 Car. Upper panel: same as Fig. 4, but with the best-fit LTTE solution with a linear baseline (O–C)lin (Eq. (8)). The linear baseline a0 + a1t is plotted separately as a dashed orange line. Lower panel: residuals from the best-fit model in seconds.

5 Discussion and conclusions

We applied the PT technique to two δ Scuti stars observed by TESS in a large number of sectors, detecting a periodic modulation on their O–C that is consistent with the presence of companions in both cases: a BD (m_p sin i = 43.1 ± 4.7 M_jup) on a P ≃ 82 days orbit around Chang 134, and a more massive body (m_p sin i ≃ 0.10–0.17 M_⊙) around V393 Car, whose orbital parameters are still not firmly constrained by the limited phase sampling of the signal. This double detection on the very first two targets analyzed by our project may appear to be very lucky. Still, the binary fraction of δ Scuti primaries as estimated by a previous PT search on Kepler data is rather high: 15.4 ± 1.4% (Murphy et al. 2018). In addition, unlike the quoted work, we initially restricted our sample in order to target only the most favorable stars in terms of sensitivity to LTTE signals. The overall binary fraction of A-type stars is estimated to be significantly larger than 50% by most authors (Duchêne & Kraus 2013), making our results less surprising. Any further statistical implication is prevented by the small size of our sample and by the heavy selection effects at play. As a side note, however, we mention that Borgniet et al. (2019) constrained the BD frequency within 2–3 au around AF stars to be below 4% (1σ) based on 225 targets observed with SOPHIE and/or HARPS, and that other authors hypothesized the presence of a brown dwarf desert extending to early-type stars (Murphy et al. 2018). This would make the companion of Chang 134 an uncommon object, worthwhile to be followed up with more data and/or different techniques. The future TESS observations (at least five sectors of Cycle 5 are already scheduled, starting in March 2023) will also possibly help in constraining the eccentricity of the LTTE orbit, forced to zero in our analysis to avoid overfitting given our limited phase coverage and small number of O–C points.

In addition to our results for these two specific stars, our analysis was based on an independent and improved implementation of the PT method (Murphy et al. 2014). It is also intended as a pilot study that is a prelude to a systematic search of LTTE companions around pulsating stars on a much larger scale, exploiting the huge sample of δ Scuti targets that are and will be monitored through the availability of TESS FFIs at 10-min cadence. The results for Chang 134, in particular, demonstrates that the detection of bodies in the substellar regime is perfectly feasible from TESS light curves, and that in some particularly favorable configuration, candidates with planetary masses (m_p sin i ≲ 13 M_jup) might be already within reach. A perturber with the same LTTE amplitude as was found on Chang 134 (α ≃ 7 s), for instance, but on a larger orbit at P ≃ 800 days, would imply a companion with a mass of ≃2 M_jup detected at high significance.

On the other hand, our analysis highlighted several limiting factors that are worth discussing further. The first factor, largely self-evident in the case of V393 Car (Sect. 4.2), is related to the sparse sampling of the TESS light curves, an unavoidable consequence of its fixed scanning law (Ricker et al. 2015). Even for targets within or close to the CVZ at ecliptic latitudes |β| ≳ 78°, that is, targets that are observed in up to 13 contiguous sectors, there is always a large gap every other year, which make the detection of signals at longer periods (P ≳ 1 yr) difficult or ambiguous. Unfortunately, this is also the region in the parameter space where the PT technique is more sensitive to low-mass perturbers, owing to the Eq. (9) relation. As TESS is continuing its mission through extension phases, the accumulated data will mitigate this issue by filling out all or most of the orbital phases on the folded O–C diagram.

A second, subtler limiting factor is the excess scatter in our O–C data, as can be seen from the very high ${\chi }_{\text{r}}^2$ values reported in Table 2 and discussed Sect. 4.2. Our residuals from the best-fit LTTE model have an rms of ~3 s in both fits, regardless of the average nominal error bars, which are about ten times lower on V393 Car than on Chang 134. The most plausible cause for this are systematic errors in the absolute calibration of the time stamp. This issue has been anticipated by the TESS team, and in particular by the TESS Asteroseismic Science Consortium (TASC), who gave accurate timing requirements in the SAC/TESS/0002/6 document⁷, including RS-TASC-05: “[…] the time given for each exposure should be accurate over a period of 10 days to better than 1 second.” No requirement is given at longer timescales, however, and the TESS data release notes (DRN⁸) have reported systematic offsets in the time stamps for several sectors on the order of seconds, that is, consistent with the excess O–C scatter we measured. While some of the reported offsets have been already corrected in subsequent update data releases, others are still waiting to be implemented. A ground-based campaign by von Essen et al. (2020) tried to independently confirm the absolute time calibration of TESS by comparing observations of a sample of eclipsing binaries in common, measuring a global offset of 5.8 ± 2.5 s; the errors of the measurements at individual epochs, however, are too large (≳ 10 s) to confirm or disprove the systematic errors we see. We emphasize that offsets of a few seconds are usually completely negligible when investigating transiting exoplanets, but they are crucial in extending the sensitivity of the PT technique to planetmass companions. An important outcome of a large-scale PT analysis of TESS data will allow us also to identify and correct these systematic errors by comparing the timing residuals of tens or hundreds of targets, following a self-calibration approach.

The same technique as applied in this work, and in general, the expertise gained through this project, will also help the scientific exploitation of PLATO (Rauer et al. 2014), for which the PT analysis could be a compelling case of ancillary science. Unlike TESS, no FFIs will be downloaded from PLATO during the nominal observing phase, and no stars earlier than F5 are to be included in the main target samples, which are focused on transit search around solar-type stars (Montalto et al. 2021). Nevertheless, about 8% of the PLATO science data rate will be allocated to the General Observer (GO) program through ESA calls open to the whole astrophysical community. Any δ Scuti star (and, more in general, pulsating star) allocated as GO within a long-pointing field (Nascimbeni et al. 2022) would result in a mostly uninterrupted light curve on a baseline of 2–4 yr, with a cadence⁹ and photometric precision far better than what is achieved by TESS, overcoming the sparse phase sampling affecting the latter. Even more important, the time stamps of each PLATO data point will be accurate within 1 s in an absolute sense by formal scientific requirement because of the specifically designed calibration processes that also include the monitoring of a preselected set of detached eclipsing binaries (Nascimbeni et al. 2022).

Acknowledgements

The authors wish to thank the referee, Dr. Roberto Silvotti, for the valuable and fruitful comments and suggestions which significantly improved our manuscript. This research has made use of the SIMBAD database (operated at CDS, Strasbourg, France;Wenger et al. 2000), the VAR-TOOLS Light Curve Analysis Program (version 1.39 released October 30, 2020, Hartman & Bakos 2016), TOPCAT and STILTS (Taylor 2005, 2006). Valerio Nascimbeni and Giampaolo Piotto recognize support by ASI under program PLATOASI/INAF agreements 2015-019-R.1-2018.

Note added in proof

After the acceptance of this paper, we were contacted by Dr. D. Hey, who let us know that he found a similar orbital solution for V393 Car using the maelstrom code (Hey et al. 2020). In particular, his results suggest that the linear-baseline model yields a better fit the timing data.

Appendix A O–C data tables

Appendix B Detected frequencies

References

All Tables

All Figures

	Fig. 1 TESS photometry of Chang 134 = TYC 9158-919-1. Upper panel: one-day section of the light curve from TESS sector 1. Lower panel: GLS periodogram of the whole light curve including nine TESS sectors, stitched together and filtered as explained in Sect. 3.
In the text

	Fig. 2 TESS photometry of V393 Carinae = HD 66260. Upper panel: one-day section of the light curve from TESS sector 1. Lower panel: GLS periodogram of the whole light curve including 14 TESS sectors, stitched together and filtered as explained in Sect. 3.
In the text

	Fig. 3 Orbital solution for Chang 134. Upper panel: O–C diagram for the main pulsation mode of Chang 134 in seconds. The best-fit LTTE solution (including a quadratic trend) is overplotted as a continuous cyan line, and the quadratic baseline is plotted as a dashed orange line (see Sect. 3 for details). Lower panel: residuals from the best-fit model in seconds.
In the text

Fig. 4 Quadratic orbital solution for V393 Car. Upper panel: O–C diagram for the main pulsation mode of V393 Car in seconds. The best-fit LTTE solution with a quadratic baseline (O–C)quad (Eq. (7)) is overplotted as a continuous cyan line (see Sect. 3 for details). The quadratic baseline a0 + a1t + a2t2 is plotted separately as a dashed orange line. The nominal error bars (~0.16 s on average) are smaller than the point size. Lower panel: residuals from the best-fit model in seconds.

In the text

	Fig. 5 Linear orbital solution for V393 Car. Upper panel: same as Fig. 4, but with the best-fit LTTE solution with a linear baseline (O–C)lin (Eq. (8)). The linear baseline a0 + a1t is plotted separately as a dashed orange line. Lower panel: residuals from the best-fit model in seconds.
In the text