11institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117
Heidelberg, Germany 22institutetext: Leiden Observatory, Leiden University,
Postbus 9513, 2300 RA Leiden, The Netherlands 33institutetext: Institut für
Astro- und Teilchenphysik, Universität Innsbruck, Technikerstraße 25, A-6020
Innsbruck 44institutetext: Institute for Astronomy, University of Hawai‘i,
Honolulu, HI 96822, USA 55institutetext: Centre for Astrophysics, University
of Southern Queensland, Toowoomba, QLD 4350, Australia 66institutetext:
Institut für Physik, Karl-Franzens Universität Graz, Universitätsplatz 5/II,
NAWI Graz, 8010 Graz, Austria 77institutetext: Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, France
88institutetext: Italian National Agency for New Technologies, Energy and
Sustainable Economic Development (ENEA), via Anguillarese 301, Rome, Italy.
99institutetext: Institut polaire français Paul Émile Victor and Programma
Nazionale di Ricerche in Antartide, Concordia Station, Antarctica
1010institutetext: Istituto Nazionale di Astrofisica, Osservatorio Astronomico
di Cagliari, Italy 1111institutetext: Department of Physics & Astronomy,
University of Rochester, Rochester, NY 14627, USA 1212institutetext:
Department of Physics, University of California at Santa Barbara, Santa
Barbara, CA 93106, USA 1313institutetext: Department of Physics and Astronomy,
Astronomy and Space Physics, Uppsala University, 751 20 Uppsala, Sweden
1414institutetext: Institut de Recherche sur les Exoplanètes, Département de
Physique, Université de Montréal, Montréal, QC H3C 3J7, Canada
1515institutetext: unaffiliated 1616institutetext: Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Drive, M/S321-100,
Pasadena, CA 91109, USA 1717institutetext: LESIA, Observatoire de Paris,
Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules
Janssen, 92195 Meudon, France 1818institutetext: Research School of Astronomy
and Astrophysics, Australian National University, Canberra, ACT 2611,
Australia 1919institutetext: South African Astronomical Observatory,
Observatory Rd, Observatory Cape Town, 7700 Cape Town, South Africa
2020institutetext: Department of Astronomy, University of Cape Town,
Rondebosch, 7700 Cape Town, South Africa 2121institutetext: Astronomy
Department, University of California, Berkeley, CA 94720, USA
2222institutetext: SETI Institute, Carl Sagan Center, 189 Bernardo Ave.,
Mountain View CA 94043, USA 2323institutetext: Institute of Astrophysics,
FORTH, GR-71110 Heraklion, Greece 2424institutetext: Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA) and
Department of Physics and Astronomy, Northwestern University, Evanston, IL
60208, USA 2525institutetext: Johns Hopkins APL, Laurel, MD, USA
2626institutetext: Astrophysics Research Centre, Queen’s University Belfast,
Belfast BT7 1NN, UK 2727institutetext: School of Physical Sciences and Centre
for Astrophysics & Relativity, Dublin City University, Glasnevin, Dublin 9,
Ireland 2828institutetext: IPAG, Univ. Grenoble Alpes, CNRS, IPAG, F-38000
Grenoble, France 2929institutetext: IMCCE - Observatoire de Paris, 77 Avenue
Denfert-Rochereau, F-75014 PARIS 3030institutetext: Institut d’Astrophysique
de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98 bis boulevard
Arago, 75014 Paris, France 3131institutetext: Institute of Astronomy,
Madingley Road, Cambridge CB3 0HA, UK 3232institutetext: Department of
Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK
3333institutetext: Centre for Exoplanets and Habitability, University of
Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK 3434institutetext: Shanghai
Observatory, Chinese Academy of Sciences, PR China 3535institutetext: Purple
Mountain Observatory, Chinese Academy of Science, Nanjing 210008, PR China

# The $\beta$ Pictoris b Hill sphere transit campaign

II. Searching for the signatures of the $\beta$ Pictoris exoplanets
through time delay analysis of the $\delta$ Scuti pulsations

Sebastian Zieba 112233 Konstanze Zwintz 33 Matthew Kenworthy 22 Daniel Hey 44
Simon J. Murphy 55 Rainer Kuschnig 66 Lyu Abe 77 Abdelkrim Agabi 77 Djamel
Mekarnia 77 Tristan Guillot 77 François-Xavier Schmider 77 Philippe Stee 77
Yuri De Pra 88 Marco Buttu 991010 Nicolas Crouzet 22 Samuel Mellon 1111 Jeb
Bailey III 1212 Remko Stuik 22 Patrick Dorval 1313 Geert-Jan J. Talens 1414
Steven Crawford 1515 Eric Mamajek 11111616 Iva Laginja 1717 Michael Ireland
1818 Blaine Lomberg 19192020 Rudi Kuhn 1919 Ignas Snellen 22 Paul Kalas
222223232424 Jason J. Wang 2525 Kevin B. Stevenson 2626 Ernst de Mooij
27272828 Anne-Marie Lagrange 292930303131 Sylvestre Lacour 1717 Mathias Nowak
3131 Paul A. Strøm 32323333 Zhang Hui 3434 Lifan Wang 3535

(Received August 18, 2023; accepted March 12, 2024)

The $\beta$ Pictoris system is the closest known stellar system with directly
detected gas giant planets, an edge-on circumstellar disc, and evidence of
falling sublimating bodies and transiting exocomets. The inner planet, $\beta$
Pictoris c, has also been indirectly detected with radial velocity (RV)
measurements. The star is a known $\delta$ Scuti pulsator, and the long-term
stability of these pulsations opens up the possibility of indirectly detecting
the gas giant planets through time delays of the pulsations due to a varying
light travel time. We search for phase shifts in the $\delta$ Scuti pulsations
consistent with the known planets $\beta$ Pictoris b and c and carry out an
analysis of the stellar pulsations of $\beta$ Pictoris over a multi-year
timescale. We used photometric data collected by the BRITE-Constellation,
bRing, ASTEP, and TESS to derive a list of the strongest and most significant
$\delta$ Scuti pulsations. We carried out an analysis with the open-source
python package maelstrom to study the stability of the pulsation modes of
$\beta$ Pictoris in order to determine the long-term trends in the observed
pulsations. We did not detect the expected signal for $\beta$ Pictoris b or
$\beta$ Pictoris c. The expected time delay is 6 seconds for $\beta$ Pictoris
c and 24 seconds for $\beta$ Pictoris b. With simulations, we determined that
the photometric noise in all the combined data sets cannot reach the
sensitivity needed to detect the expected timing drifts. An analysis of the
pulsational modes of $\beta$ Pictoris using maelstrom showed that the modes
themselves drift on the timescale of a year, fundamentally limiting our
ability to detect exoplanets around $\beta$ Pictoris via pulsation timing.

###### Key Words.:

Asteroseismology

##  1 Introduction

$\beta$ Pictoris is a nearby southern hemisphere star visible with the naked
eye for which $\delta$ Scuti -like pulsations were discovered by Koen (2003).
The planetary-mass companion $\beta$ Pictoris b was detected using the
VLT/NaCo instrument with direct imaging (Lagrange et al., 2009b, 2010).
Evidence of a second planet in the $\beta$ Pictoris system was published by
Lagrange et al. (2019b) using the radial velocity (RV) method and recently
independently confirmed by Nowak et al. (2020) and Lagrange et al. (2020)
using VLTI/GRAVITY observations.

The lifetime and frequency stability of $\delta$ Scuti pulsations make them
astronomical “stellar clocks,” and therefore, they are great targets for
applying timing techniques (Compton et al., 2016). The common orbital motion
of a star together with a companion around the barycenter of a system results
in a periodic early or late arrival of the stellar pulsational signals. This
principle led to the first detection of planets orbiting a pulsar outside the
Solar System (Wolszczan & Frail, 1992; Wolszczan, 1994). The periodic
variation of the arrival times can be seen as either a frequency modulation
(FM; Shibahashi & Kurtz, 2012; Shibahashi et al., 2015) or a phase modulation
(PM; Murphy et al., 2014; Murphy & Shibahashi, 2015; Murphy et al., 2016b).
The latter perspective works better for companions in wider orbits.

By applying the PM method on Kepler data, Murphy et al. (2016a) discovered a
massive planet ($m\sin i\approx 12M_{J}$) with an orbital period of about 840
days around a $\delta$ Scuti star. In addition to the discovery of this
planet, the PM method has led to the detection of 341 binaries and hundreds of
more candidates (Murphy et al., 2018). Furthermore, it has provided us with
the eccentricity, period, and mass function of these companions orbiting
stars, just as the RV method does (e.g. Nesvold & Kuchner, 2015). Applying the
same method to pulsating stars observed by the TESS mission will lead to many
more binary systems with full orbital solutions.

In this work, we use the data collected by the TESS satellite in its primary
mission and data collected by the Hill sphere111The Hill sphere is the region
around a planet where masses, such as moons and planetary rings, are
gravitationally bound to the planet. transit campaign, which was an
international effort of space-based (e.g. through the BRITE-Constellation) and
ground-based (e.g. through bRing, ASTEP) observations that searched for
signatures of material around the giant planet $\beta$ Pictoris b (Kalas et
al., 2019; Kenworthy, 2017). We analyzed this photometric data by searching
for phase variations (and therefore time delays) caused by orbital motion in
the pulsational signals. $\beta$ Pictoris was observed for approximately four
months from October 2018 to February 2019 during the primary mission of TESS.
A second visit occurred during TESS’ extended mission from November 2020 to
February 2021 (see Table 3 for a summary of all visits).

In Section 2, we describe the properties of the different components in the
$\beta$ Pictoris system. Section 3 has a summary of all observational
instruments and a frequency analysis for the photometry collected by TESS. The
theory and equations for this paper can be found in Section 4, and our results
and conclusions follow in Sections 5 and 6.

##  2 The $\beta$ Pictoris system

$\beta$ Pictoris (HD 39060; HR 2020) is one of the most studied and intriguing
star-planet systems. The Infrared Astronomical Satellite (IRAS) discovered an
infrared excess (Aumann et al., 1984) for this bright and close southern star
that was attributed to the presence of a circumstellar disc. This disc was
first imaged by Smith & Terrile (1984) and clearly showed the edge-on geometry
of this system. The gas and dust in this disc is mostly “second generation,”
that is, constantly replenished by collisions of comets and asteroids
(Lagrange et al., 2000). A warp in this disc (Augereau et al., 2001; Mouillet
et al., 1997; Nesvold & Kuchner, 2015) and signatures of evaporating exocomets
(also called falling evaporating bodies, or FEBs) in spectroscopy (Ferlet et
al., 1987; Beust & Morbidelli, 2000) were attributed to an exoplanet orbiting
the star interacting dynamically with its environment (for more information on
exocomets observed around $\beta$ Pictoris detected in photometry, see Zieba
et al., 2019; Strøm et al., 2020; Pavlenko et al., 2022; Lecavelier des Etangs
et al., 2022).

###  2.1 The star

Koen (2003) discovered $\delta$ Scuti -type pulsations at the millimagnitude
level originating from $\beta$ Pictoris. Further analysis by Mékarnia et al.
(2017), Zwintz et al. (2019), and Zieba et al. (2019) showed dozens of
additional frequencies between 20 and 80 cycles per day. An asteroseismic
large spacing, $\Delta\nu$, has been measured for $\beta$ Pictoris (Bedding et
al., 2020), which might facilitate a precise asteroseismic age in the future.
The pulsations also induce intrinsic variations in the RV at $\lesssim$ 1
$\mathrm{km\,s}^{-1}$ peak to peak (Lagrange et al., 2009a, 2012; Galland et
al., 2006), which hampers the search for planets with the RV method in this
system. A selection of the fundamental properties of the star $\beta$ Pictoris
are listed in Table 1.

Table 1: Various stellar parameters of the star $\beta$ Pictoris.

Parameter | Value | Reference
---|---|---
RA (J2000.0) | 05h 47m 17.09s | 1
DEC (J2000.0) | -51h 03m 59.41s | 1
V (mag) | 3.86 | 2
TESS (mag) | 3.696 | 1
age (Myr) | 23 $\pm$ 3 | 3
parallax (mas) | 50.93 $\pm$ 0.15 | 4,5,6
Distance (pc) | 19.63 $\pm$ 0.06 | 4,5,6
Spectral class | A6V | 7
Radius ($R_{\odot}$) | 1.497 $\pm$ 0.025 | 8
Mass ($M_{\odot}$) | 1.75${}_{-0.02}^{+0.03}$ | 9
Teff (K) | 8090 $\pm$ 59 | 8

(1) Stassun et al. (2019); (2) Cousins (1971); (3) Mamajek & Bell (2014); (4)
Gaia Collaboration et al. (2016); (5) Gaia Collaboration et al. (2023); (6)
Lindegren et al. (2021); (7) Gray et al. (2006); (8) Zwintz et al. (2019); (9)
Lacour et al. (2021).

###  2.2 The planets: $\beta$ Pictoris b and c

The warp of the inner disc of $\beta$ Pictoris observed by the Hubble Space
Telescope and in ground-based observations was one of the indirect lines of
evidence for a massive substellar companion orbiting the star (Burrows et al.,
1995; Mouillet et al., 1997; Heap et al., 2000; Golimowski et al., 2006).
Signatures of infalling exocomets in the spectra of the star also needed a
“perturber” to scatter them onto eccentric inner system-bearing orbits. The
planet, $\beta$ Pictoris b, was then discovered using the VLT/NaCo instrument
data in 2003 (Lagrange et al., 2009b) and was later confirmed by Lagrange et
al. (2010). A transit-like event was observed in 1981 and attributed to a
planet (Lecavelier Des Etangs et al., 1995). However, a better orbit
determination with the VLT/SPHERE instrument ruled out $\beta$ Pictoris b as
the cause of that event (Lagrange et al., 2019a). Using data from the Gemini
Planet Imager, Wang et al. (2016) were able to rule out a transit of the
planet $\beta$ Pictoris b during the conjunction in 2017 at a 10$\sigma$
level. However, a Hill sphere transit was predicted for the time between late
2017 and early 2018 (Lecavelier des Etangs & Vidal-Madjar, 2016; Wang et al.,
2016). Various observational campaigns were initiated by PicSat (Nowak et al.,
2018; Menegaldo et al., 2022); bRing (Kenworthy, 2017); and the BRITE-
Constellation (Weiss et al., 2014) in order to photometrically observe
possible material around the planet; however, there was no significant
detection (Kenworthy et al., 2021). The mass and the orbital solution of the
outer planet $\beta$ Pictoris b are listed in Table 2.

Table 2: Various parameters of the planets $\beta$ Pictoris b and c based on
Lacour et al. (2021).

Parameter | Unit |  $\beta$ Pictoris b |  $\beta$ Pictoris c
---|---|---|---
Mass | $M_{J}$ | $11.90^{+2.93}_{-3.04}$ | $8.89^{+0.75}_{-0.75}$
a | au | $9.93^{+0.03}_{-0.03}$ | $2.68^{+0.02}_{-0.02}$
e | — | $0.103^{+0.003}_{-0.003}$ | $0.32^{+0.02}_{-0.02}$
i | ∘ | $89.00^{+0.00}_{-0.01}$ | $88.95^{+0.09}_{-0.10}$
$\varpi$(a)(a)footnotemark: $(a)$ | ∘ | $199.3^{+2.8}_{-3.1}$ | $66.0^{+1.8}_{-1.7}$
$\tau$(b)(b)footnotemark: $(b)$ | — | $0.719^{+0.008}_{-0.010}$ | $0.724^{+0.006}_{-0.006}$
P | years | $23.61^{+0.09}_{-0.09}$ | 3.34 $\pm$ 0.04
P | days | $8623^{+31}_{-32}$ | 1221 $\pm$ 15
$t_{p}$(c)(c)footnotemark: $(c)$ | MJD | 65243 | 59888
222
$a$$a$footnotetext: The argument of periastron $\varpi$ reported in Lacour et
al. (2021) follows the definition in Blunt et al. (2020) and therefore refers
to the orbit of the companion and not to the star. In this paper, however, we
refer to the orbit of the star when we use $\varpi$.
$b$$b$footnotetext: reference epoch MJD 59000 (31 May 2020).
$c$$c$footnotetext: derived from $\tau$.

Evidence of an additional planet in the $\beta$ Pictoris system was published
by Lagrange et al. (2019b). Over 6000 spectra of the star taken between 2003
and 2018 by the HARPS instrument at the ESO La Silla 3.6 m telescope have been
analyzed, and they showed a hint of a planetary signal. $\beta$ Pictoris c was
then ultimately directly detected by Nowak et al. (2020) and Lagrange et al.
(2020) using VLTI/GRAVITY observations. A list of parameters for the planet
can be found in Table 2.

##  3 Observations

Due to the 2017-2018 Hill Sphere Transit of $\beta$ Pictoris b, an
international campaign of space- and ground-based observations was launched in
order to search for signatures of material around the giant planet (Kalas et
al., 2019; Kenworthy et al., 2021). Table 3 summarizes various properties of
the different light curves. Changes to those light curves other than the
Gaussian high-pass procedure, which is explained in Section A.2, are noted in
the following subsections. The data provided by BRITE-Constellation was left
unchanged. A detailed analysis of the photometry of $\beta$ Pictoris collected
by BRITE-Constellation and bRing was published by Zwintz et al. (2019). As we
wanted to measure the periodic motion of a star around a barycenter, we also
wanted to correct for the motion the Earth in the Solar System. All the
observations in this work were therefore converted to the Barycentric Julian
Date in the barycentric dynamical time (BJD${}_{\text{TBD}}$) standard using
the Python tool BARYCORRPY (Kanodia & Wright, 2018), which is based on the IDL
code BARYCORR (Wright & Eastman, 2014).

Table 3: Summary of the properties of the various instruments and
corresponding light curves.

Observation |  | Wavelength
---
(nm)
| Observation
---
start
| Observation
---
end
| $T$
---
(days)
| $1/T$
---
($10^{-3}$ d-1)
| $f_{\text{Ny.}}$
---
( d-1)
| cadence
---
(s)
| DC
---
(%)
BHr | 550 - 700 | 16 March 2015 | 2 June 2015 | 78.32 | 12.77 | 4167 | 10.37 | 6.78
BTr + BHr | 550 - 700 | 4 Nov. 2016 | 17 June 2017 | 224.6 | 4.453 | 2128 | 20.30 | 7.07
BHr | 550 - 700 | 9 Nov. 2017 | 25 April 2018 | 167.3 | 5.976 | 2128 | 20.30 | 7.48
bRing | 463 - 639 | 2 Feb. 2017 | 1 Sept. 2018 | 575.5 | 1.738 | 135.4 | 319.1 | 27.0
ASTEP17 | 695 - 844 | 28 March 2017 | 14 Sept. 2017 | 170.0 | 5.881 | 495.8 | 87.13 | 18.9
ASTEP18 | 695 - 844 | 28 March 2018 | 15 July 2018 | 109.3 | 9.150 | 502.8 | 85.92 | 29.2
TESS | 600 - 1000 | 19 Oct. 2018 | 1 Feb. 2019 | 105.2 | 9.507 | 360.0 | 120.0 | 85.3
TESS | 600 - 1000 | 20 Nov. 2020 | 8 Feb. 2021 | 79.8 | 12.53 | 360.0 | 120.0 | 90.2

333
The term $T$ denotes the time base of the observations, the reciprocal value
$1/T$ corresponds to the Rayleigh criterion, $f_{\text{Ny.}}$ is the Nyquist
frequency, and DC is the duty cycle. BRITE Lem (BLb) is equipped with a blue
filter and observed $\beta$ Pictoris from December 2016 until June 2017, but
due to significantly higher noise in the time series, the data was disregarded
from the analysis. See Zwintz et al. (2019) for an analysis of the BLb
observations.

Figure 1: Full light curve of all available observations of the star $\beta$
Pictoris used in this work.

###  3.1 BRITE-Constellation

The BRITE-Constellation (Weiss et al., 2014) consists of five nanosatellites
collecting photometry for the brightest stars on the sky. In this work, we
analyzed data collected by three of the satellites: BRITE-Heweliusz (BHr),
BRITE-Toronto (BTr), and BRITE-Lem (BLb). Being in a low-earth orbit, the
orbital periods of the satellites are all around 100 minutes. A minimum of 15
minutes per orbit is dedicated to observations. Three different runs were
conducted in the constellations around Pictor and Vela, which also included
the star $\beta$ Pictoris. A summary of the durations and various properties
of the observations can be found in Table 3. The pipeline for the photometry
reduction is described in Popowicz et al. (2017). An analysis of all BRITE
observations was conducted in Zwintz et al. (2019). For the three runs by BHr,
BTr+BHr, and BHr, which all used the red BRITE filter, six, 13, and eight
significant frequencies were extracted, respectively. The only run with a blue
filter by BLb suffered from higher noise compared to the other BRITE
observations. Zwintz et al. (2019) has reported four frequencies in the
collected BLb photometry. The blue observations were discarded from this
analysis, as the data quality was not good enough to provide additional
information.

###  3.2 bRing

The bRing project (which stands for “the $\beta$ Pictoris b Ring project”) was
initiated in order to collect photometry of $\beta$ Pictoris during the Hill
sphere transit of $\beta$ Pictoris b at the end of 2017 (Stuik et al., 2017).
To that end, two stations in South Africa and Australia were built, each
consisting of two wide-field cameras. Their design is based on the Multi-Site
All-Sky CAmeRA (MASCARA) (Snellen et al., 2012; Talens et al., 2017). The
capability of bRing to monitor bright stars and find previously unknown
variables has been shown by Mellon et al. (2019). More information on the
observing strategy and design of bRing can be found in Stuik et al. (2017).
The reduction pipeline for the MASCARA and bRing instruments is described in
Talens et al. (2018). With a passband of 463 – 639 nm, bRing collected the
shortest wavelengths of all observatories considered in this work. We expected
to see the highest pulsational amplitudes in these data, as $\beta$ Pictoris
is a star of spectral type A6 (Zwintz et al., 2019) and has its energy maximum
in the blue optical wavelengths. Due to some evident outliers in the data, one
5$\sigma$ clip with respect to the median of the data set was applied. This
significantly weakened the one-day aliases in the spectral window. An
iterative sigma clipping procedure was not conducted due to noticeable changes
in the amplitudes of the pulsations in this case (see Hogg et al. (2010) for a
discussion of sigma clipping in order to remove outliers). The observations by
bRing were separated into two equally sized segments in order to gain more
time delay measurements while also maintaining a precision in frequency and
phase comparable to the ASTEP observations. Zwintz et al. (2019) found six
significant frequencies in the photometry collected by bRing. All of them are
also identified in the data collected by BRITE, ASTEP, and TESS.

###  3.3 ASTEP

The Antarctic Search for Transiting ExoPlanets, or ASTEP, is an automated
telescope with an aperture of 40 cm located at the Concordia station at Dome C
in Antarctica (Abe et al., 2013; Guillot et al., 2015; Mékarnia et al., 2017).
It uses a Sloan i’ filter (centered at 763 nm). We only used measurements with
a sun elevation lower than -18∘. Notably, data points where the centroid of
the star did not fall on the central pixel suffer from strong outliers. The
removal of these outliers and a 5$\sigma$ clip with respect to the median
weakened aliases significantly but without noticeable changes in the amplitude
of the strongest pulsational frequencies. Mékarnia et al. (2017) conducted a
frequency analysis of the $\beta$ Pictoris photometry collected by the ASTEP
observatory, and they are consistent with the ones seen in the TESS data.

###  3.4 TESS

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2015) was
launched in April 2018 in order to find transiting exoplanets around nearby
bright stars. The data of $\beta$ Pictoris (TIC 270577175, T = 3.696 mag) was
collected from 19 October 2018 to 1 February 2019 in sectors four through
seven and from 20 November 2020 to 8 February 2021 in the sectors 32 through
34. The data for first four sectors were obtained during TESS’ primary
mission, and data for the three other sectors were obtained approximately two
years later as part of the first extended mission. $\beta$ Pictoris is one of
the preselected targets for which short-cadence (2 minutes) data are provided.
Due to the high-cadence data, the high photometric precision of TESS, its high
duty cycle, and the long baseline, $\delta$ Scuti pulsations can be resolved
and identified with high precision. The photometric data of $\beta$ Pictoris
as observed by TESS was accessed and modified with the Python package
lightkurve (Lightkurve Collaboration et al., 2018), which retrieves the data
from the MAST archive.444https://archive.stsci.edu/tess/ For this analysis, we
used the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP;
Smith et al., 2012; Stumpe et al., 2012) light curves, which are produced from
the Science Processing Operations Center (SPOC) pipeline (Jenkins et al.,
2016; Jenkins, 2017; Jenkins et al., 2010). These PDCSAP light curves were
corrected for systematics by the SPOC pipeline. We also visually inspected the
target pixel files (TPF) in order to rule out various instrumental and
astrophysical effects, such as Solar System asteroids or comets crossing the
field of view. A comparison of the Lomb-Scargle periodogram (Lomb, 1976;
Scargle, 1982) of the raw simple aperture photometry and PDCSAP light curves
showed a significant change in the noise at low frequencies. This is due to
the systematic effects present in sector four. The lowest noise in the low-
frequency range can be found for the PDCSAP light curve with a completely
removed fourth sector. This light curve was then used for the main frequency
analysis.

The individual sectors were normalized by dividing each of the sectors by
their respective median flux, and the sectors were combined into one light
curve. Furthermore, every measurement with a non-zero quality flag (see Sect.
9 in the TESS Science Data Products Description
Document)555https://archive.stsci.edu/missions/tess/doc/EXP-TESS-ARC-ICD-
TM-0014.pdf was removed. Such anomalies as cosmic ray events or instrumental
issues were marked by these quality flags.

The frequency analysis was conducted using the Python package SMURFS (Müllner,
2020) and checked with the software package Period04 (Lenz & Breger, 2005).
Following Breger et al. (1993), all pulsation frequencies down to a signal-to-
noise ratio of four were extracted. The frequency range analyzed is between
zero and the Nyquist frequency of 360 cycles per day. Following the procedure
described in Zieba et al. (2019), 37 significant p-modes in frequencies
ranging from 34 to 76 $\text{d}^{-1}$ were identified. As we are only
interested in the strongest pulsational frequencies for this time delay
analysis, we did not try to further recover any of the lower amplitude modes.
A list of the extracted frequencies can be found in Table 5.

##  4 Theory and methodology

In this chapter, we discuss the theory behind time delays and the methods used
in order to finally arrive at the time delay plot. Importantly, this plot can
be used to search for companions around pulsating stars.

###  4.1 $\delta$ Scuti stars

$\delta$ Scuti stars can be found at the intersection region between the main
sequence and the instability strip on the Hertzsprung-Russel diagram. Thanks
to the nearly uninterrupted, high-precision photometry of Kepler’s primary
four-year mission, the general understanding of pulsating stars has been
revolutionized. $\delta$ Scuti stars have masses between approximately 1.5 and
2.5 $M_{\odot}$. They pulsate in radial and non-radial low-degree, low-order
pressure (p) modes that are primarily driven by an opacity mechanism (also
called a $\kappa$-mechanism) in their HeII zone with contributions from
turbulent pressure (Houdek, 2000; Antoci et al., 2014) and the edge-bump
mechanism (Murphy et al., 2020b). The oscillations have periods between 18
minutes and 8 hours respectively 80 and three cycles per day (Aerts et al.,
2010). Linear combinations of those oscillations can, however, create peaks at
lower frequencies (Breger & Montgomery, 2014). Besides main-sequence and more-
evolved stars, $\delta$ Scuti pulsations were observed in pre-main-sequence
stars, thus giving us the possibility to learn about early stellar evolution
(Zwintz et al., 2014; Murphy et al., 2021; Steindl et al., 2022).

###  4.2 The ephemeris equation

The search for time delays in certain astrophysical signals requires a
(quasi-)periodic process in space. A review on this and the related equations
can be found in Hermes (2018). There are different processes that are “clock-
like” under the assumption of a closed system, including the exceptionally
stable signals of pulsars, the eclipse time of binary stars, or certain
pulsating stars, as in our case. Deviations from periodic signals can be used
to analyze the spin-down of pulsars or to discover companions around pulsars
(Wolszczan & Frail, 1992; Wolszczan, 1994), eclipsing binaries (Barnes &
Moffett, 1975), or pulsating stars (Silvotti et al., 2007) (for a general
review of pulsating stars in binary systems see Murphy, 2018). To do so, one
creates O-C (observed minus calculated) diagrams (see e.g. Sterken, 2005) in
order to search for deviations from the predicted ephemeris in the
observations. O-C diagrams work the best if the star is pulsating in only a
single mode and if the maxima are narrow and well defined, as they are easy to
track in that case. However, these diagrams especially struggle with multi-
mode pulsators.

###  4.3 Frequency modulation and phase modulation: The state of the art

Building on the established methods of O-C diagrams, two new and complementary
techniques have emerged for finding companions around pulsating stars. The FM
method (Shibahashi & Kurtz, 2012) searches and analyzes the variations in the
frequency of a pulsating star induced by a companion. The periodic FM creates
multiples around every pulsation peak in the frequency spectrum. Their
frequencies, relative amplitudes, and phases can be used to get a full orbital
solution, as described in Shibahashi et al. (2015). The effectiveness of the
FM method was validated through a comparison with an eclipsing binary system
(Kurtz et al., 2015), and it is best suited for data sets with a baseline that
exceeds the orbital period of the companion.

The PM method is more sensitive to companions in wider orbits. It was
developed by Murphy et al. (2014), Murphy & Shibahashi (2015), and Murphy et
al. (2016b). Compton et al. (2016) showed that $\delta$ Scuti stars and white
dwarfs are best suited for this method. Its effectiveness was demonstrated by
Schmid et al. (2015) by showing the binary nature of KIC10080943 using the PM
method and attributing certain pulsations to the corresponding star in the
binary due to the antiphase modulation in time delays. Such a system with
observable time delays in both components is called a PB2, analogous to
spectroscopic terminology, where binary star systems are called SB2s if both
stars show observable RVs. Other proof of the functionality of the PM method
was shown by Derekas et al. (2019) by comparing the orbital parameters derived
from RV with those from PM.

An additional advantage of the PM method is its easier automation for many
stars. When applying this method to 2224 main-sequence A/F stars in the four-
year main Kepler data, Murphy et al. (2018) were able to find 317 PB1 systems,
where only one component is pulsating and showing time delays, and 24 PB2
systems, where two stars are pulsating. It is worth noting here that
determining orbital solutions using spectra and generating RV curves for the
same number of stars would be much more time intensive.

Other methods were developed by Koen (2014) and Balona (2014) to search for
binary systems by tracing the $\delta$ Scuti pulsations of stars. In contrast
to the FM and PM methods, these methods are not able to provide a full orbital
solution, which is usually gained by analyzing RV curves of spectroscopic
binaries.

###  4.4 Time delays

Time delays arise when a signal (in our case always an electromagnetic wave
with the propagation velocity defined by the speed of light) has to travel
different distances at different epochs. Following Smart (1977) and Balona
(2014), the distance $r$ between the pulsating star and the center of gravity
of its system can be described by

$r=\frac{a_{1}\left(1-e^{2}\right)}{1+e\cos f},$ (1)

where $a_{1}$ denotes the semi-major axis of the star, $e$ is its
eccentricity, and $f$ is the true anomaly. The distance to the star varies
relative to the Earth by

$z=r\sin(f+\varpi)\sin i,$ (2)

with $\varpi$ being the argument of periapsis, that is, the angle between the
nodal point and the periapsis,666The argument of periapsis is usually denoted
with $\omega$. This symbol, however, is used in asteroseismology to denote the
angular oscillation frequency. Also, one should not confuse $\varpi$ with the
longitude of periapsis, which is the sum of the longitude of the ascending
node $\Omega$ and the argument of periapsis. and $i$ as the inclination of the
system.

At this point, Equation 1 can be substituted into Equation 2. The time delay
$\tau=-z/c$ is then completely described by the following equation:

$\tau(t,\textbf{x})=-\frac{a_{1}\sin i}{c}(1-e^{2})\frac{\sin f\cos\varpi+\cos
f\sin\varpi}{1+e\cos f}.$ (3)

The set $\textbf{x}=(\Omega=2\pi/P,a_{1}\sin i/c,e,\varpi,t_{p})$ in Equation
3 includes all of the system-specific parameters needed to describe the time
delay for a given time $t$. The term $P$ is the orbital period of the system,
or equivalently $1/P=\nu_{\text{orb}}$ the orbital frequency, and thus
$\Omega$ is the angular orbital frequency. The projected semi-major axis of
the pulsating star is described by $a_{1}\sin i$. Dividing this quantity by
the speed of light $c$ gives us the size of the orbit for the pulsating star
in light seconds. The argument of periapsis is described by $\varpi$ and the
time of periapsis passage by $t_{p}$. (For a graphical visualization of the
orbital parameters, see Murphy & Shibahashi (2015).)

The two trigonometric functions of the true anomaly, $\sin f$ and $\cos f$,
can be expressed in terms of series expansions and Bessel functions:

$\cos f=-e+\frac{2\left(1-e^{2}\right)}{e}\sum_{n=1}^{\infty}J_{n}(ne)\cos
n\Omega\left(t-t_{\mathrm{p}}\right),$ (4)

$\sin f=2\sqrt{1-e^{2}}\sum_{n=1}^{\infty}J_{n}^{\prime}(ne)\sin
n\Omega\left(t-t_{\mathrm{p}}\right),$ (5)

with $J_{n}^{\prime}(x)=\text{d}J_{n}(x)/\text{d}x$ (the derivation of
Equation 4 and 5 can be found in Appendix A of Shibahashi et al., 2015). The
changing distances between us and the clock in space are fundamentally
connected with varying radial velocities, $v_{\text{rad }}$:

$v_{\text{rad}}=c\frac{\text{d}\tau}{\text{d}t}.$ (6)

By substituting Equation 3 into Equation 6, we obtain

$v_{\text{ rad }}=-\frac{\Omega a_{1}\sin
i}{\sqrt{1-e^{2}}}[\cos(f+\varpi)+e\cos\varpi].$ (7)

Given Equation 6 and the convention that a positive RV corresponds with a
receding object and a negative RV with an approaching one, we could deduce the
following: a negative time delay is due to an early arrival of the signal,
that is, the star is closer to us, and vice versa777This convention for the
time delays was established with Murphy & Shibahashi (2015), while Murphy et
al. (2014) used reversed signs;their plots are therefore mirrored around the
vertical axis. (see Table 4).

Table 4: Sign convention for the RV $v_{\text{rad}}$ and the time delays $\tau$. | positive sign (+) | negative sign (-)
---|---|---
$v_{\text{rad}}$ | moving away | approaching
$\tau$ | farther away / late arrival | closer / early arrival

One can see in Equations 3 and 7 that the time delay as well as the RV of a
system can be completely described by the orbital parameters. If we obtain
those parameters by one method, we can predict what we should observe with the
other one. Furthermore, if we generate the time delay plot from our
observations, we can apply a chi-squared minimization technique in order to
get the parameters in set x. This concept was introduced with Murphy &
Shibahashi (2015) and is a major improvement to Murphy et al. (2014), where
the time delay measurements were numerically differentiated in order to derive
the parameters from the obtained RV curve.

Finally, by using two of the derived orbital parameters, $a_{1}\sin i/c$ and
$P_{\text{orb}}$, we can calculate the mass function $f\left(m_{1},m_{2},\sin
i\right)$ for the binary system:

$f\left(m_{1},m_{2},\sin i\right):=\frac{\left(m_{2}\sin
i\right)^{3}}{\left(m_{1}+m_{2}\right)^{2}}=\frac{4\pi^{2}c^{3}}{G}v_{\text{
orb }}^{2}\left(\frac{a_{1}\sin i}{c}\right)^{3},$ (8)

where $m_{2}$ is the mass of the (usually non-pulsating) companion and $G$ is
the gravitational constant.

###  4.5 Phase modulation method: Methodology

Before we could create the time delay plot, we had to analyze the change in
the phase of the various pulsation modes with time. The basic equations for
that can be found in Murphy et al. (2014) and are summarized in the following.
We started by dividing the light curve into $n$ equally sized segments. Then,
we calculated the phase in every segment for each frequency. This left us with
a series of phases $\Phi_{j}$ for every segment $(1,2,\ldots,n)$ for a fixed
frequency $\nu_{j}$:

$\Phi_{j}=\left[\phi_{1j},\phi_{2j},\ldots,\phi_{ij},\ldots,\phi_{nj}\right].$
(9)

Numerically, the phase in a segment is derived by calculating the argument of
the Fourier Transformation in the respective segment:

$\Phi(t;\nu)=\tan^{-1}\left(\frac{\operatorname{Im}(\mathrm{F}(t;\nu,\delta
t))}{\operatorname{Real}(\mathrm{F}(t;\nu,\delta t))}\right),$ (10)

where $\mathrm{F}(t;\nu,\delta t)$ is the value of the Fourier Transformation
of the time series for frequency $\nu$ in segment $\delta t$.

As phases are frequency dependent, the resulting phase shifts have different
amplitudes for different frequencies. To get rid of this effect, we converted
them into time delays by first calculating the relative phase shifts:

$\Delta\phi_{ij}=\phi_{ij}-\overline{\phi}_{j},$ (11)

with $\overline{\phi}_{j}$ as the mean phase of frequency $\nu_{j}$:

$\overline{\phi}_{j}=\frac{1}{n}\sum_{i=1}^{n}\phi_{ij}.$ (12)

The time delay $\tau_{ij}$ for segment $i$ and frequency $\nu_{j}$ is thus
simply the relative phase shift divided by the angular pulsation frequency:

$\tau_{ij}=\frac{\Delta\phi_{ij}}{2\pi\nu_{j}}.$ (13)

A planet with an orbital frequency of $v_{\mathrm{orb}}$ in a circular orbit
will induce time delays that can be described by a sine function with phase
$\psi$:

$\tau(t)=A\sin\left(2\pi v_{\mathrm{orb}}t+\psi\right).$ (14)

The amplitude $A$ can be simply derived by inserting the center of mass
equation $m_{1}a_{1}=m_{2}a_{2}$ into $\tau=a_{1}\sin i/c$:

$\tau=\frac{a_{\text{P}}\sin i}{c}\frac{M_{\text{S}}}{M_{*}},$ (15)

with $a_{\text{S}}$ as the semi-major axis of the companion. The terms
$M_{\text{S}}$ and $M_{*}$ are the mass of the companion and the star,
respectively. Using Equation 15 and the mass of around 1.8 $M_{\odot}$ given
in Table 1, we could calculate the expected time delay for a given period.
This is visualized in Figure 2. The time delay is around 24 seconds for
$\beta$ Pictoris b and 6 seconds for $\beta$ Pictoris c. For comparison, the
smallest time delay detected in the main Kepler data is 7 seconds (Murphy et
al., 2016a).

Figure 2: Time delays for the $\beta$ Pictoris system. The colors indicate the
expected time delays for an edge-on planet in a circular orbit. The
uncertainties in the orbital period for the planets are smaller than the
marker size.

For more eccentric orbits, the pulsation time plot is described by a sum of
harmonics with amplitudes $A_{k}$ and phases $\phi_{k}$ corresponding to order
$k$:

$\tau(t)=\sum_{k=1}^{N}A_{k}\sin\left(2\pi
kv_{\mathrm{orb}}t+\psi_{k}\right).$ (16)

The height of the first harmonic relative to the one of the orbital frequency
is a measure of the eccentricity. The theory behind this is described in
Appendix A of Murphy et al. (2014). A visualization of that can be seen in
Figure 3. An increase in eccentricity also influences the amplitude of the
time delay. This is given by the following equation:

$\frac{a_{1}\sin
i}{c}=\frac{(\tau_{\text{max}}-\tau_{\text{min}})}{2}\left(1-e^{2}\cos^{2}\varpi\right)^{-1/2}.$
(17)

The maximum time delay is therefore reached in the case of $\varpi=\pm\pi/2$
or for the simple circular orbit case.

Figure 3: Simulation of a companion in a circular ($e=0$; left column) and
eccentric ($e=0.9$; right column) orbit as if observed by the Kepler Space
Telescope. The following parameters were used: $P=1$ year, $\varpi=0$,
$M_{\text{pulsating}}=1.8$ $M_{\odot}$, and $M_{\text{companion}}=0.1$
$M_{\odot}$. This led to a semi-amplitude of around 34 seconds in the circular
case (using Equation 15) and around 15 seconds in the eccentric case (using
Equation 17). Upper panel: Simulated time delay plot. Lower panel: Fourier
transformation of the time delays. One can clearly see the relative increase
of the first harmonic at two cycles/year for the eccentric case.

The larger the ratio between the orbital size $a_{1}\sin i/c$ and the
pulsation period $1/\nu$, the higher the sensitivity of the method (Murphy et
al., 2016b).

Due to the size of the segments, one has to make a trade-off between time or
frequency resolution. Using a shorter segment size has the advantage of a
finer sensitivity at periastron; however, the uncertainties are simultaneously
increased because of a poorer frequency resolution in the Fourier transform.

Under the assumption of Gaussian noise, increasing the cadence of an
observation by a factor of $N$ decreases the uncertainties in the measured
phases by a factor of $\sqrt{N}$ (Murphy, 2012). The phase errors also scale
inversely with amplitude, which means that the most valuable frequencies are
the ones with the highest amplitudes.

###  4.6 Intrinsic amplitude and phase variations

Amplitude modulations in $\delta$ Scuti stars have been observed in the past
and thoroughly analyzed in the four-year main Kepler data by Bowman et al.
(2016). Additionally, due to intrinsic reasons for those modulations (such as
the coupling of pulsational modes or pairs of close unresolved frequencies
leading to a beating effect), binarity can cause variability.

$\beta$ Pictoris is known to show amplitude variation in certain pulsational
frequencies, as reported by Zwintz et al. (2019) and Mékarnia et al. (2017).
However, PMs have not been observed yet (Zwintz et al., 2019).

###  4.7 Light curve reduction

Following Murphy et al. (2016a), unused frequencies were pre-whitened from our
light curves, as their presence adds unwanted variance to the data.
Furthermore, a high-pass filter was applied to the light curve to remove any
remaining instrumental signal and low-frequency oscillations, preserving all
content at frequencies above 5 d-1. The effect of a high-pass filter on low
frequencies can be seen in Figure 4.

Figure 4: Comparison of the amplitude spectra of the “raw” PDCSAP light curve
(upper panel) and the Gaussian high-pass filtered light curve (lower panel).
The power of the peaks below five d-1 are significantly weakened without
influencing the $\delta$ Scuti pulsations.

##  5 Results

To track the PM over all data sets, we started by determining which
frequencies have a signal-to-noise ratio greater than four in all
observations. This is the case for the four strongest frequencies in the TESS
data (the first four frequencies listed in Table 5). The stability of those
frequencies over the different observations is analyzed in Section 5.1. We
then looked at time delay curves created from simulated light curves. For
this, the “best-case scenario” of a four-year Kepler observation of $\beta$
Pictoris is studied in Section 5.2. The time delays of the real observations
and a comparison to a simulated re-creation can be found in Section 5.3.
Finally, we analyze the pulsational stability of the $\delta$ Scuti pulsations
of $\beta$ Pictoris using TESS data in Section 5.5.

###  5.1 Frequency stability between the different observations

As mentioned in Section 4.5, the PM method derives a time delay from the
observed PM at fixed frequencies. The precision with which pulsational
frequencies can be determined depends on the quality of the data (cadence,
timebase, precision, etc.). The photometry collected by the TESS mission has
the smallest uncertainties in frequency of all data sets (see Fig. 5). We
therefore used TESS as a “gold standard” for the frequencies used in the PM
method. The uncertainties in the frequencies were calculated following
Montgomery & Odonoghue (1999). However, as noted in their publication, these
errors are a lower limit of the true values. Keeping in mind that the actual
error bars are probably bigger, one can see that the frequencies are in
agreement with each other across the different data sets (Fig. 5).

Figure 5: Frequencies and their uncertainties over all data sets for the four
modes that are visible in all observations from the four different
observatories. The dashed red line marks the frequency determined by the TESS
mission, which has the smallest uncertainties. The uncertainties were
calculated following Montgomery & Odonoghue (1999).

###  5.2 Simulation based on Kepler data

Figure 6 shows the expected time delays for $\beta$ Pictoris caused by $\beta$
Pictoris b and $\beta$ Pictoris c using the full orbital solution given in
Table 2. The addition of the individual contributions on the phases gives the
total time delay curve as seen by the solid line in Figure 6.

Figure 6: Expected time delays for two planets in the $\beta$ Pictoris system.
The dashed line is for $\beta$ Pictoris b, the dashed-dotted line is for
$\beta$ Pictoris c, and the solid line shows both. The blue shaded region
marks the time span of Kepler’s four-year main mission; we note that Kepler
did not observe $\beta$ Pictoris.

In order to see how such properties as photometric precision, cadence, and
gaps in the observations influence the derived time delays, we simulated light
curves of $\beta$ Pictoris. We used the actual time stamps of the short-
cadence observations of Kepler (Borucki et al., 2010), which have a cadence of
around one minute. The simulations consist of a multi-sine of the frequencies
listed in Table 5. Using even more frequencies increases the computational
time without influencing the results of the simulations due to their low
amplitudes. The time stamps were then modulated by the expected time delay at
a given time using Equation 3 and assuming a two-planet configuration in this
system. Further, Gaussian noise on the order of 30 ppm was added to every data
point, which is comparable to the noise floor of TESS. Following the procedure
explained in Section 4.5, the light curve was separated into 20-day segments,
and the time delays were calculated from the phases in every segment with a
fixed frequency. Finally, we calculated weighted time delay values and their
corresponding uncertainties using the first three frequencies with the highest
amplitudes.

The first simulation (Fig. 7) used the Kepler short-cadence one-minute time
stamps. The measured time delays follow the prediction for a two-planet case.
Removing every second data point, which effectively reduces the number of
measurements by 50%, does not change the result (Fig. 8). One can, however,
observe a small increase in the uncertainties for the time delays. This is
expected, as a decrease of data points by a factor of two increases the
uncertainty by a factor of $\sqrt{2}$, assuming Gaussian noise (Murphy, 2012).
Finally, we induced gaps into the light curve, effectively simulating ground-
based observations by having data only for half of the day. Due to a worse
spectral window, there are many more peaks present in the amplitude spectrum.
These new peaks influence the phases for the observed frequencies, as they
cannot be resolved anymore. The scatter in the time delays increases
significantly (Fig. 9), and one cannot conclusively distinguish between a one-
planet solution (only $\beta$ Pictoris b) or a two-planet solution. A bigger
segment size mitigates this effect, as expected by the Rayleigh criterion.

This clearly illustrates that gaps influence the time delays the strongest, as
the uncertainties in phase only scale with the square root of the cadence
factor. One should also remove identified frequencies that are not used in the
time delay analysis to get rid of their aliases.

Figure 7: Derived time delay curve using 20-day segments by simulating Kepler
observations of $\beta$ Pictoris. Properties of the simulated light curve:
one-minute cadence, continuous observations, 20 ppm noise in flux. The
strongest frequencies (f1, f2, f3, …) are shown here with their uncertainties.
The weighted average of the measurements is shown in black. Figure 8: Derived
time delay curve using 20-day segments by simulating Kepler observations of
$\beta$ Pictoris. Properties of the simulated light curve: two-minute cadence,
continuous observations, 20 ppm noise in flux. Figure 9: Derived time delay
curve using 20-day segments by simulating Kepler observations of $\beta$
Pictoris. Properties of the simulated light curve: one-minute cadence, 0.5-day
gaps every day, 20 ppm noise in flux.

###  5.3 Time delay analysis of the photometry

Figure 10 shows the predicted time delays caused by the planets in the $\beta$
Pictoris system during the times when the observatories BRITE, ASTEP, bRing,
and TESS collected photometry for the star. The available observations have
been introduced in Section 3. The semi-amplitude of the predicted time delays
for $\beta$ Pictoris b and c is around 24 and 6 seconds, respectively.

As seen in Section 5.1, the TESS observations show the smallest uncertainties
in frequency and were therefore used as a “gold standard” in this analysis.
The frequency was thus fixed to the TESS values, as the PM method observes the
phase shifts at a constant frequency (see Section 4.5). The time delay
predictions (blue lines in Figure 10, 11 and 12) were also normalized to the
midpoint time of TESS. As a time delay is a relative measure and not an
absolute one, we set the time delay for TESS to zero. The evaluated time
delays shown in Figures 11 and 12 are therefore relative to the TESS values.

The code used to calculate the time delays was written for this analysis and
is heavily based on existing ones, namely
timedelay888https://github.com/danhey/timedelay and
maelstrom999https://github.com/danhey/maelstrom (Hey et al., 2020a). The
equations that were needed in order to evaluate the time delays are given in
Section 4.4 and 4.5.

The phases were calculated by subtracting the midpoint time of the full data
set. As discussed in Section 5.1, there are only four frequencies that are
significant in all observations. The phases for each data set were then
calculated based on a least-squares routine and their uncertainties from the
respective covariance matrices. Equation 13 gives the conversion between the
phase of a frequency and the respective time delays.

Figure 11 shows the derived time delays for the four different frequencies.
They are clearly not consistent with each other. As discussed in Section 4.6,
this rules out an extrinsic cause for the modulations (e.g., a companion), as
all frequencies would show a similar behavior (examples of this are shown in
Figure 3, 7, 8 and 9). A change in frequency was ruled out in Section 5.1.
Furthermore,Zwintz et al. (2019) showed no significant phase change for our
four frequencies in the BHr 2018 data set (therein, these four frequencies
have the designations F8, F11, F13, and F15).

Next, we attempted to reproduce the different data sets as faithfully as
possible and compare them with the time delay values shown in Figure 11. For
that, we first determined the frequency, amplitude, and phase of the four
pulsation modes visible in all observations. We calculated the residual noise
for the pre-whitened data sets, which was then used to estimate the
uncertainties following Montgomery & Odonoghue (1999). As before, we fixed the
frequency to the TESS value. As we do not know the exact “true” frequency of
the pulsations with infinite precision, we introduced an offset between the
true pulsational frequencies and the TESS data set in the simulations. The
uncertainty in frequency for these four strongest frequencies is on the order
of $10^{-5}$ d-1 (see Table 5). This offset explains the linear trend for
every frequency that is visible in Figure 11. The time delays of the simulated
data set are shown in Figure 12.

This linear trend was discussed on simulated data in Section 3.2 of Murphy et
al. (2016b) and “almost certainly” explains the observed trend in the WASP
data of Murphy et al. (2013). A way to correct for it is to evaluate the slope
between two maxima or minima of the sinusoidal variations. This is not a
possibility in our case, as we would have had to further segment the data sets
to identify the position of the maxima or minima, leading to even higher
scatter in the phase. The lower panel of Figure 11 also shows that the
uncertainties in the derived time delays for the data sets other than TESS are
too big to differentiate between a one-planet or two-planet scenario, even
without this linear trend. We therefore discuss the possibility of a second
$\beta$ Pictoris observation by TESS in the extended mission in the next
section.

Figure 10: Time delay predictions for $\beta$ Pictoris b (dashed line), c
(dashed-dotted line), and both planets (solid line). Times when the star was
observed are marked with lines. Figure 11: Time delay plot calculated from the
phases of four different frequencies for all available observations by BRITE,
bRing, ASTEP and TESS. Each color represents a frequency (f1, f2, f3, and f4)
listed in Table 5. The blue lines indicate time delay predictions for $\beta$
Pictoris b (dashed line), c (dashed, dotted line), and both planets (solid
line). The lower panel is a zoom-in of the upper panel. The uncertainties in
the time delays were derived from the covariance matrices given by the least-
squares procedure, which was used in order to calculate the phases of the
respective frequencies. The ticks at the top of the plot denote the various
observatories: B15, B16, and B17 for the BRITE observations in the years 2015,
2016, and 2017, respectively; A17 and A18 for the ASTEP observations in 2017
and 2018; and R1 and R2 for the first and second part of the bRing data.
Figure 12: Time delay plot for the simulated data set showing a high
similarity to the time delays of the real data set presented in Figure 12. The
colored points represent the simulated time delays for the four strongest
pulsational frequencies. The blue lines indicate time delay predictions for
$\beta$ Pictoris b (dashed line), c (dashed-dotted line), and both planets
(solid line). The lower panel is a zoom-in of the upper panel. A description
of the ticks at the top of the plot can be found in the caption of Figure 11.

###  5.4 Detection limits for $\beta$ Pictoris

Figure 13:  Detection limits for the $\beta$ Pictoris system based on the
calculations presented in Hey et al. (2020b) assuming 1.8 $M_{\odot}$ for the
star and an inclination of $90^{\circ}$ for both planets. The detectable
parameter space for companions depending on their orbital period
$P_{\textrm{orb}}$ and mass is shaded red, assuming $\beta$ Pictoris is a non-
ideal $\delta$ Scuti pulsator (see Section 5.5 for a discussion on the
pulsational modes of $\beta$ Pictoris as seen by TESS). $\beta$ Pictoris c is
generally not detectable due to the intrinsic noise of the pulsations, and the
time delay caused by $\beta$ Pictoris b has a period that is too long in
comparison to the baseline of our observations. The masses corresponding to
brown dwarfs (13 M${}_{\textrm{J}}$ $\lesssim$ $M$ $\lesssim$ 80
M${}_{\textrm{J}}$) are shaded in yellow, and the stellar regime ($M$ $>$ 80
M${}_{\textrm{J}}$) is in magenta, with M${}_{\textrm{J}}$ being a Jupiter
mass. The figure is adapted from Hey et al. (2020b).

Here, we determine the detection limits for companions in the $\beta$ Pictoris
system using the PM method. In a prior study by Hey et al. (2020b), the
authors tried to estimate the detection limits of companions around $\delta$
Scuti stars. Their approach involved simulating time-series observations of
these pulsating stars and adding white noise to the data. This allowed them to
find a relationship between the S/N of the stellar pulsations and the observed
scatter in the resulting time delay ($a\,\sin i/c$). The established
relationship could be directly converted into a parameter space defining
detectable companion masses and their orbital periods. The authors found that
the detectability of companions strongly depends on the S/N of the stellar
pulsations. Given the comparably low observed S/N of the $\delta$ Scuti
pulsations observed in $\beta$ Pictoris and the instability of the pulsational
modes, we opted for a more conservative detection limit than what
(84${}^{\textrm{th}}$ percentile in Hey et al., 2020b) presented in their
prior work.

Our findings are presented in Figure 13. We determined that the intrinsic
variability of $\beta$ Pictoris is too high to detect planet c. The other
companion, $\beta$ Pictoris b, is primarily not detectable at the moment due
to the short baseline of observations compared to its long orbital period of
approximately 24 years. It is worth noting that these calculations assume that
the pulsational modes stay stable during the time of observations. However, we
show in Section 5.5 by using the TESS data that this assumption is generally
not met, with some modes appearing and disappearing during the observations.
This further complicates any efforts to detect companions around the star.

###  5.5 Analysis of pulsational stability using TESS data

$\beta$ Pictoris was observed in seven individual sectors between October 2018
and February 2021 (see Table 3). We performed a frequency analysis of these
TESS sectors using maelstrom. We find that $\beta$ Pictoris is seemingly
undergoing significant frequency and amplitude modulation, which buries any
signal induced by planetary companions. The periodograms of the stellar
pulsations clearly show significant amplitude modulation by the rotational
signal (see Fig. 14). We also find that the star is showing modes that are
appearing and disappearing on short timescales. Figure 15 shows a mode that
seems to be just appearing during the second half of the TESS observations. In
summary, we find that $\beta$ Pictoris’ modes are not stable enough to probe
to the necessary phase precision for the planetary companions.

Figure 14: Two-dimensional periodogram showing the frequency region between 45
and 55 d-1. One can see that most of the modes are significantly amplitude
modulated by the rotational signal. The mode at 54 d-1 goes much faster, which
is probably due to beating with nearby modes. The mode at 50.5 d-1 itself
undergoes incoherent FM. Figure 15: Periodogram showing a pulsational mode
around 47.3 d-1 that is just appearing at the end of the observations.

###  5.6 Comparison to KIC 7917485

Here we compare the $\beta$ Pictoris system to another A star with a planet
detected through pulsation timing, KIC 7917485 (Murphy et al., 2016a), and
evaluate the differences between the host stars that have affected the
detectability of their corresponding planets. Compared to other Kepler
$\delta$ Scuti stars, including those with binary star companions, (e.g.,
Murphy et al., 2018, 2020a), and Murphy et al. (2016a) found that KIC 7917485
had an exceptionally low time delay noise (sixth lowest of the 2040 $\delta$
Scuti stars in the Kepler primary mission sample). Similar planetary-mass
companions might exist around the second and ninth lowest noise stars in the
Kepler sample (KIC 9700322 and KIC 8453431). Unlike KIC 7917485, however, the
observations of these systems and their potential companions do not cover a
full orbital period, thus precluding determination of their orbital
parameters.
The additional time-delay jitter seen in $\beta$ Pictoris might arise from
mode interaction or from other effects causing the observed changes in mode
amplitude. It is this jitter that hinders the ability to detect planetary-mass
objects around pulsating stars. This is comparable to the intrinsic RV jitter
of a star, which impacts the ability to detect planets using RV measurements.
The intrinsic variability is evident in the $\delta$ Scuti pulsations of
$\beta$ Pictoris and generally renders white dwarfs and subdwarfs as less
ideal targets for the analysis of their time delays using the PM method
(Murphy, 2018). This limitation persists despite the high oscillation
frequencies exhibited by these pulsating stars that would otherwise make them
promising targets for detecting time delays caused by a companion (Compton et
al., 2016).

##  6 Conclusions

In this work, we have analyzed the time delays derived from the phases of the
$\delta$ Scuti pulsations of $\beta$ Pictoris. The photometric data of the
star were collected over a time period of approximately four years by four
different observatories: the BRITE-Constellation, bRing, ASTEP, and TESS. In
contrast to previous studies, we did not segment the observations into smaller
sets (e.g. ten-day bins). This would have caused high uncertainties in the
phases and therefore also in the time delays. Nevertheless, we could not see
the influence of $\beta$ Pictoris b or c in the data due to this time delay
scatter. The uncertainty in the pulsational frequencies leads to a linear
trend in the time delays and has also been seen in a previous study by Murphy
et al. (2013) and in simulations by Murphy et al. (2016b). We performed a
frequency analysis using the open-source tool maelstrom. We find that $\beta$
Pictoris does not have the needed stability to detect planetary companions
using the time delay method. The stellar pulsations clearly show strong
amplitude modulation caused by the rotational signal and identify modes that
seemingly appeared during our observations.

Previous studies have used the PM method on Kepler data, finding many binary
star systems (Murphy et al., 2016b, ; and references therein) and a planet
(Murphy et al., 2016a). This work is the first to use the PM method with so
many different data sets that have a precision significantly lower than the
Kepler mission. However, the PM method remains valuable technique, as it is
able to find planets and stars in a parameter space that is poorly covered by
other methods, such as the RV method (see e.g. Murphy, 2018).

###### Acknowledgements.

This work includes data collected by the TESS mission, which are publicly
available from the Mikulski Archive for Space Telescopes (MAST). Funding for
the TESS mission is provided by the NASA Explorer Program. This work has made
use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. We made use of the software
package Period04 (Lenz & Breger, 2005), the Python programming language
(Rossum, 1995), and the open-source Python packages numpy (van der Walt et
al., 2011), matplotlib (Hunter, 2007), astropy (Astropy Collaboration et al.,
2013), lightkurve (Lightkurve Collaboration et al., 2018),
timedelay101010https://github.com/danhey/timedelay,
maelstrom111111https://github.com/danhey/maelstrom and SMURFS (Müllner, 2020).
This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France. This research has made use of the NASA Exoplanet Archive,
which is operated by the California Institute of Technology, under contract
with the National Aeronautics and Space Administration under the Exoplanet
Exploration Program. Part of this research was carried out at the Jet
Propulsion Laboratory, California Institute of Technology, under a contract
with the National Aeronautics and Space Administration (80NM0018D0004). This
research has made use of data collected by the BRITE-Constellation satellite
mission, designed, built, launched, operated and supported by the Austrian
Research Promotion Agency (FFG), the University of Vienna, the Technical
University of Graz, the University of Innsbruck, the Canadian Space Agency
(CSA), the University of Toronto Institute for Aerospace Studies (UTIAS), the
Foundation for Polish Science & Technology (FNiTP MNiSW), and National Science
Centre (NCN). The bRing observatory at Siding Springs, Australia was supported
by a University of Rochester University Research Award. The field activities
at Dome C for ASTEP benefit from the support of the French and Italian polar
agencies IPEV and PNRA in the framework of the Concordia station programme.
The PicSat team thanks funding from the European Research Council (ERC) under
the Horizon 2020 research and innovation programme (Grant agreement No.
639248, LITHIUM). SJM was supported by the Australian Research Council (ARC)
through Future Fellowship FT210100485. Finally, we thank the anonymous referee
for a detailed report, which helped us to improve the quality of this
manuscript.

## References

  * Abe et al. (2013) Abe, L., Gonçalves, I., Agabi, A., et al. 2013, A&A, 553, A49
  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Springer Dordrecht)
  * Antoci et al. (2014) Antoci, V., Cunha, M., Houdek, G., et al. 2014, ApJ, 796, 118
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Augereau et al. (2001) Augereau, J. C., Nelson, R. P., Lagrange, A. M., Papaloizou, J. C. B., & Mouillet, D. 2001, A&A, 370, 447
  * Aumann et al. (1984) Aumann, H. H., Gillett, F. C., Beichman, C. A., et al. 1984, ApJ, 278, L23
  * Balona (2014) Balona, L. A. 2014, MNRAS, 443, 1946
  * Barnes & Moffett (1975) Barnes, T. G., I. & Moffett, T. J. 1975, AJ, 80, 48
  * Bedding et al. (2020) Bedding, T. R., Murphy, S. J., Hey, D. R., et al. 2020, Nature, 581, 147
  * Beust & Morbidelli (2000) Beust, H. & Morbidelli, A. 2000, Icarus, 143, 170
  * Blunt et al. (2020) Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Bowman et al. (2016) Bowman, D. M., Kurtz, D. W., Breger, M., Murphy, S. J., & Holdsworth, D. L. 2016, MNRAS, 460, 1970
  * Breger & Montgomery (2014) Breger, M. & Montgomery, M. H. 2014, ApJ, 783, 89
  * Breger et al. (1993) Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482
  * Burrows et al. (1995) Burrows, C. J., Krist, J. E., Stapelfeldt, K. R., & WFPC2 Investigation Definition Team. 1995, in American Astronomical Society Meeting Abstracts, Vol. 187, American Astronomical Society Meeting Abstracts, 32.05
  * Compton et al. (2016) Compton, D. L., Bedding, T. R., Murphy, S. J., & Stello, D. 2016, MNRAS, 461, 1943
  * Cousins (1971) Cousins, A. W. J. 1971, Royal Observatory Annals, 7
  * Derekas et al. (2019) Derekas, A., Murphy, S. J., Dálya, G., et al. 2019, MNRAS, 486, 2129
  * Ferlet et al. (1987) Ferlet, R., Hobbs, L. M., & Madjar, A. V. 1987, A&A, 185, 267
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Galland et al. (2006) Galland, F., Lagrange, A. M., Udry, S., et al. 2006, A&A, 447, 355
  * Golimowski et al. (2006) Golimowski, D. A., Ardila, D. R., Krist, J. E., et al. 2006, AJ, 131, 3109
  * Gray et al. (2006) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161
  * Guillot et al. (2015) Guillot, T., Abe, L., Agabi, A., et al. 2015, Astronomische Nachrichten, 336, 638
  * Heap et al. (2000) Heap, S. R., Lindler, D. J., Lanz, T. M., et al. 2000, ApJ, 539, 435
  * Hermes (2018) Hermes, J. J. 2018, Timing by Stellar Pulsations as an Exoplanet Discovery Method (Springer Cham), 6
  * Hey et al. (2020a) Hey, D., Murphy, S., Foreman-Mackey, D., et al. 2020a, The Journal of Open Source Software, 5, 2125
  * Hey et al. (2020b) Hey, D. R., Murphy, S. J., Foreman-Mackey, D., et al. 2020b, AJ, 159, 202
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686
  * Houdek (2000) Houdek, G. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 210, Delta Scuti and Related Stars, ed. M. Breger & M. Montgomery, 454
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Jenkins (2017) Jenkins, J. M. 2017, Kepler Data Processing Handbook: Overview of the Science Operations Center, Tech. rep., NASA Ames Research Center
  * Jenkins et al. (2010) Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010, ApJ, 713, L87
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, 99133E
  * Kalas et al. (2019) Kalas, P., Zwintz, K., Kenworthy, M., et al. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 218.03
  * Kanodia & Wright (2018) Kanodia, S. & Wright, J. T. 2018, Barycorrpy: Barycentric velocity calculation and leap second management
  * Kenworthy (2017) Kenworthy, M. 2017, Nature Astronomy, 1, 0099
  * Kenworthy et al. (2021) Kenworthy, M. A., Mellon, S. N., Bailey, J. I., et al. 2021, A&A, 648, A15
  * Koen (2003) Koen, C. 2003, MNRAS, 341, 1385
  * Koen (2014) Koen, C. 2014, MNRAS, 444, 1486
  * Kurtz et al. (2015) Kurtz, D. W., Hambleton, K. M., Shibahashi, H., Murphy, S. J., & Prša, A. 2015, MNRAS, 446, 1223
  * Lacour et al. (2021) Lacour, S., Wang, J. J., Rodet, L., et al. 2021, A&A, 654, L2
  * Lagrange et al. (2000) Lagrange, A. M., Backman, D. E., & Artymowicz, P. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell, 639
  * Lagrange et al. (2019a) Lagrange, A. M., Boccaletti, A., Langlois, M., et al. 2019a, A&A, 621, L8
  * Lagrange et al. (2010) Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57
  * Lagrange et al. (2012) Lagrange, A. M., De Bondt, K., Meunier, N., et al. 2012, A&A, 542, A18
  * Lagrange et al. (2009a) Lagrange, A. M., Desort, M., Galland, F., Udry, S., & Mayor, M. 2009a, A&A, 495, 335
  * Lagrange et al. (2009b) Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009b, A&A, 493, L21
  * Lagrange et al. (2019b) Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019b, Nature Astronomy, 421
  * Lagrange et al. (2020) Lagrange, A.-M., Rubini, P., Nowak, M., et al. 2020, A&A
  * Lecavelier des Etangs et al. (2022) Lecavelier des Etangs, A., Cros, L., Hébrard, G., et al. 2022, Scientific Reports, 12, 5855
  * Lecavelier Des Etangs et al. (1995) Lecavelier Des Etangs, A., Deleuil, M., Vidal-Madjar, A., et al. 1995, A&A, 299, 557
  * Lecavelier des Etangs & Vidal-Madjar (2016) Lecavelier des Etangs, A. & Vidal-Madjar, A. 2016, A&A, 588, A60
  * Lenz & Breger (2005) Lenz, P. & Breger, M. 2005, Communications in Asteroseismology, 146, 53
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M. a., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2
  * Lomb (1976) Lomb, N. R. 1976, Astrophysics and Space Science, 39, 447
  * Mamajek & Bell (2014) Mamajek, E. E. & Bell, C. P. M. 2014, MNRAS, 445, 2169
  * Mékarnia et al. (2017) Mékarnia, D., Chapellier, E., Guillot, T., et al. 2017, A&A, 608, L6
  * Mellon et al. (2019) Mellon, S. N., Mamajek, E. E., Stuik, R., et al. 2019, ApJS, 244, 15
  * Menegaldo et al. (2022) Menegaldo, C. G., de O. Fialho, F., Janot-Pacheco, E., Pait, F. M., & Lapeyrère, V. 2022, PASP, 134, 034501
  * Montgomery & Odonoghue (1999) Montgomery, M. H. & Odonoghue, D. 1999, Delta Scuti Star Newsletter, 13, 28
  * Mouillet et al. (1997) Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896
  * Müllner (2020) Müllner, M. 2020, MarcoMuellner/SMURFS, Zenodo
  * Murphy (2012) Murphy, S. J. 2012, MNRAS, 422, 665
  * Murphy (2018) Murphy, S. J. 2018, arXiv e-prints, arXiv:1811.12659
  * Murphy et al. (2020a) Murphy, S. J., Barbara, N. H., Hey, D., Bedding, T. R., & Fulcher, B. D. 2020a, MNRAS, 493, 5382
  * Murphy et al. (2016a) Murphy, S. J., Bedding, T. R., & Shibahashi, H. 2016a, ApJ, 827, L17
  * Murphy et al. (2014) Murphy, S. J., Bedding, T. R., Shibahashi, H., Kurtz, D. W., & Kjeldsen, H. 2014, MNRAS, 441, 2515
  * Murphy et al. (2021) Murphy, S. J., Joyce, M., Bedding, T. R., White, T. R., & Kama, M. 2021, MNRAS, 502, 1633
  * Murphy et al. (2018) Murphy, S. J., Moe, M., Kurtz, D. W., et al. 2018, MNRAS, 474, 4322
  * Murphy et al. (2013) Murphy, S. J., Pigulski, A., Kurtz, D. W., et al. 2013, MNRAS, 432, 2284
  * Murphy et al. (2020b) Murphy, S. J., Saio, H., Takada-Hidai, M., et al. 2020b, MNRAS, 498, 4272
  * Murphy & Shibahashi (2015) Murphy, S. J. & Shibahashi, H. 2015, MNRAS, 450, 4475
  * Murphy et al. (2016b) Murphy, S. J., Shibahashi, H., & Bedding, T. R. 2016b, MNRAS, 461, 4215
  * Nesvold & Kuchner (2015) Nesvold, E. R. & Kuchner, M. J. 2015, ApJ, 798, 83
  * Nowak et al. (2018) Nowak, M., Lacour, S., Crouzier, A., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10698, Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, 1069821
  * Nowak et al. (2020) Nowak, M., Lacour, S., Lagrange, A.-M., et al. 2020, A&A
  * Pavlenko et al. (2022) Pavlenko, Y., Kulyk, I., Shubina, O., et al. 2022, A&A, 660, A49
  * Popowicz et al. (2017) Popowicz, A., Pigulski, A., Bernacki, K., et al. 2017, A&A, 605, A26
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rossum (1995) Rossum, G. 1995, Python Reference Manual, Tech. rep., Centrum voor Wiskunde en Informatica (CWI), Amsterdam, The Netherlands, The Netherlands
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Schmid et al. (2015) Schmid, V. S., Tkachenko, A., Aerts, C., et al. 2015, A&A, 584, A35
  * Shibahashi & Kurtz (2012) Shibahashi, H. & Kurtz, D. W. 2012, MNRAS, 422, 738
  * Shibahashi et al. (2015) Shibahashi, H., Kurtz, D. W., & Murphy, S. J. 2015, MNRAS, 450, 3999
  * Silvotti et al. (2007) Silvotti, R., Schuh, S., Janulis, R., et al. 2007, Nature, 449, 189
  * Smart (1977) Smart, W. M. 1977, Textbook on Spherical Astronomy, 6th edn. (Cambridge University Press)
  * Smith & Terrile (1984) Smith, B. A. & Terrile, R. J. 1984, Science, 226, 1421
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, Publications of the Astronomical Society of the Pacific, 124, 1000
  * Snellen et al. (2012) Snellen, I. A. G., Stuik, R., Navarro, R., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8444, Proc. SPIE, 84440I
  * Stassun et al. (2019) Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138
  * Steindl et al. (2022) Steindl, T., Zwintz, K., & Müllner, M. 2022, A&A, 664, A32
  * Sterken (2005) Sterken, C. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 335, The Light-Time Effect in Astrophysics: Causes and cures of the O-C diagram, ed. C. Sterken, 3
  * Strøm et al. (2020) Strøm, P. A., Bodewits, D., Knight, M. M., et al. 2020, PASP, 132, 101001
  * Stuik et al. (2017) Stuik, R., Bailey, J. I., Dorval, P., et al. 2017, A&A, 607, A45
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, Publications of the Astronomical Society of the Pacific, 124, 985
  * Talens et al. (2017) Talens, G. J. J., Albrecht, S., Spronck, J. F. P., et al. 2017, A&A, 606, A73
  * Talens et al. (2018) Talens, G. J. J., Deul, E. R., Stuik, R., et al. 2018, A&A, 619, A154
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Wang et al. (2016) Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97
  * Weiss et al. (2014) Weiss, W. W., Rucinski, S. M., Moffat, A. F. J., et al. 2014, Publications of the Astronomical Society of the Pacific, 126, 573
  * Wolszczan (1994) Wolszczan, A. 1994, Science, 264, 538
  * Wolszczan & Frail (1992) Wolszczan, A. & Frail, D. A. 1992, Nature, 355, 145
  * Wright & Eastman (2014) Wright, J. T. & Eastman, J. D. 2014, PASP, 126, 838
  * Zieba et al. (2019) Zieba, S., Zwintz, K., Kenworthy, M. A., & Kennedy, G. M. 2019, A&A, 625, L13
  * Zwintz et al. (2014) Zwintz, K., Fossati, L., Ryabchikova, T., et al. 2014, Science, 345, 550
  * Zwintz et al. (2019) Zwintz, K., Reese, D. R., Neiner, C., et al. 2019, A&A, 627, A28

##  Appendix A TESS frequency analysis

###  A.1 Frequency list

Table 5: Pulsational frequencies, amplitudes in instrumental millimagnitudes, and normalized flux in parts per million, phases, and signal-to-noise ratio sorted by the pre-whitening sequence. # | Freq. ($d^{-1}$) | Ampl. (mmag) | Ampl. (ppm) | Phase | S/N
---|---|---|---|---|---
1 | 47.43895(6) | 1.029(9) | 948(9) | 0.9071(14) | 20.6
2 | 53.69166(7) | 0.948(9) | 873(9) | 0.2782(16) | 19.4
3 | 50.49168(7) | 0.926(9) | 852(9) | 0.5567(16) | 23.4
4 | 54.23716(12) | 0.553(9) | 509(9) | 0.982(3) | 22.8
5 | 39.06315(15) | 0.442(9) | 407(9) | 0.699(3) | 22.5
6 | 46.54259(16) | 0.415(9) | 382(9) | 0.391(4) | 18.7
7 | 48.9192(3) | 0.230(9) | 212(9) | 0.950(6) | 17.3
8 | 43.5283(3) | 0.214(9) | 197(9) | 0.029(7) | 19.9
9 | 47.2853(4) | 0.182(9) | 168(9) | 0.186(8) | 16.2
10 | 57.4525(4) | 0.164(9) | 151(9) | 0.577(9) | 18.0
11 | 34.7605(5) | 0.143(9) | 131(9) | 0.754(10) | 23.7
12 | 38.1297(5) | 0.131(9) | 121(9) | 0.980(11) | 20.2
13 | 45.2698(5) | 0.120(9) | 110(9) | 0.411(12) | 12.8
14 | 51.4969(6) | 0.118(9) | 109(9) | 0.397(13) | 14.7
15 | 47.2686(7) | 0.093(9) | 85(9) | 0.419(16) | 12.3
16 | 50.8310(8) | 0.086(9) | 79(9) | 0.630(17) | 12.9
17 | 49.7131(8) | 0.085(9) | 78(9) | 0.290(17) | 11.3
18 | 53.8545(8) | 0.085(9) | 78(9) | 0.566(17) | 9.5
19 | 44.6833(8) | 0.084(9) | 77(9) | 0.297(18) | 11.3
20 | 65.1356(8) | 0.083(9) | 76(9) | 0.350(18) | 17.9
21 | 43.8292(8) | 0.082(9) | 76(9) | 0.555(18) | 16.3
22 | 49.5595(8) | 0.079(9) | 73(9) | 0.980(19) | 13.4
23 | 42.0365(9) | 0.077(9) | 71(9) | 0.327(19) | 11.1
24 | 54.2269(9) | 0.073(9) | 67(9) | 0.51(2) | 10.0
25 | 41.6498(9) | 0.071(9) | 65(9) | 0.59(2) | 13.0
26 | 48.1381(10) | 0.064(9) | 59(9) | 0.23(2) | 11.2
27 | 45.8998(10) | 0.064(9) | 59(9) | 0.73(2) | 12.4
28 | 50.2689(12) | 0.054(9) | 50(9) | 0.90(3) | 12.1
29 | 75.6780(13) | 0.052(9) | 48(9) | 0.68(3) | 12.2
30 | 58.3469(13) | 0.050(9) | 46(9) | 0.61(3) | 11.3
31 | 45.4375(14) | 0.047(9) | 44(9) | 0.00(3) | 12.3
32 | 54.4625(14) | 0.047(9) | 43(9) | 0.16(3) | 8.1
33 | 53.6827(15) | 0.042(9) | 39(9) | 0.16(3) | 7.3
34 | 53.5521(16) | 0.040(9) | 37(9) | 0.48(4) | 7.7
35 | 42.1735(16) | 0.040(9) | 37(9) | 0.70(4) | 9.5
36 | 58.2515(17) | 0.039(9) | 36(9) | 0.81(4) | 10.5
37 | 42.3963(17) | 0.039(9) | 36(9) | 0.68(4) | 11.5

###  A.2 Gaussian high-pass filter

By applying a Gaussian high-pass filter on the TESS light curve (see Fig. 16),
the long-term variations such as systematics and the exocomets are
significantly weakened. At the same time however, the higher frequencies
related to the $\delta$ Scuti pulsations are preserved.

Figure 16: Comparison of the PDCSAP light curve (black in the background) and
the Gaussian high-pass filter (red dots) of it. The Gaussian high-pass filter
clearly shows less long-term variations in the light curve.