# Classifying Be Star Variability With TESS. I. The Southern Ecliptic

Jonathan Labadie-Bartz, Alex C. Carciofi, Tajan Henrique de Amorim, Amanda Rubio, André Luiz Figueiredo, Pedro Ticiani dos Santos, Keegan Thomson-Paressant

► To cite this version:

Jonathan Labadie-Bartz, Alex C. Carciofi, Tajan Henrique de Amorim, Amanda Rubio, André Luiz Figueiredo, et al.. Classifying Be Star Variability With TESS. I. The Southern Ecliptic. The Astronomical Journal, 2022, 163, 10.3847/1538-3881/ac5abd . insu-03713324

HAL Id: insu-03713324

https://insu.hal.science/insu-03713324

Submitted on 4 Jul 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Classifying Be Star Variability With TESS. I. The Southern Ecliptic

Jonathan Labadie-Bartz<fnref n="1" />, Alex C. Carciofi<fnref n="1" />, Tajan Henrique de Amorim<fnref n="1" />, Amanda Rubio<fnref n="1" />, André Luiz Figueiredo<fnref n="1" />, Pedro Ticiani dos Santos<fnref n="1" />, and Keegan Thomson-Paressant<fnref n="2" />

## Abstract

TESS photometry is analyzed for 430 classical Be stars observed in the first year of the mission. The often complex and diverse variability of each object in this sample is classified to obtain an understanding of the behavior of this class as a population. Ninety-seven percent of the systems are variable above the noise level, with timescales spanning nearly the entire range of what is accessible with TESS, from tens of minutes to tens of days. The variability seen with TESS is summarized as follows. Nearly every system contains multiple periodic signals in the frequency regime between about 0.5 and 4 day⁻<fnref n="1" />. One or more groups of closely spaced frequencies is the most common feature, present in 87% of the sample. Among the Be stars with brightening events that are characteristic of mass ejection episodes (18% of the full sample, or 31% of early-type stars), all have at least one frequency group, and the majority of these (83%) show a concurrent temporary amplitude enhancement in one or more frequency groups. About one-third (34%) of the sample is dominated by longer-term trends (timescales >2 day). Low-frequency stochastic signals are prominent in about 25% of the sample, with varying degrees of intensity. Higher-frequency signals (6 < *f* < 15 day⁻<fnref n="1" />) are sometimes seen (in 15% of the sample) and in most cases likely reflect p-mode pulsation. In rare cases (~3%), even higher frequencies beyond the traditional p-mode regime (*f* > 15 day⁻<fnref n="1" />) are observed.

*Unified Astronomy Thesaurus concepts: Be stars (142); Pulsating variable stars (1307); Stellar mass loss (1613); Photometry (1234); Light curves (918); Early-type emission stars (428); Circumstellar disks (235); Multi-periodic variable stars (1079)*

*Supporting material: machine-readable table*

## 1. Introduction

Classical Be stars have been studied for over 150 yr, yet key aspects of their nature remain veiled. Since their discovery in 1866 (Secchi 1866), it has been established that classical Be stars (here after simply Be stars) are nonsupergiant B-type stars with rapid, near-critical rotation, and are in general nonradial pulsators, which nonradiatively eject mass to form a viscous, near-Keplerian circumstellar “decretion” disk from which characteristic spectral emission features arise (Rivinius et al. 2013, and references therein). While significant progress has been made in the past several decades, the following questions regarding Be stars remain outstanding in general: How do they acquire such fast rotation? Why do they pulsate the way they do? How does rapid rotation influence internal processes such as angular momentum transport and chemical mixing? How are matter and angular momentum ejected at sufficiently high amounts? How does this ejected matter organize itself around the star on a geometrically thin and mostly Keplerian disk?

Mass, rotation, binarity, metallicity, and magnetic fields are the primary factors that dictate the life of a star. Be stars span the entire spectral type range from late O to early A, do not host large-scale magnetic fields down to a detection limit of about 50–100 Gauss (Wade et al. 2016, in a sample of 85 Be stars), are more common in lower-metallicity environments

(Maeder et al. 1999; Wisniewski & Bjorkman 2006; Peters et al. 2020), and are the most rapidly rotating nondegenerate class of objects known, on average rotating at or above 80% of their critical rotation rate (Frémat et al. 2005; Rivinius et al. 2013). The binary fraction and binary parameters of Be stars as a population are somewhat uncertain, but there is evidence that the binary parameters are similar to those of B-type stars in general (Oudmaijer & Parr 2010), as well as suggestions that a very high fraction of Be stars exist in binaries (Klement et al. 2019), often with evolved companions and especially hot subdwarf stars (Wang et al. 2021). Be stars then represent critical test beds, especially for theories that explain the role of rotation in stars, which to date remain insufficiently developed for rapid rotators.

Space photometry has led to significant advances in the field of Be stars in recent years. Analysis of space-based photometry has revealed that pulsation is ubiquitous among classical Be stars, and that they pulsate primarily in low-order g-modes where gravity is the restoring force (Rivinius et al. 2003), similar to the class of slowly pulsating B (SPB) stars (De Cat 2002). Higher-frequency p-modes (pressure being the restoring force, exemplified in the β Cephei pulsators; Stankov & Handler 2005) are also observed in some Be stars, but are less common than g-modes. Rossby waves (r-modes, where centrifugal forces dominate) may also be present in Be stars (Saio 2013; Saio et al. 2018b). Typical observed frequencies of p-modes are between approximately 6 and 15 day⁻<fnref n="1" />, and g- and r-modes are typically observed with frequencies less than 6 day⁻<fnref n="1" />. Analysis of photometry from the MOST (Walker et al. 2003), BRITE (Weiss et al. 2014), Kepler (Koch et al. 2010),

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, Rua do Matão 1226, Cidade Universitária, B-05508-900 São Paulo, SP, Brazil; jbartz@usp.br</footnote>

<footnote n="2">² LESIA, Paris Observatory, PSL University, CNRS, Sorbonne University, Université de Paris, 5 place Jules Janssen, F-92195 Meudon, France
*Received 2020 October 22; revised 2022 February 4; accepted 2022 February 7; published 2022 April 21*</footnote>
and CoRoT (Baglin et al. 2006) satellites has shown that the frequency spectra of Be stars are often complex relative to other B-type main-sequence pulsators (the β Cephei and SPB stars), typically exhibiting multiperiodicity, groups of closely spaced frequencies (as well as isolated frequencies), signatures of low-frequency stochastic variability, and long-term (usually aperiodic) trends (Walker et al. 2005a; Rivinius et al. 2016; Baade et al. 2017, 2018a; Semaan et al. 2018). The fact that virtually all Be stars pulsate suggests that pulsation is an important aspect of these systems, and is also a potentially useful probe of their interiors via asteroseismology. High-precision space photometry therefore represents a valuable tool to study the physics of Be stars and to learn about the role of rapid rotation in stellar structure and evolution.

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) mission, launched in 2018, has opened a new window into OB star variability. The Kepler spacecraft dramatically advanced the state of the art of stellar variability with its unprecedented photometric precision and long (4 yr) observational baseline of a single field (Borucki et al. 2010). However, due to its observing strategy, only a small number of relatively faint OB stars were observed. While TESS has the same general goal as Kepler (the discovery of transiting exoplanets), its observing strategy is markedly different, at great benefit to the study of OB stars. In its prime mission, TESS covered ~74% of the sky in two years, with a large field of view that shifts approximately every 27 days. Unlike Kepler, the TESS sectors have significant overlap with the galactic plane, where the vast majority of OB stars are found. Whereas Kepler observed only three known Be stars (Rivinius et al. 2016), TESS is observing over 1000. Another benefit is that the TESS mission was designed for brighter stars ($V \lesssim 12$), which means that the systems viewed by TESS have more comprehensive historical data sets and are more practical to observe from the ground.

Capitalizing on the strengths of TESS, the main goal of this work is to provide an overview of the variability seen in the population of over 400 Be stars observed in the first year of the TESS mission at a high precision and at short timescales (tens of minutes to weeks). Variability characteristics are ascribed to every star in the sample, providing insight to the behavior of the population as a whole and bringing to light patterns that exist according to spectral type and correlations between the different variability characteristics. In this context, the most typical signals are those of stellar pulsation that manifest as periodic signals in the flux of a given system. As has been reported in many studies, such periodic signals often form “frequency groups” in the observed power spectra (e.g., Walker et al. 2005a; Rivinius et al. 2016; Baade et al. 2018a; Semaan et al. 2018). These frequency groups are in general the most characteristic signature of Be stars observed with space photometry. In addition to photospheric signals, TESS is also sensitive to changes in the circumstellar environment close to the star, which can be associated with episodes of mass ejection (with the matter perhaps being inhomogeneously distributed in azimuth).

In Section 2, the TESS satellite and its data products are introduced, and methods for data extraction are described. Section 3 describes the analysis methods and shows example results of these methods for artificial light curves. Section 4 introduces the characteristic features seen in the TESS data of Be stars, which are then used to describe each star in the

sample. The results of the analysis of these signals are presented in Section 5, including discussion of each characteristic and the relevant astrophysical context. In Section 6 a broad overview is given with an emphasis on correlations between the different variability classifications presented in Section 5, followed by concluding remarks in Section 7.

## 2. Data

NASA TESS (Ricker et al. 2015) is a photometric mission performing wide-field photometry over nearly the entire sky. The four identical cameras of TESS cover a combined field of view of $24^\circ \times 96^\circ$. During the first year of TESS operations, nearly the entire southern ecliptic sky was observed in 13 sectors, with each sector being observed for 27.4 days. Some regions of the sky are observed in multiple sectors. TESS records red optical light with a wide bandpass spanning roughly 600–1000 nm, centered on the traditional Cousins I band. For optimal targets, the noise floor is approximately 60 ppm h$^{-1}$.

Full-frame images (FFIs) from TESS are available at a 30 minutes cadence for the entire field of view, allowing light curves to be extracted for all objects that fall on the detector. Certain high-priority targets were preselected by the TESS mission to be observed with 2 minutes cadence, some of which were chosen from guest investigator programs. For the sample studied in this work, light curves were extracted from the FFIs for all systems, and 2 minutes cadence light curves were also used whenever available. Generally the 30 minutes and 2 minutes light curves contain the same signals, which bolsters our confidence in the methods used for the FFI light-curve extraction. However, there are some subtle differences that are further explained in the following subsections.

Because of the large pixel size of TESS (~20"), the broad PSF, and the fact that the vast majority of OB stars lie in or near the Galactic plane, it is important to consider the effect of blending, where contaminating flux from neighboring sources can contribute to the light curve extracted for the target star. For each star in this sample, the “contamination ratio” (defined as the amount of flux from all neighboring sources divided by the flux of the target system) was calculated in the same fashion as described in Section 3.3.3 of Stassun et al. (2018) using the Python code TIC_INSPECT.<fnref n="3" /> Additionally, when analyzing each system, images from the digital sky survey were compared against the TESS images (which include the aperture that was used) to gauge potential blending issues in a more qualitative way. With moderately blended objects (where one or more neighboring stars of roughly comparable brightness contribute to the total aperture flux), different light-curve extraction methods using different aperture sizes were compared (see the following subsection), and in extreme cases, heavily blended systems (where one or more neighboring stars dominate the total aperture flux) were not included in further analysis.<fnref n="4" /> Additionally, for moderately blended cases, Simbad was queried to check for known variables in the vicinity of the target. In practice, for the majority of objects, source confusion due to blending is of little concern in the present study (except

<footnote n="3">³ https://github.com/mpaegert/tic_inspect</footnote>

<footnote n="4">⁴ e.g., TIC 434254190 = BQ Cru with a contamination ratio of 3.8173, TIC 261862960 = mu.02 Cru with a contamination ratio of 2.6170, TIC 207043035 = HD 258782 with a contamination ratio of 14.6662, TIC 380224370 = Cl* NGC 4103 with a contamination ratio of 3.8180, and TIC 709917565 = β Mon A with a contamination ratio of 2.7740.</footnote>
perhaps for the higher-frequency signals in some cases) with the main effect being a mild suppression of amplitude in the 30 minutes cadence data. The median (mean) value for the contamination ratio of the Be star sample is 4.6% (12.0%), where typically dozens to hundreds of much fainter stars contribute to this contamination (such that any one source is unlikely to contribute to the light curve in a significant way).

## 2.1. Extracting Light Curves from TESS FFIs

Light curves (LCs) are extracted from the TESS FFIs using three different methods. This is done to increase our confidence in the results when these methods produce LCs that agree, and also decreases the incidence of false-positive detections of signals if they exist in only one version of the extracted LC, being perhaps caused by imperfect removal of systematic trends or blending from neighboring objects.

The first method used to extract LCs begins by using the LIGHTKURVE package (Lightkurve Collaboration et al. 2018) and TESSCUT (Brasseur et al. 2019) to download a target pixel file with a 50 × 50 pixel grid centered on the target star's coordinates for every available TESS sector. An aperture threshold of $10\sigma$ relative to the median flux level is used as a first step to automatically determine the aperture mask for the target star. The size of the target pixel mask is allowed to scale with stellar brightness, ranging from a radius of 2 pixels for the fainter stars ($T_{mag} \ge 10$) to 6.5 pixels for the brightest ($T_{mag} \le 5$). All pixels outside of a 15 × 15 pixel exclusion zone (centered on the target star) are then used as regressors in a principal component analysis (PCA) to remove common trends across this region of the CCD (using five PCA components). The result is a PCA detrended light curve that is largely free from systematic trends. At the same time, an alternate version of the light curve is produced (second method), using only background removal, since in some instances a PCA detrending method will remove astrophysical variations in the target star with timescales of many days or longer.

The third method uses the ASTROQUERY (Ginsburg et al. 2019) routine CATALOGS to identify five stars of similar brightness on the same CCD as the target star within an annulus of an inner radius of $0.1$ and an outer radius of $0.35$ (although this is allowed to vary if there are too few stars of sufficient brightness within the original annulus). TESSCUT is then used to extract a light curve for the target star and all of the identified neighboring stars using a 3 × 3 grid of pixels. The trend filtering algorithm in the VARTOOLS light-curve analysis package (Hartman 2012) is then used to identify and remove trends that are common to the set of the target star and its neighbors. This aperture is almost always smaller than that used in the PCA method, and is therefore less susceptible to blending (at the cost of achieving a lower signal-to-noise ratio (S/N) due to excluding some flux-containing pixels).

In some cases, one of the above methods will fail partially or completely (most often when a target is close to the edge of a detector, or in the most crowded fields). Using different extraction methods allows us to determine when this is the case so that improperly reduced data is flagged as such. The first method using LIGHTKURVE tends to be more reliable and produces a light curve with a higher S/N, and all plots of 30 minutes cadence data shown in this work use this version of the light curve. The sampling rate of 30 minutes allows for the detection of frequencies up to $24 \text{ day}^{-1}$.

## 2.2. Two Minutes Cadence Light Curves

About 65% of the sample has 2 minutes cadence light curves available from the TESS office. When these LCs exist, the same analysis was performed as was done for the LCs extracted from the FFIs. In nearly all cases, the results are virtually identical between different versions of the LC. We use the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux from the TESS LCs, which is calibrated in a way that often removes long-term trends. This generally yields a cleaner version of the LC (compared to the simple aperture photometry, SAP, flux) that is well suited for frequency analysis. The disadvantage is that the detrending process may remove astrophysical variability on relatively long timescales. These, however, can easily be identified in the 30 minutes data. For that reason, both sets of data were always analyzed for all stars for which they are available. The higher cadence 2 minutes LCs allow for the detection of much higher frequencies relative to the data extracted from the FFIs, up to the Nyquist limit of $360 \text{ day}^{-1}$.

## 2.3. Sample Selection

The Be Star Spectra (BeSS) database<fnref n="5" /> (Neiner et al. 2011) was used to create a list of Be stars between V-mag of 4 and 12. The sample of Be stars from Chojnowski et al. (2015) was also included. Only targets in the southern ecliptic hemisphere were chosen, in order to match the fields observed by TESS in its first year of operations. Of this initial list of 543 stars, 64 were either not observed by TESS or the data are problematic (i.e., when a star falls onto the very edge of a detector, or systematic effects or unusually high noise levels severely hamper our ability to analyze a given LC) and therefore not analyzed further. A further 49 systems are rejected from the sample, with 35 of these being confirmed as something other than a classical Be star (see Appendix B.2), and 14 are strongly suspected of being some other type of system (see Appendix B.3).

The remaining 430 classical Be stars make up the sample studied here. Of these, 218 (51%) are of early-type (B3 and earlier), 79 (18%) are mid-type (B4, B5, and B6), 110 (26%) are late-type (B7 and later), and 23 (5%) are of unknown spectral type.

An important aspect of the analysis undertaken in this work was determining if a given object is not a classical Be star. To this end, a literature search was conducted, any available BeSS or APOGEE spectra were inspected, the spectral energy distribution (SED) was inspected,<fnref n="6" /> and any unusual features in the light curve were noted (e.g., indications of a short-period binary) for each object. Systems that masquerade as Be stars include OB stars with strong magnetic fields, chemically peculiar stars, B[e] stars, interacting and/or close binaries, Herbig Ae/Be stars, supergiants, and OB stars embedded in nebulae.

The 49 objects thus rejected from the sample suggest a “contamination fraction” of ~10%. Some of the objects rejected here have been included in works that aim to describe Be star populations (e.g., Frémat et al. 2005; Labadie-Bartz et al. 2017), underscoring the need to carefully vet samples of supposed classical Be stars against contamination from other types of astrophysical systems. This contamination fraction

<footnote n="5">⁵ http://basebe.obspm.fr</footnote>

<footnote n="6">⁶ Via the “Photometry viewer” hosted on Vizier (http://vizier.u-strasbg.fr/vizier/sed/).</footnote>
may be slightly underestimated if additional systems are
revealed to be impostor Be stars, like the 20 systems in
Appendix B.4 for which there is insufficient evidence to reject
their classical Be stars designation, despite having character-
istics that may cast some doubt as to their nature.

The majority of these targets also have years-long light
curves from the ground-based Kilodegree Extremely Little
Telescope (KELT; Pepper et al. 2007, 2012), where the data for
some of the TESS sample are published and available in
Labadie-Bartz et al. (2017, 2018). While the KELT data was
not directly used in the analysis of these systems in this work,
in some instances it was used to corroborate long-term trends or
periodic signals seen in TESS, or to confirm systems that are
not Be stars (e.g., binaries with ellipsoidal variation). Further
works studying these systems should take advantage of the
more comprehensive suite of time-series data provided by
TESS, KELT, BeSS, APOGEE, ASAS (Bernhard et al. 2018),
and perhaps other sources.

3. Analysis

A variety of methods were used to analyze the sample. Standard Fourier methods lie at the core of the analysis, but the diversity of signals expressed by Be stars demands a careful study, including visual inspection of the light curves, separately considering different frequency regimes, measuring aperiodic variability, tracing variable amplitudes of periodic signals, and documenting correlations in time between different signals. Consideration is also given to the degree of blending, and potential issues from saturation, systematic effects, and noise. Details about the data analysis are described in the following subsections.

3.1. *Removing Outliers and Bad Data*

Because of the variety of signals in the data, especially in
systems exhibiting longer-term trends of relatively high ampl-
itude, outliers cannot be removed through standard methods
such as sigma clipping. Instead, outliers and sections of poor-
quality data were removed in the following way. After
reconstructing the light curve with a sum of Fourier terms, the
calculated fit was subtracted from the observed data, and points
greater than five times the median absolute deviation from the
median of the residuals were automatically identified as outliers.
At this stage, the observed and calculated data and the residuals
were manually examined, allowing for the possibility of
removing additional sections of poor-quality data. These cleaned
light curves, with outliers and sections of poor-quality data
removed, are used in all following analysis steps.

3.2. *Low-frequency Detrending and Recovery of Periodic Signals*

Because of the generally complex nature of photometric
variability of Be stars, there are signals with a wide range of
timescales and behavior that we are interested in measuring,
which requires extra analysis steps. For example, in a system
with significant aperiodic low-frequency variability, it is
prudent to detrend against these signals when analyzing the
data for higher-frequency signals.

A Fourier-based detrending method was employed, whereby
each sector of data for each object was detrended against all
signals lower than 0.5 day⁻¹ up to a false-alarm probability of
10⁻² (in a similar fashion as Rivinius et al. 2016). This cutoff is
simply a convenient choice for removing the most significant

low-frequency signals to prepare the light curve for a frequency
analysis above 0.5 day⁻¹. This Fourier-based method has the
advantage of being applicable in a uniform fashion to every
light curve regardless of its complexity and is robust in tracing
all low-frequency behavior, whether periodic or not. Other
detrending options, such as moving averages or polynomial fits
either require some fine-tuning on a case-by-case basis (e.g., on
the order of the polynomial or size of the moving average
window), and/or can introduce spurious signals when used for
detrending (especially near the edges of data substrings for
moving averages, and when high-order polynomials are
required to trace more complex low-frequency signals). The
Fourier-based detrending method does not introduce spurious
signals at higher frequencies, which would otherwise hamper
the subsequent frequency analysis steps.

A frequency analysis was then performed for the entire set of
TESS observations for each star, separately for the original and
low-frequency detrended versions of the light curves. The
TIMESERIES.LOMBSCARGLE package (VanderPlas et al. 2012;
VanderPlas & Ivezić 2015) of ASTROPY (Astropy Collabora-
tion et al. 2013, 2018) was used to compute these modified
Fourier transforms. The above procedure was applied to both
the 30 and 2 minutes cadence data, when available.

To determine the individual frequencies present in the data,
the VARTOOLS light-curve analysis software with the LOMB-
SCARGLE routine (Press et al. 1992; Zechmeister & Kürster
2009) was used to detect and iteratively pre-whiten the data
against each recovered signal (up to a false-alarm probability of
10⁻²), recording the frequency, phase, amplitude, S/N, and
false-alarm probability. The frequency and amplitude of these
signals are used to plot the pre-whitened periodograms (with
discrete values for each detected frequency), and are also used
to reconstruct each light curve based on the recovered signals,
which are visually compared to the original photometric data to
ensure a good fit from the Fourier analysis. These methods are
similar to what is typically done for space photometry of Be
stars, e.g., with MOST (Walker et al. 2005a), Kepler (Rivinius
et al. 2016), and CoRoT (Semaan et al. 2018).

Note that the choice to stop the iterative fitting routine down
to the significance of a false-alarm probability of 10⁻² is
arbitrary but reasonable and convenient. Since the VARTOOLS
software assumes white noise and cannot be forced to
simultaneously fit the red noise profile, a “continuum” of
relatively low-amplitude frequencies is often identified by the
routine (apparent in most of the plotted frequency spectra).

3.3. *Wavelet Plots*

It is common for Be stars to show photometric signals that
vary in amplitude over time. Wavelet plots are a convenient
way to visualize this, as they depict the frequency spectrum as a
function of time (at the cost of a degraded frequency
resolution). The Python package SCALEGRAM<fnref n="7" /> was used to
perform a wavelet analysis for each star, considering the
original and low-frequency detrended (signals with frequencies
<0.5 day⁻¹ removed) versions of the data separately.

3.4. *Testing Artificial Light Curves*

In order to better understand the signals that are recovered
from TESS light curves of Be stars, tests were performed where

<footnote n="7">⁷ https://github.com/alsauve/scaleogram</footnote>
artificial signals are injected into the data and then recovered. This process begins by taking a TESS light curve of a typical Be star (TIC 14088298, $T_{mag} = 8.125$, using 2 minutes PDCSAP data) from a single sector (sector 5), and removing all signals through iterative pre-whitening, leaving a light curve that contains only noise. Signals are then injected, and the light curve is analyzed using the same methods that are applied to the real data. The injected signals were motivated by common features seen in the data for Be stars.

In the first trial, nine sinusoidal signals with different but constant amplitudes are injected into the pre-whitened light curve in two groups, with four signals centered around $2\ \text{day}^{-1}$, and five signals centered around $3.3\ \text{day}^{-1}$ to imitate the frequency groups that are commonly seen in Be stars (the amplitudes of these signals are between 0.2 and 1 ppt).

The artificial light curve is then analyzed through iterative pre-whitening. All of the injected signals are recovered, with the mean difference between the injected and recovered frequencies being 0.06%, and the recovered amplitudes are precise to within 10%. However, other signals are also found in the analysis. These spurious signals are all close to the injected frequencies and lie within their respective groups. Panel (A) of Figure 1 shows the frequency spectrum and the signals recovered through iterative pre-whitening of this first artificial light curve, with the location of the injected signals also marked. Only one frequency group is shown for clarity, as the behavior of the other group is qualitatively the same. In the following trials, only this group near $3.3\ \text{day}^{-1}$ is modified, and there is no discernible influence on regions of the frequency spectrum that are not shown in Figure 1.

Next, a single frequency at $3.5\ \text{day}^{-1}$ is injected into the light curve with an amplitude that decreases linearly from 1 ppt to reach an amplitude of 0 at the end of the light curve. Panel (B) of Figure 1 shows this trial. The injected frequency is recovered to within 0.02% with an amplitude indistinguishable from that of the injected signal at its strongest. Other signals are also found close to the injected frequency and with amplitudes decreasing with distance from this frequency. These signals are weaker than the injected single frequency, but are still significant relative to the noise level. The periodogram peak is also wider than it would be if the amplitude were not modulated.

The original nine signals are again injected, with the strongest (at $3.5\ \text{day}^{-1}$) being modulated in amplitude in the same way as the second trial. Again, there are spurious peaks near $3.5\ \text{day}^{-1}$ caused by the amplitude modulation, but now there are additional peaks in the group that did not exist in the first trial, and one of the injected signals (the second in the group, at $3.282\ \text{day}^{-1}$) is not properly recovered. Panel (C) of Figure 1 shows this.

This test is repeated, but with the amplitude of the $3.5\ \text{day}^{-1}$ signal being constant and the amplitude of the middle frequency of the group, at $3.32\ \text{day}^{-1}$, linearly decreasing at the same rate as the second and third trials. Spurious peaks appear, and the precision with which the five injected signals of the group are recovered is degraded, as shown in panel (D) of Figure 1.

A trial similar to the second one is done, but instead of a single frequency with a linearly decreasing amplitude, the amplitude is modulated by a sinusoid with a frequency of $0.0275\ \text{day}^{-1}$. In this case, the original injected frequency at $3.5\ \text{day}^{-1}$ is not recovered, but rather two strong peaks appear in the periodogram (panel (E)) of Figure 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Analysis of an artificial light curve with known signals injected. The locations of injected signals are shown with red dotted vertical lines, with amplitudes of constant signals marked by a red x, and the maximum amplitude of signals that vary in strength with time is marked by a green circle. The frequency spectrum of the artificial light curve is shown in light blue, and the signals found through iterative pre-whitening of the artificial light curve are vertical black lines, where their height reflects the amplitude of the recovered signal. The window function (with arbitrary vertical scale), centered at $3.1\ \text{day}^{-1}$, is shown in panel (B). Note that each sector will have a slightly different window function, and for systems viewed in multiple sectors the window function becomes more narrow. Panel (A): a frequency group where all injected signals have a constant amplitude. Panel (B): a single frequency with a linearly decreasing amplitude is injected. Panel (C): the same as panel (A), but the strongest signal at $3.5\ \text{day}^{-1}$ has a linearly decreasing amplitude. Panel (D): the same as the panel (A), but the middle signal at $3.32\ \text{day}^{-1}$ has a linearly decreasing amplitude. Panel (E): the same as the panel (B), but the amplitude is modulated by a sinusoid. Panel (F): the same as panel (A), but all amplitudes are smaller by a factor of ten.</figcaption></figure>

This can be explained by the trigonometric identity $\sin(\alpha) \times \sin(\beta) = [\sin(\alpha - \beta) - \cos(\alpha + \beta)]/2$. In other words, with a standard frequency analysis alone, it is impossible to distinguish between a single frequency, $\beta$, whose amplitude is modulated by some lower frequency, $\alpha$, and two signals of constant amplitude located at $\beta - \alpha$ and $\beta + \alpha$. In principle, this
can be remedied by also considering phase information. However, such an analysis is not performed in this work.

Finally, the first trial was repeated but with all amplitudes being smaller by a factor of 10 (panel (F) of Figure 1). With these amplitudes approaching the noise floor of the data, the injected signals are poorly recovered. This demonstrates that a practical, but approximate, lower limit on the amplitude of signals that can be reliably recovered in a typical TESS light curve is ~0.05 ppt.<fnref n="8" />

The first trial, where the injected signals that make up two groups all have a constant amplitude, is shown in more detail in Figure 2. The wavelet analysis of this artificial light curve shows that the two groups are variable in power over time. However, this is solely a consequence of the beating phenomenon between signals in a given group, since each injected signal is known to be constant in amplitude. The beating can also be seen as a corresponding change of the overall amplitude in the light curve itself.

Numerous other similar tests were performed, but the above examples serve to demonstrate the overall results of attempting to recover signals in groups of closely spaced frequencies when amplitudes are allowed to vary (as is often the case with Be stars). The main conclusions of these tests are that the methods used to detect signals in the TESS data will almost always result in some spurious signals whenever there are groups of frequencies and that signals with varying amplitudes can, sometimes dramatically, compound this effect. While the frequency groups themselves, and usually the strongest frequencies they comprise, are reliably recovered, some degree of spurious detections is inevitable. Complex beating patterns can exist within a frequency group, causing apparent modulation of the strength of the group over time despite all individual signals having a constant amplitude. Therefore, caution must be exercised when analyzing TESS data for Be stars, and these limitations must be kept in mind when considering the light variability seen in this sample. While these limitations are well known in general, the complex nature of Be star variability and the low-frequency resolution of most TESS light curves exacerbates such issues.

### 3.5. Interpreting Observations of Be Star Systems

Certain observed variations can be confidently attributed to either the star or the disk. In terms of photometry, brightening or fading events that occur on timescales of months or years are understood to be due to disk growth or dissipation, while coherent, stable periodic signals on timescales of around 1 day and less are best attributed to stellar pulsation. There are, however, many cases where photometric signals are ambiguous in origin. The stellar rotation period, orbital period in the close circumstellar environment, and possible pulsational periods are all very similar. Since many factors can influence the total brightness of the system (and other observables) in often complex and time-variable ways, care must be taken in interpreting photometric data. In the following, examples are provided of variability that can be firmly connected to either the disk or star, but debatable cases are also discussed, in order to

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Artificial light curve (top), modified Fourier periodogram (left), and wavelet plot (middle) for the first trial discussed in Section 3.4 with two frequency groups, with each signal having a constant amplitude. The wavelet plot shows both groups having apparently variable power, which is solely a consequence of the beating phenomenon. In the wavelet plot, horizontal dashed lines mark the injected frequencies, and regions occupied by hash marks denote gaps in data and regions where wavelet signals are otherwise unreliable due to edge effects.</figcaption></figure>

highlight both the great potential, as well as the complications ensuing from interpreting space photometry alone.

## 4. Characteristic Features of Be Stars in TESS

The main goal of this work is to assign characteristic variability features to each star in the sample, and to then use these characteristics to describe the sample as a population. This section introduces these features and shows examples.

### 4.1. Light Curves and Their Frequency Spectra

Most of the characteristic features in the TESS data of Be stars can be inferred by inspecting the light curve and the frequency spectrum. While each light curve and frequency spectrum is unique, there are certain features that are common among many members of the sample. Figure 3 shows examples of Be stars that show these characteristic features, which are introduced and described here.

In what follows, a convention is adopted where low frequencies are those less than $0.5 \text{ day}^{-1}$, mid frequencies are between $0.5 \text{ and } 6 \text{ day}^{-1}$, high frequencies are between $6 \text{ and } 15 \text{ day}^{-1}$, and very high frequencies are greater than $15 \text{ day}^{-1}$. While these distinctions are somewhat arbitrary, they are physically motivated. Low-frequency signals are generally below those of typical g-mode pulsation observed in Be stars (Rivinius et al. 2003), mid (high) frequencies span the typical range of g-mode (p-mode) pulsation in Be stars, and very high frequencies are above the typical p-mode regime in Be stars (e.g., Handler 2013; Bowman 2020, and references therein). While the rapid rotation (and evolutionary stage) of Be stars can complicate this simplified scheme, it is useful to choose

<footnote n="8">⁸ A true lower limit depends on the brightness of the target, how successfully systematic trends are removed, the details of the frequency spectrum in the vicinity of a given signal (including red noise), the number of sectors in which the target was observed, the choice of aperture used to extract the light curve, and perhaps other factors.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** TESS light curves (left) and frequency spectra (right) for a representative selection of Be stars that show certain characteristic features, as described in Section 4. For panel (B), the periodogram is re-calculated after removing the low-frequency (<0.5 day⁻¹) signals (which are shown in a lighter gray color). Panel (F) uses 2 minutes cadence data, which better emphasizes the highest-frequency signals. All other panels use 30 minutes cadence data. The frequency axis of the periodogram in panel (F) is extended to include the high frequencies. Signals at frequencies higher than 10 day⁻¹ are absent in all other stars shown here. The TIC ID, common ID, spectral type, and TESS sectors are printed in the periodogram plots.</figcaption></figure>

cutoffs to delineate these categories in order to classify the observed variability.

The characteristic features adopted are as follows:

1. *Frequency groups:* Many closely spaced frequencies often form groups in the frequency spectra of Be stars. The system in panel (A) of Figure 3 shows three well-defined groups near 0.1, 0.95, and 1.9 day⁻¹, and panel (E) shows two groups near 3 and 6 day⁻¹. Panels (B) and (F) also show two frequency groups each. There is one group in panel (D) (near 2 day⁻¹) plus stochastic variation at lower frequencies. Panel C is more ambiguous as plotted, but shows two prominent groups near 0.6 and 1.4 day⁻¹, which are more easily identified through iterative pre-whitening and acknowledging that the periodogram peaks are wider than the window function, suggesting they contain multiple unresolved signals. It is common for a harmonic series that begins with the mid-frequency groups to extend to higher frequencies.

2. *Longer-term trends dominate:* Longer-term trends with timescales >2 days are the most prominent features in the data. Panel (A) in Figure 3 is an example of this characteristic where oscillatory behavior with a timescale of ~10 days is evident (and is the highest-amplitude signal). Likewise, the high-amplitude longer-term trends

in panel (B) belong to this category. This is usually apparent from the light curve, but can be more quantitatively determined if the strongest periodogram signals are in the low-frequency regime.

(a) *Flickers:* Flickers are a specific type of longer-term trend that are characterized by an initial increase (or decrease) in brightness of a few percent over a few days, followed by a return toward baseline. The largest-amplitude features in panel (B) of Figure 3 are examples of this (e.g., the variability in the first 10 days is characteristic of a single flicker event). Flickers are not an oscillation around the mean brightness (like in panel (A)), but are rather a marked departure from the baseline brightness.

3. *Low-frequency stochastic variation:* Nonperiodic variability is a significant feature of the data. Low-frequency stochastic signals can appear as extra “noise” in the frequency spectrum that is strongest at the lowest frequencies, and decreases toward higher frequencies. However, this “red noise” is astrophysical, and arises from genuine variability. Panel (D) in Figure 3 shows an example that includes prominent low-frequency stochastic variability. The forest of signals between 0 and 2 day⁻¹ in the periodogram is stochastic in nature, while the frequency group just above 2 day⁻¹ is a distinct (periodic) feature (i.e., a frequency group). Panel (B)
likewise includes low-frequency stochastic variability (manifesting as an underlying “continuum” of signals at low frequencies in the periodogram) in addition to frequency groups that clearly stand above the local noise.

4. *Isolated frequencies:* In contrast to groups, some frequencies are isolated. There are many isolated frequencies in the periodogram of panel (F) in Figure 3 (including at 1.7, 3.4, 6.7, 13.6, and 22.5 day⁻¹), and also some in panel (D; with low amplitudes between 6.0 and 7.5 day⁻¹, and being more apparent in the 2 minutes cadence data after removing low-frequency trends).

(a) *High- and very-high-frequency signals:* Systems that exhibit periodic signals in the high- (6 < *f* < 15 day⁻¹) and very-high- (*f* > 15 day⁻¹) frequency regime are recorded.<fnref n="9" /> Panel (F) in Figure 3 shows a star with many of these high- and very-high-frequency signals, while panel (D) also meets the criteria of having high-frequency signals. Harmonics alone are not considered here, nor are near-harmonic sequences of groups that begin at mid frequencies. For example, if there is a signal at 4.0 day⁻¹ with an exact harmonic at 8.0 day⁻¹, it is not classified as a high-frequency signal since the presence of this harmonic may simply indicate that the fundamental signal (at 4.0 day⁻¹) is not perfectly sinusoidal. In panel (E), the groups near 6 day⁻¹ and 9.5 day⁻¹ are the second and third groups in a series that begins with the strongest group near 3 day⁻¹, and thus this star does not meet the criteria for having high-frequency signals.

(b) *Harmonics of isolated signals:* Some frequency spectra show clear harmonics, where a signal is found at an integer number times the frequency of another signal. In some cases these harmonics are exact, while in others they are approximate. An exact harmonic is seen in panel (F) of Figure 3, where the lowest frequency, $f_0 = 1.684$ day⁻¹, has a first harmonic at $2 \times f_0 = 3.368$ day⁻¹ (and increasingly smaller-amplitude second, third, and fourth harmonics). Frequency groups often have (multiple) harmonics, but those are considered separately.

5. *No detected signals:* A small fraction of the stars in the sample shows no variability above the TESS noise level. This is generally restricted to later spectral types, where it is well known that amplitudes are relatively low.

## 4.2. Time-variable Signals: Temporarily Enhanced Frequency Group Strength

While traditional light-curve and frequency analysis provide valuable information about the signals present in the data, it is clear from Figure 3 that some signals are variable in time. The most notable aspect of this is seen when flickers coincide with enhancement of the power of one or more frequency groups (the group enhancement may occur during or after the rising branch of the flicker). This is clearly seen, for example, in Figure 4 where both main groups (near 2 day⁻¹ and 4 day⁻¹) are strongest during the two flicker events. We limit our consideration of time-resolved signal analysis to only this situation, recording instances where groups are enhanced coincident with flickers. Further analysis is possible (e.g.,

<figure><img src="image_7.png" /><figcaption>**Figure 4.** A system where flickers correspond to the enhancement of the two main frequency groups. First: 2 minutes cadence LC (black), with the low-frequency signals in red. Second: the LC after removing low-frequency signals. Third: wavelet analysis (after removing the low frequencies). Fourth: frequency spectrum (after removing the low frequencies). The top three panels share the same x-axis (TJD).</figcaption></figure>

quantifying correlations between signal amplitudes over time) and will be explored in future works.

In principle, the alternative could be considered, where decreases in group strength are instead noted. Generally, variations in group strength most obviously manifest as increases in amplitude that are highly localized in time and thus represent a departure from their “typical” strength (averaged over the full TESS observing baseline). This is especially evident when viewing systems where groups are relatively strong during flicker events. However, the reverse (where frequency group amplitude suddenly drops beyond what is expected for beating patterns) is rarely seen.

## 4.3. Interpreting Results of Our Analysis

After identifying the characteristic features of interest, plots for the sample were visually inspected in order to determine which of the above characteristics can be attributed to each star and other information about the signals that are present (e.g., the location and relative strength of frequency groups).

There is some degree of subjectivity in assigning variability classifications. Each object was analyzed in detail independently by three authors of this work. If a consensus was reached regarding a given classification, then that classification was assigned to the system. If there was disagreement or uncertainty, then the object was inspected in more detail and a final decision was made regarding the classification in question.

<footnote n="9">⁹ Frequencies up to the Nyquist limit (360 day⁻¹ for 2 minutes cadence, and 24 day⁻¹ for 30 minutes cadence) are probed.</footnote>
<table><caption>**Table 1**
Percentages Showing Variability Classifications</caption><thead><tr><th>Variability<br>Characteristic</th><th>All</th><th>Early</th><th>Mid</th><th>Late</th><th>Unknown</th><th>All<br>(Cont.<br>Ratio &lt;5%)</th></tr><tr><td></td><td>(430)</td><td>(B3 and Earlier)<br>(218)</td><td>(B4, B5, and B6)<br>(79)</td><td>(B7 and Later)<br>(110)</td><td>(23)</td><td>(217)</td></tr></thead><tbody><tr><td>freq. groups</td><td>87.4%</td><td>91.3%</td><td>93.7%</td><td>77.3%</td><td>78.3%</td><td>88.5%</td></tr><tr><td>typical group configuration</td><td>83.5%</td><td>85.4%</td><td>85.1%</td><td>77.6%</td><td>83.3%</td><td>84.9%</td></tr><tr><td>longer-term trends dominate</td><td>34.0%</td><td>51.4%</td><td>24.1%</td><td>9.1%</td><td>21.7%</td><td>31.3%</td></tr><tr><td>flickers</td><td>17.9%</td><td>30.7%</td><td>7.6%</td><td>0.9%</td><td>13.0%</td><td>18.0%</td></tr><tr><td>flickers + enhanced freq. groups</td><td>83.1%</td><td>82.1%</td><td>100.0%</td><td>100.0%</td><td>66.7%</td><td>82.1%</td></tr><tr><td>low-frequency stochastic</td><td>24.7%</td><td>33.0%</td><td>12.7%</td><td>15.5%</td><td>30.4%</td><td>23.5%</td></tr><tr><td>isolated freqs.</td><td>29.3%</td><td>22.5%</td><td>34.2%</td><td>36.4%</td><td>43.5%</td><td>26.3%</td></tr><tr><td>harmonics of isolated freqs.</td><td>21.4%</td><td>6.1%</td><td>25.9%</td><td>30.0%</td><td>50.0%</td><td>22.8%</td></tr><tr><td>high freq.</td><td>14.7%</td><td>15.6%</td><td>16.5%</td><td>12.7%</td><td>8.7%</td><td>10.6%</td></tr><tr><td>very high freq.</td><td>3.0%</td><td>2.3%</td><td>2.5%</td><td>5.5%</td><td>0.0%</td><td>2.8%</td></tr></tbody></table>

**Note.** Fraction of stars showing each type of variability, according to their spectral type, followed by the total number of systems showing the given characteristic. The percentage of systems with flickers + enhanced frequency groups are calculated from only the systems with frequency groups. The percentage of systems with a typical group configuration are calculated from only those with frequency groups, while the percentage with harmonics of isolated frequencies are calculated from only those with isolated frequencies. The final column shows the variability fractions for only the half of the sample with the smallest contamination ratios (<5%), which are least susceptible to blending issues.

In some cases, it remained not possible to confirm or reject a given classification.

Panels (E) and (F) of Figure 3 show an apparent modulation of the amplitude of the most prominent signals. In these and all similar cases, such modulation is interpreted as a beat envelope from two or more closely spaced signals even if the observational baseline is too short to fully resolve these signals (i.e., a frequency group). The alternative view that a single frequency is varying in amplitude on these relatively short timescales is not employed. Certain other Be stars studied with much longer baseline space photometry (e.g., $StH\alpha$ 166 with 4 yr of Kepler data; Rivinius et al. 2016) show closely spaced groups with beat envelopes that have timescales of tens to hundreds of days. Although TESS may not fully resolve each signal in such a group, the apparent presence of a beat envelope is a reasonable indicator that multiple signals exist.

### 4.4. Determining the Center of Frequency Groups

A Python routine was developed to more objectively quantify groups and their properties from the pre-whitened frequency spectrum for each star. This clustering algorithm served to identify groups, and to provide numerical descriptions of their net amplitude and weighted center, thus allowing the location and relative strength of groups identified in a given light curve to be compared. Further details are provided in Appendix A.

## 5. Results and Discussion for Each Variability Type

The variability features for all Be stars in the sample were tabulated and the occurrence rates are shown in Table 1. The last columns in Table 1 are computed from the ~half of the sample with the lowest contamination ratios (<5%), which are largely free from any blending issues. This subsample has virtually the same distribution of spectral type bins as the full sample. In the following subsections, the variability types and characteristic examples are discussed. These results are discussed in a broader scope in Section 6.

### 5.1. Frequency Groups

#### 5.1.1. Overview

The existence of frequency groups in the power spectra is a common feature of Be stars observed from space (Walker et al. 2005a; Saio et al. 2007; Saio 2013; Kurtz et al. 2015; Baade et al. 2018a; Semaan et al. 2018). Eighty-seven percent (376/430) of this sample shows one or more frequency groups. According to spectral type, this percentage is 91% (199/218) for early-types (B3 and earlier), 94% (74/79) for mid-types (B4–B6), and 77% (85/110) for late-types (B7 and later).

Frequency groups in classical Be stars are clearly complex in their nature. The majority of Be stars exhibit two or more frequency groups. However, a not insignificant fraction of these systems (116, or 27%) shows only one or zero frequency groups (or have groups in an atypical configuration, which is not following an approximate harmonic series). This is contingent on the noise level (perhaps there are low-amplitude frequency groups below the detection threshold), and it is also sometimes observed that the strength of groups and/or their most prominent frequencies can significantly vary over time (Smith et al. 2006; Labadie-Bartz & Carciofi 2020a; Borre et al. 2020; Labadie-Bartz et al. 2021), so that the lack of observed groups at the time when a given star is observed by TESS does not necessarily mean that groups have not existed in the past (or will not become more prominent in the future).

Some Be stars observed with sufficiently long baselines demonstrate that it is possible for signals to be too closely spaced to resolve in just one or a few TESS sectors. For example, in one of the Be stars observed by TESS in all 13 sectors of Cycle 1, TIC 279430029 = HD 53048, it is evident from an analysis of the full year-long light curve that there are frequency groups. Yet, these signals are so closely spaced that they are completely unresolved in a single TESS sector with no hint of a beat envelope. Therefore it is possible that some objects that are here determined to have only isolated frequencies (and not groups) may instead host very narrow frequency groups.
The presence of frequency groups is the most common characteristic signal seen in the Be stars of this sample. Balona & Ozuyar (2021) analyzed TESS data for a sample of 441 classical Be stars (from sectors 1–26), and found that only 73% show short-period variability, which is considerably lower than even the fraction of stars in the present sample showing frequency groups. Of the 15 Be stars in the study of Semaan et al. (2018), 12 or 13 (80%–87%) were found to have frequency groups in CoRoT photometry (their sample skews more heavily toward early-types compared to this work). In the 15 γ Cas analogs (which are early-type classical Be stars) studied with TESS photometry in Nazé et al. (2020b), they report frequency groups in 10 systems, but it could be argued that all but two exhibit groups. Indeed, there is a consensus that frequency groups are very common in Be stars regardless of any differences in interpretation.

Multiple frequency groups like those seen in the majority of Be stars are not common in non-Be OB stars (e.g., Bowman et al. 2020; Burssens et al. 2020). Frequency groups are often a consequence of very rapid rotation in pulsators (Saio et al. 2018b; Lee & Saio 2020), and may be an important component of the Be phenomenon (where the Be phenomenon is understood to be the active ejection of material and not simply the presence of a disk, which may have been formed at some time in the past and does not necessarily imply ongoing mass ejection).

There are multiple processes known to produce frequency groups in rapid rotators. g-mode pulsation can naturally form frequency groups, as modes with different radial orders ($n$), but the same value of $l$ and $m$, will oscillate at similar but slightly different frequencies. For sectoral g-modes (i.e., $l = |m|$ and no latitudinal nodal lines exist, which are common in Be stars; Maintz et al. 2003; Rivinius et al. 2003; Neiner et al. 2012a), the frequency in the inertial (observational) frame is proportional to $|m|$ and the first group is generally found above the rotation frequency by about 20% (Cameron et al. 2008; Saio et al. 2018a, 2018b; Semaan et al. 2018). r-modes in moderate to rapidly rotating stars also form frequency groups, where the strongest of these r-mode groups ($|m| = 1$) is located at slightly below the rotational frequency, with further groups becoming weaker in observed amplitude with increasing $|m|$, and located slightly below $|m| \times f_{\text{rot}}$ (Saio et al. 2018b). Variability in the circumstellar environment, especially during active mass ejection, can cause frequency groups near and slightly below the rotational frequency, and also generate harmonics (Štefl et al. 1998; Baade et al. 2016). The characteristic frequencies of this phenomenon are dictated by the orbital timescale close to the star, which depends on the orbital radius of the recently ejected material (which generally evolves with time). Balona & Ozuyar (2021) have interpreted frequency groups in Be stars as reflecting rotation of an inhomogeneous surface where magnetic fields are invoked. In general, detailed modeling and/or spectroscopic observations are needed to distinguish between these scenarios, any or all of which may occur in a given Be star over some observational baseline.

About 10% of the sample shows one or more frequency groups that are qualitatively similar to the groups of r-modes reported in Saio et al. (2018b), with a somewhat broad “hump” and a relatively strong and narrow signal on the high-frequency side of the hump. Panel (7) in Figure 5 shows one example, but TIC IDs 237059039 and 295099096 exhibit the greatest similarity to these “hump and spike” r-mode groups. Figures 3

and 9 of Semaan et al. (2018; for the Be stars CoRoT 102686433, 102785480) show frequency spectra that are perhaps consistent with r-modes. While frequency groups are very common, and despite predictions that r-modes are expected in rapid rotators (Saio et al. 2018b), “hump and spike” configurations consistent with r-modes seem to be the exception and not the rule. If r-modes are indeed ubiquitous in Be stars, they may be somewhat hidden in the relatively complex frequency groups, such that disentangling the contributions from r-modes, g-modes, and/or circumstellar variations is a challenge. With contemporaneous time-series spectroscopy and/or polarimetry, it may be possible to distinguish between r-modes and circumstellar signals, which are both expected to form groups slightly below the stellar rotation frequency.

### 5.1.2. Typical Group Configurations and Relative Strength

Although there is a wide range in the location, width, number, and relative strength of individual frequencies making up the group, and relative strength of the groups themselves, there are some patterns that are common when considering the whole sample. The most typical configuration includes three groups (although there may be further harmonics of these with decreasing amplitude). The lowest-frequency group, $g0$, is centered at $<0.5 \text{ day}^{-1}$ (and often much lower, as is more clearly seen in many instances in similar studies with longer time baselines; e.g., Neiner et al. 2012a; Baade et al. 2016, 2017; Semaan et al. 2018). The next group, $g1$, is centered at some intermediate frequency (typically between 0.5 and $3 \text{ day}^{-1}$), and $g2$ is located at approximately twice the frequency of $g1$. A variation of this configuration is seen when $g0$ is absent, but $g1$ and $g2$ still follow the same pattern. Figure 5 shows many examples of frequency groups, where all but panel (10) exhibit this typical configuration. In what follows, the center of a typical frequency group, as determined by the clustering algorithm of Appendix A, is written as $f_{g1}$ or $f_{g2}$. Note that in some systems, an otherwise typical near-harmonic series of groups may skip what would be the third group in the series (i.e., at $\sim 3 \times f_{g1}$) as in panels (7) and (9) of Figure 5.

The frequency spectra of this sample have varying degrees of complexity, and there are some cases where multiple signals exist in the vicinity of each other. Panel (8) of Figure 5 shows one such example, where the frequency spectrum does include two groups following the typical pattern ($f_{g1} \approx 2.45 \text{ day}^{-1}$, and $f_{g2} \approx 4.94 \text{ day}^{-1}$). There is also a strong isolated signal ($f \approx 1.30 \text{ day}^{-1}$) and its subharmonic, which are seemingly unrelated to the two frequency groups. Situations like this may reflect “composite frequency spectra,” where the $1.30 \text{ day}^{-1}$ signal and its subharmonic may indicate binarity or rotation (perhaps of a companion star or system), while the pair of higher-frequency groups is formed by families of pulsation typical of Be stars. This example is purely speculative—such interpretation requires further data and analysis and is beyond the scope of this work.

All cases where the typical group configuration is realized (and when it is not) and the relative strength and central location of $g1$ and $g2$ were recorded. No attempt was made to specify the center of $g0$ because the frequency resolution in this regime is poor, owing to the short TESS observing baseline, except that, per the adopted definition, $g0$ is centered at $<0.5 \text{ day}^{-1}$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** A selection of frequency spectra showing various permutations of group configurations and morphology. Solid colored vertical lines are the signals recovered with iterative pre-whitening, and the solid black line is the modified Fourier periodogram without any pre-whitening. The inverted solid lines show the same frequencies multiplied by two, and below that, the dashed lines show the frequencies multiplied by three to aid in visualizing the locations of the first and second harmonics of the recovered signals. The center of the two main frequency groups ($g1$ and $g2$) are marked in each plot with the symbol size indicating which is the stronger of the two (or with equal symbol sizes when both groups are similar in strength). Panel (1): closely spaced, narrow groups monotonically decreasing in strength. Panel (2): $g1$ and $g2$ have similar strengths, without further harmonics. Panel (3): wider groups, with signals corresponding to the second harmonic of $g1$, but without significant signals corresponding to the first harmonic of $g2$ (or the third harmonic of $g1$). Panel (4): $g1$ and $g2$ are wide, and the region beyond $g2$ is populated with signals that form a "continuum" corresponding to harmonics of signals in $g1$ and $g2$. Panel (5): $g1$ is narrow and dominated by a single signal, and $g2$ is relatively wide and centered at a frequency less than $2 \times g1$. Panel (6): $g2$ is centered at greater than $2 \times g1$, and there are also signals located in the region occupied by the second harmonic of $g1$. Panel (7): $g1$ and $g2$ are narrow and mirror each other in structure, but there are only signals that correspond to $2 \times g2$ and not $3 \times g1$. Panel (8): two fairly typical groups, $g1$ and $g2$, being centered near 2.5 and $5 \times 10^{-1}$, plus a strong apparently isolated signal at $1.30 \times 10^{-1}$ and its weaker subharmonic at $0.65 \times 10^{-1}$. Panel (9): widely separated groups, where $g2$ is stronger than $g1$, and with power located in the region corresponding to $2 \times g2$ and not $3 \times g1$. Panel (10): only one frequency group exists, along with weak low-frequency signals.</figcaption></figure>

Among the systems with frequency groups, 84% (314/376) show the typical configuration with a well-defined $g1$ and $g2$ (that is, both groups are prominent above the noise level, and the location of $g2$ is approximately at twice the frequency of $g1$), with this fraction decreasing slightly from early to late spectral types (85%–170/199, 85%–63/74, and 78%–66/85, respectively, but this decrease is likely impacted by instrumental sensitivity, as amplitudes are low for later spectral types). When a frequency spectrum does include groups, but does not follow the typical configuration, there is often just one group, or, less commonly, there are multiple groups that are clearly not harmonically related (e.g., the second group is not located near twice the frequency of the first group). Panel (10)

of Figure 5 is an example of an atypical group configuration, with a single group centered at $2.6 \times 10^{-1}$ without a second group. The lack of a typical group configuration (including those with no groups at all) is not evidence against a given system being (or including) a classical Be star (e.g., Rivinius et al. 2020). In both Be stars located slightly beyond the main sequence in Semaan et al. (2018), no groups exist and only low-frequency stochastic variability is seen. There may also be correlations between the separation between groups and the fractional main-sequence age of a star (Semaan et al. 2018). The TESS data provide an excellent opportunity to test these ideas, since these brighter stars are more amenable to spectroscopic studies from which their stellar properties can
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The center of the two main frequency groups is shown for each system with a typical frequency group configuration. The stars are ordered according to the location of $f_{g1}$, increasing upwards. Black (red) dots mark $f_{g1}$ ($f_{g2}$), and green triangles mark $f_{g2} - f_{g1}$. The gray lines are centered at $2 \times (f_{g1} \pm 5\%)$. In the main panel, marker size is proportional to the relative group strength.</figcaption></figure>

be measured (including age, rotation rate, effective temperature, and mass).

For all systems with groups in the typical configuration, the weighted central frequency of the two primary groups ($f_{g1}$ and $f_{g2}$), and their differences, are shown in Figure 6, organized according to $f_{g1}$. Figure 7 also presents information about these groups, showing the ratio of the central frequency of g2 and g1 ($f_{g2}/f_{g1}$), and the fraction of systems having groups of similar strength ($A_{g1} \approx A_{g2}$), g1 being stronger than g2 ($A_{g1} > A_{g2}$), or g2 being stronger than g1 ($A_{g1} < A_{g2}$).

From these plots, it is clear that the ratio $f_{g2}/f_{g1}$ is centered at 2.0 (or equivalently, $f_{g2} - f_{g1} \approx f_{g1}$), having a higher (lower) scatter for earlier (later) spectral subtypes. The distribution of the relative strength of g2 and g1 differs across the range of spectral subtypes, where later-type stars are preferentially found to have a relatively strong g1 ($A_{g1} > A_{g2}$), with only 7% having $A_{g1} < A_{g2}$. A similar, but more mild trend is seen in the mid-type stars, and there is a roughly flat distribution for early-type stars. Among the 71 systems with flickers and a typical group configuration, a different trend is seen where only 12 (17%) have $A_{g1} > A_{g2}$. Comparing the strength of g1 and g2 is somewhat subjective, but this decision is aided by the algorithm used to identify groups and their central frequency and relative strengths described in Section 4.4 and Appendix A.

## 5.2. Light Curves Dominated by Longer-term Variability

The light curves of 34% (146/430) of this sample are dominated by longer-term trends ($f < 0.5$ day$^{-1}$, or characteristic timescales >2 days), which may or may not be periodic or cyclic in nature. Like most other types of variability, the light curves of early-type stars are more likely to display this behavior (51%; 112/218) compared to mid- (24%; 19/79) and late- (9%; 10/110) types. Flickers stand out as the most remarkable type of behavior that causes a light curve to be dominated by longer-term trends, and in these cases, the low frequencies are (at least in part) tied to occasional episodes of mass ejection.

However, not all systems with longer-term trends exhibit flickers, and systems in which low frequencies dominate are

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Top: histogram of the ratio of the weighted central frequency of g2 and g1 ($f_{g2}/f_{g1}$). Bottom: comparison of the relative strength of g1 and g2, with the fraction of each spectral subtype bin falling into the three categories printed above its respective bar in the plot. A comparison of the relative strengths of g1 and g2 is also made for all systems with flickers (the majority being early-type).</figcaption></figure>

not necessarily ejecting mass during TESS observations. On the other hand, longer-term trends may reflect a variable density in the innermost disk (i.e., from variable mass ejection rates) even if the photometric signal(s) do not follow the typical pattern of a flicker event. These slow signals can also reflect pulsation—in particular the nonlinear coupling of two or more pulsation modes, as a simple beating pattern of two frequencies will not produce power at their beating period in a frequency spectrum (Kurtz et al. 2015; Baade et al. 2018b).
In the notation of this study and other works, often multiple
signals form a low-frequency group ($g0$) below $0.5$ day$^{-1}$. In
space photometry with longer observational baselines (e.g.,
BRITE, CoRoT, Kepler, and SMEI), it is sometimes found that
one or more of the signals in $g0$ are related to prominent signals
in $g1$ through difference frequencies. That is, a signal in $g0$ is
found at the difference between two signals in $g1$. Since there
are often many frequencies in $g1$, $g0$ can be populated by
multiple difference frequencies. Unfortunately, the TESS
observational baseline is usually too short to confidently
measure the location of signals in $g0$; so in the present work, no
attempt is made to relate the frequencies in $g0$ to those in $g1$.

Other studies have found photometric variability in this
range of timescales (days to weeks) at a similar incidence rate.
In their sample of 15 Be stars observed with CoRoT (10 of
which are B3 and earlier), Semaan et al. (2018) found
variability with timescales of days and longer in nine systems,
although in two of these (their stars number 17 and 22), the
outburst amplitudes are very low and would not have been
considered as having substantial longer-term trends (or flickers)
based on the criteria of the present analysis of TESS data.
Balona & Ozuyar (2020) analyzed the TESS light curves of 57
TESS Be stars, but use the PDCSAP flux, which in many cases
removes longer-term trends, such that comparisons cannot be
made in regards to the longer-term variability. Ground-based
photometric studies of Be star variability generally have much
longer observational baselines compared to TESS, such that
"long-term variability" (seen from the ground) may have
timescales of hundreds of days or longer but with low rates of
change such that these systems may not appear significantly
variable in TESS, and it is therefore difficult to make direct
comparisons to this study. Nevertheless, a sample of ~500 Be
stars observed with the KELT survey showed that 38% of Be
stars (and 51% of early spectral types) exhibited periodic or
cyclic variability with periods between 2 and 100 days
(Labadie-Bartz et al. 2017), and studies with both KELT and
ASAS photometry (Bernhard et al. 2018) show a large number
of systems with outbursts that occur at a high rate (up to 20
outbursts per year). Systems similar to those that exhibit the
aforementioned behavior are reasonably likely to be variable on
timescales of days and longer in TESS.

5.3. Flickers

5.3.1. Background and Overview

Outbursts (discrete episodes of mass ejection) are a well-
known and common feature of Be stars, which manifest in
photometry as an increase in brightness as the growing quantity
of circumstellar material emits and reprocesses stellar light
(Carciofi et al. 2012; Haubois et al. 2012; Sigut & Patel 2013;
Labadie-Bartz et al. 2018). In the case of systems viewed at
high inclination angles (i ≳ 75°), the growth of a disk instead
causes a net dimming, as the relatively cool circumstellar
material obscures the stellar photosphere. The timescales of
outbursts range from days to many years or even longer
(Rivinius et al. 2016; Bernhard et al. 2018; Ghoreyshi et al.
2018; Labadie-Bartz et al. 2018; Rímulo et al. 2018; Semaan
et al. 2018). The term “flicker” can be used to describe
outbursts with short timescales (days to weeks; Keller et al.
2002). With the relatively short baseline of TESS light
curves, it is only possible to probe these flickers on short
timescales. Due mostly to a lack of a sufficiently large sample

with high-precision near-continuous photometry, flickers with
timescales of days have so far been poorly studied.

As exemplified in the Be star HD 49330 observed with CoRoT (Huat et al. 2009), a typical flicker event can be broken up into four phases: the relative quiescence phase (where the light curve is flat), the precursor phase (where the brightness gradually and mildly decreases), the outburst phase (where the brightness increases), and a relaxation phase (where the brightness returns toward the baseline level). These phases also have spectroscopic counterparts (e.g., Rivinius et al. 1998a). Precursor phases are not always seen.

5.3.2. Main Results from TESS

Flickers were found in 18% (77/430) of the sample. There is a strong dependence on spectral type, with 31% (67/218) of early systems, and 8% (6/79) and 1% (1/110) of mid- and late-type systems showing flickers. This is consistent with other observational studies that find early Be stars to be significantly more variable and active, and with relatively large disk densities (and therefore higher-amplitude observables; Vieira et al. 2017; Labadie-Bartz et al. 2018). Figure 8 shows the light curve and wavelet plot for five examples of systems with flickers.

In this work, the vast majority of flickers are identified by a
net brightening. Dimming flickers in shell stars are relatively
difficult to confirm with photometry alone on the short
timescales of TESS observations, largely because of the
relatively small amplitude and naturally longer timescale of
the change in optical continuum flux compared to equivalent
events viewed at lower inclination angles (Haubois et al. 2012).
Nevertheless, there are a few instances in this sample where
dimming flickers are apparent. The second panel of Figure 8
shows one example.

Systems viewed at an inclination angle of i ~ 70° may exhibit no detectable change in their optical flux associated with the growth or dissipation of a disk-optical, and near-IR photometry is largely blind to disk events at this intermediate inclination angle (Haubois et al. 2012). Some Be stars build up disks over timescales much longer than a TESS observing sector (e.g., Labadie-Bartz et al. 2017; Rímulo et al. 2018), where the rate of change in brightness is gradual and difficult to detect in TESS. Mass loss happening on top of a strong preexisting disk will change the photometric excess only if it has sufficient strength and duration to alter the size of the continuum emitting region (Vieira et al. 2015). Weak and/or too short mass-loss events can occur, which may still contribute to the disk mass budget, but without a significant detectable photometric signature. For these (and possibly other) reasons, the incidence rate of flickers in our sample is underestimated, and a lack of flickers in a given system does not imply the absence of ongoing mass ejection.

5.3.3. Photometric Flickers as Tracers of Mass Ejection

Figure 9 demonstrates an example where the TESS bright-
ness and Hα line profile evolve together during a flicker for the
Be star HD 194779, which was observed multiple times
spectroscopically over the TESS observing baseline. In this
system, the Hα line indicates a transition from a disk-less to a
disk-possessing state. This transition unambiguously indicates
the ejection of material from the star into the circumstellar
environment, which coincides with a flicker as seen in the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Examples of systems with flickers. Below each light curve is the wavelet plot showing the evolution of the frequency spectrum over time after subtracting the low-frequency signals (the red curve), mainly demonstrating how one or both of the frequency groups typically reach higher amplitudes during flickers.</figcaption></figure>

TESS photometry. Although this system was observed in TESS Cycle 2, and is thus not included in the present sample, it serves to illustrate the connection between a photometric flicker and the ejection of stellar material. A spectroscopic campaign including HD 194779 and many other Be stars contemporaneous with TESS observations is ongoing and will be the subject of forthcoming works, with preliminary results reinforcing the interpretation of photometric flickers as mass ejection events.

5.3.4. Flickers and Frequency Groups

Other works have presented evidence of the frequency
spectrum changing before, during, and after outbursts (e.g.,
Huat et al. 2009; Baade et al. 2018a; Semaan et al. 2018). This

is also seen in the majority of flickering Be stars in this sample,
as there are many cases where flickers roughly coincide with
the emergence or enhancement of frequency groups; about 83%
(64/77) of stars with flickers show this behavior. Under the
interpretation that flickers correspond to mass ejection events,
there is some ambiguity in interpreting the enhancement of
frequency groups in these instances. Frequencies about 10%
lower than the dominant pulsation frequency were first
discovered in the line profiles of Be stars actively ejecting
material (and later seen in photometry of active Be stars), with
the preferred interpretation being that they arise in the
circumstellar environment due to an inhomogeneous distribu-
tion of recently ejected material orbiting the star (Štefl
frequencies; Štefl et al. 1998; Baade et al. 2016). The
enhancement of the overall strength of a frequency group
during times of outburst may then be related to increased
pulsational amplitude and/or circumstellar variability. Circum-
stellar signals are expected to be stronger in systems viewed at
higher inclination angles (Baade et al. 2016).

All systems exhibiting flickers also have one or more
frequency groups, suggesting a strong link between these
features. In addition to the net change in brightness, which is
the hallmark of a flicker, one or more frequency groups are
often (but not always) enhanced near in time to the flicker
event. Figure 8 demonstrates examples of this common trait. In
systems with typical group configurations and flickers, it is
more common for g2 to be more strongly variable than g1.

The temporary in-phase superposition of two or more modes
in a frequency group has been in some instances found to
correspond to mass ejection episodes, for example, in the Be
stars µ Cen (Rivinius et al. 1998b), StHα 166 (Rivinius et al.
2016), η Cen (Baade et al. 2016), 28 Cygni (Baade et al.
2018a), and 25 Ori (Baade et al. 2018b). In some of these cases
(25 Ori, µ, and η Cen), the resultant amplitude is greater than
the sum of the base mode amplitudes, perhaps indicating
nonlinear amplification (as opposed to beating). The enhance-
ment in frequency group strength seen to commonly coincide
with flickers in the TESS data may reflect this; i.e., nonlinear
amplification (or in the simplest case, linear beating) in the
amplitude of frequency groups may be triggering some of the
flickers. However, the exact timing of frequency group
enhancement may happen during the rising branch of a flicker
(e.g., Figure 4 and the second and fourth panels of Figure 8),
may happen near the time of peak brightness (e.g., the first
panel of Figure 8), or may be more ambiguous (e.g., in the third
and fifth panels of Figure 8).

All six of the stars with flickers in the 15 star sample of Semaan et al. (2018) also have frequency groups, and for these stars, g1 and g2 (using our notation) are relatively strong during the brightening phase (or the initial fading phase in the case of shell stars); but in one case (their star number 6), the amplitude of g1 is relatively low during the precursor phase. On the other hand, in some cases, a decrease in amplitude of frequency groups and/or isolated frequencies are reported to coincide with flicker events. In HD 49330, Huat et al. (2009) reported a decrease in amplitude of two isolated signals near 5.03 and 16.89 day⁻¹ during the rising branch of the flicker (while at the same time the two lower-frequency groups increase in amplitude). In the TESS sample, we do not find any convincing cases where the amplitude of any frequency groups decrease during some phase of a flicker (beyond what is expected from
<figure><img src="image_12.png" /><figcaption>**Figure 9.** An example of a star with a typical flicker that was monitored spectroscopically during TESS observations. A portion of the TESS light curve is shown (top-left panel), along with the Hα line from the spectroscopic observations (right panel). The colored triangles in the top-left panel indicate the spectroscopic epochs, with the corresponding line profiles in the right panel (following the same color code). Hα emission begins to appear near BJD-2457000 = 1720, indicating that material is being ejected into the circumstellar environment. The Hα equivalent width (EW) is plotted in the bottom-left panel, as measured between the vertical dashed lines in the right panel (with smaller values corresponding to more emission). Large symbol sizes in the bottom-left panel correspond to the thicker lines in the right panel to highlight spectra taken near the beginning of the flicker, near peak brightness, and late in the dissipation phase.</figcaption></figure>

beating for groups with many frequencies). Likewise, we do not see variations in the amplitude of isolated frequencies coincident with flickers. This is not to say that isolated frequencies cannot vary in amplitude, especially as the evidence for this in HD 49330 is convincing. Instead, this may be due in part to a relative dearth of flicker stars that have isolated frequencies in our sample, and also our analysis methods, which did not explicitly search for this phenomenon.

## 5.4. Low-frequency Stochastic Variability

### 5.4.1. Overview and Main Results from TESS

Low-frequency stochastic variability is sometimes a prominent feature in the frequency spectra of Be stars, as has been pointed out for some Be stars observed from space (e.g., *η* Cen, where the star is apparently nearly continually losing mass; Baade et al. 2016). Twenty-five percent (106/430) of this sample has low-frequency stochastic variability as a prominent feature. Low-frequency stochastic variability is noticeable among early (33%, 72/218), mid (13%, 10/79), and late (15%; 17/110) spectral types, with amplitudes generally being highest in early-type stars. Figure 10 shows six systems where low-frequency stochastic variability is a prominent feature.

In this work, the designation of low-frequency stochastic variability is applied qualitatively. An important difference between this work and the work and terminology of, e.g., Bowman et al. (2019a, 2019b, 2020), is that a system is considered to show low-frequency stochastic variability only if it is obviously a dominant aspect of its light curve and frequency spectrum (viewed in a linear scale). Otherwise, without more in-depth analysis, this label becomes meaningless as virtually every star in this sample is expected to have some degree of low-frequency stochastic variability (i.e., astrophysical correlated red noise) as it is defined in other works, such as Bowman et al. (2019b, 2020).

A proper analysis and parameterization of the stochastic low-frequency excess of the TESS Be star sample can and should be undertaken and compared to the results for the more slowly rotating OB star sample of Bowman et al. (2020). However, so

far, there at least do not seem to be any glaring differences between the red noise properties of Be stars compared to more slowly rotating B-type dwarfs. For instance, the eight early-type Be stars in Nazé et al. (2020b) with measured red noise parameters all fit into the range of parameter values determined for the B dwarfs in Bowman et al. (2020). Large samples seem necessary to reveal any correlations since there is considerable scatter in the fitted red noise parameters of OB stars, but there are already hints that evolutionary status may shape the red noise parameters to some degree. The two post-main-sequence Be stars analyzed in Semaan et al. (2018) are dominated by low-frequency stochastic variability and are without coherent periodic signals, and in general more evolved stars have higher red noise amplitudes and probably also longer characteristic timescales (Bowman et al. 2020). As a further point of comparison, the massive evolved Wolf-Rayet (WR) stars and luminous blue variables (LBVs) also have similar red noise properties compared to the main-sequence OB stars, but perhaps with slightly higher amplitudes and longer timescales (Nazé et al. 2021).

### 5.4.2. Underlying Mechanism(s) and Relevance to the Be Phenomenon

Low-frequency stochastic variability is a ubiquitous feature of OB stars, with a range in amplitude (typically between 0.05 and 5 ppt) and characteristic timescale (typically between 0.2 and 2 days), with more evolved massive stars typically having higher amplitudes and longer timescales (e.g., Bowman et al. 2019b, 2020). There are many physical processes that can plausibly give rise to low-frequency stochastic variability in classical Be stars, including inhomogeneities in the stellar surface or wind combined with rotation (Moffat 2008; Aerts et al. 2018; Simón-Díaz et al. 2018), subsurface convection layers (Cantiello et al. 2009, 2011; Cantiello & Braithwaite 2011, 2019; Cantiello et al. 2021), and internal gravity waves generated at the interface between the convective core and the radiative envelope, which propagate outwards (Rogers et al. 2013; Bowman et al. 2019a, 2019b; Edelmann et al. 2019; Horst et al. 2020). Although low-frequency stochastic
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Light curves (left) and frequency spectra (right) of a few example systems exhibiting low-frequency stochastic variability. In the frequency spectra plots, the lighter gray curve is computed from the raw data (shown in the light-curve plots), and the black curve is computed after removing signals with frequencies below $0.5 \text{ day}^{-1}$.</figcaption></figure>

variability does not contain any coherent periodic signals, there is still valuable diagnostic potential in the profile and amplitude of the frequency spectrum of these stochastic signals. For example, Bowman et al. (2020) analyzed TESS data for 70 OB stars with TESS, concluding that there is strong evidence for internal gravity waves by comparing measurements of the stochastic frequency spectrum to models of wave propagation in stellar interiors, which originate at the interface of the convective core and radiative envelope.

These internal gravity waves may be an important aspect of Be stars, since they are efficient at transporting angular momentum from the stellar interior outward (Rogers et al. 2013; Edelmann et al. 2019). While it is beyond the scope of this work to quantitatively analyze the stochastic variability seen in the TESS Be star sample, it is important to note that these features, although appearing random and incoherent, can still be used as a diagnostic to infer the physical origin(s) of the signals.

Furthermore, internal gravity waves can drive g-mode pulsation in Be star envelopes. This is especially important in early-type Be stars, which are generally too hot for the $\kappa$ mechanism to excite g-mode pulsations (Dziembowski et al. 1993; Neiner et al. 2012b, 2020). However, even in the hottest Be stars, the $\kappa$ mechanism may still act in some limited capacity to drive g-modes if certain conditions are met (Pamyatnykh 1999; Neiner et al. 2012b), such as enhanced metallicity or if the star is evolved to near the terminal age main sequence (TAMS), but this alone cannot explain the observed

distribution of frequency groups among the earliest stars in this sample.

Since virtually all of the Be stars with periodicity have prominent signals in the traditional low-frequency g-mode regime<fnref n="10" /> (including the early spectral types; see Figure 6), the excitation of these modes by internal gravity waves may even be fundamental to the Be phenomenon in these systems. As outlined in Neiner et al. (2020), internal gravity waves can serve the dual purpose of transporting angular momentum to the envelope (causing it to spin up and lowering the barrier to mass ejection) and can also excite groups of g-mode pulsation, where the activity and processes in the surface layers may then meet the conditions required to trigger an outburst whereby mass and angular momentum is transported to the disk and ultimately out of the system.

## 5.5. Isolated Signals and Possible Harmonics

Many of the frequency spectra (29%, 126/430) of this sample contain isolated, individual signals that do not obviously belong to a group. These signals favor cooler stars, as 23% (49/218), 34% (27/79), and 36% (40/110) of early-, mid-, and late-type stars, respectively, show these isolated signals. Again, this may be related to instrumental sensitivity, where a frequency group may exist, but only the strongest signal rises above the noise level (an effect expected to be more

<footnote n="10">¹⁰ Exceptions include TIC 216158265, 216875138, 338738989 (see the Appendix B).</footnote>
pronounced in later-type stars with lower intrinsic amplitudes). The prototypical Be star, γ Cas, despite having a spectral type of B0.5 IVe and being intensively observed for well over 100 yr, suffered from this observational bias where only isolated signals were identified (e.g., Smith et al. 2006; Borre et al. 2020) until being observed by TESS, which finally revealed the presence of groups (Nazé et al. 2020a; Labadie-Bartz et al. 2021). Nevertheless, it is clear that in many cases, isolated signals rise well above the noise floor and are apparently unrelated to groups.

Harmonics of isolated signals appear in about 6% (27/430) of cases (or, in 21% of stars with isolated signals; 27/126), again being more common toward later spectral types (where 6% (3/49), 26% (7/27), and 30% (12/40) of early-, mid-, and late-type stars with isolated signals have clear harmonics). These signals may reflect individual pulsation modes, or, at lower frequencies, may be related to rotation (either of the Be star, or perhaps a companion), a close binary (possibly in a hierarchical triple configuration with a relatively distant Be star), or are simply the strongest signal in an otherwise undetected group.

In the adopted characterization scheme, isolated signals must be roughly constant in amplitude over the observing baseline of TESS. Otherwise, they would appear as two or more signals in the frequency spectrum (not necessarily of equal amplitude) as demonstrated in Section 3.4. It is, however, certainly possible that seemingly isolated signals can vary in strength over longer time baselines. The case of γ Cas exemplifies this, where an apparently isolated signal at $0.82 \text{ day}^{-1}$ was known for many years, but with a decreasing amplitude to the point where it is no longer present in space photometry of the system (Henry & Smith 2012; Borre et al. 2020). The converse is also possible, where a seemingly isolated signal in TESS reflects a pair of frequencies with a beat period greatly exceeding the observational baseline.

## 5.6. High and Very High Frequencies

### 5.6.1. Main Results for High Frequencies ($6 < f < 15 \text{ day}^{-1}$)

A subset of the sample (15%, 63/430) has high-frequency signals that are not simply harmonics of the typical frequency groups (or isolated signals at lower frequencies), but are rather individual frequencies (or they can exist in groups, but with a different morphology or configuration than the typical lower-frequency groups). It should be kept in mind, however, that the short baseline of TESS implies a frequency resolution of $\approx 0.04 \text{ day}^{-1}$, so it is possible that what is here detected as single frequencies might be two or more unresolved frequencies. High-frequency signals are the only variability classification that differs when comparing the relatively uncontaminated subsample (last columns in Table 1) to the full sample. Therefore, it is possible that in some cases, high-frequency signals are due to blending and are not intrinsic in a given Be star.

High-frequency signals are seen in three of the 15 Be stars (20%) in Semaan et al. (2018), and in ~10% of the γ Cas analogs analyzed in Nazé et al. (2020a, 2020b), most notably in γ Cas and π Aqr. Similarly, ~10% of the Be stars in Balona & Ozuyar (2020) show these sort of signals. However, these fractions may be higher if the detection thresholds are relaxed. Some caution should be exercised since genuine higher-frequency signals may be (near) harmonically related to

lower-frequency groups. Nonetheless, the agreement seems good between the aforementioned and current studies, with high-frequency signals being detected in space photometry of ~10%–20% of Be stars.

β Cephei stars are early B-type stars (roughly B0–B2.5) that pulsate in p-modes with typical frequencies between 3 and $12 \text{ day}^{-1}$ (Stankov & Handler 2005). While Be stars most typically show low-frequency pulsation similar to the SPB stars (De Cat 2002), some classical Be stars are also observed to have β Cephei pulsation (e.g., Walker et al. 2005b; Huat et al. 2009; Nazé et al. 2020a; Labadie-Bartz et al. 2020). It is likely that some fraction of the Be stars that show these high-frequency signals in TESS also belongs to the class of Be stars that are SPB/β Cephei hybrid pulsators.

However, signals in the traditional β Cephei regime are not necessarily p-modes, because a pulsational signal in the SPB regime in the corotating frame can be pushed to higher observed frequencies in rapidly rotating stars via the equation $f_{\text{obs}} = f_{\text{co-rotating}} - m\Omega$, where $\Omega$ is the rotational frequency of the star, and $m$ is the azimuthal order of the pulsation mode (being negative for prograde modes, and positive for retrograde modes). However, such a situation usually appears in the frequency spectrum as a harmonic series of groups (e.g., panels (7) and (9) of Figure 5, with observed amplitudes usually decreasing with larger $|m|$), and these signals would therefore not be considered as (isolated) high frequencies. There is also the possibility of combination frequencies as described in Kurtz et al. (2015), which may be an important ingredient in the observed frequency spectra. We caution that any studies of Be stars that investigate high-frequency signals must also consider the mid- and lower-frequency signals. While a given system may have signals in the high-frequency regime, in many cases, these are related to a near-harmonic series of frequency groups (even if the lower frequencies are removed with iterative pre-whitening).

Perhaps surprisingly, the presence of high-frequency signals is not limited to early-type stars. This is counterintuitive if these signals reflect β Cephei pulsation, since the β Cephei phenomenon (i.e., p-modes excited by the $\kappa$ mechanism) is restricted to early B-type stars. We find high-frequency signals in 16% (34/218) of early-type, 17% (13/79) of mid-type, and 13% (14/110) of late-type stars. Of the 217 relatively uncontaminated systems (where <5% of the measured flux is from neighboring sources), 11.4% (13/114) of the early-type, 12.5% (5/40) of the mid-type, and 7.1% (4/56) of the late-type stars exhibit high-frequency variations.

### 5.6.2. Main Results for Very High Frequencies ($f > 15 \text{ day}^{-1}$)

Very high frequencies are uncommon, with only 3% (13/430) of this sample have signals beyond the typical β Cephei regime with frequencies greater than $15 \text{ day}^{-1}$. The highest-frequency signals detected are near $75 \text{ day}^{-1}$ in the system TIC 427400331 = HD 290662. Among the 13 Be stars with frequencies >$15 \text{ day}^{-1}$, only four have signals above $24 \text{ day}^{-1}$ (see Appendix B.1). Signals at the lower end of this very-high-frequency regime may be typical p-mode pulsation shifted to higher observed frequencies due to rapid rotation, as is suggested for the signals at 17.27 and $19.31 \text{ day}^{-1}$ in the Oe star ζ Ophiuchi (Walker et al. 2005b). Or, in some cases, combination frequencies can produce signals at very high frequencies, and therefore some very-high-frequency signals may not be independent, but rather a sum of two lower-frequency signals (or perhaps a more complex combination; Kurtz et al. 2015; Burssens et al. 2020). To search for combination frequencies in
these stars, a few of the highest-amplitude signals were selected to compute potential combination frequencies of the form $nf_i + mf_j$ with orders $n, m \in [-3, 3]$ (as done in Kurtz et al. 2015; Burssens et al. 2020). Indeed, in some (but not all) of the systems with very high frequencies, many or all of these signals are combinations of signals with $f < 15$ day$^{-1}$. For example, in TIC 140001327 (see Section B.1), the very-high-frequency signals at $16.586$ day$^{-1}$, $18.087$ day$^{-1}$, $25.568$ day$^{-1}$, and $30.871$ day$^{-1}$ are linear combinations of lower-frequency signals (namely at 1.501, 3.093, 5.303, 7.481, 9.831, and 11.283 day$^{-1}$). A thorough analysis of this sort is beyond the scope of this paper, but has the potential to add valuable diagnostic information since the presence of combination frequencies implies that the parent modes oscillate in the same star and not from separate unresolved companions.

However, toward the higher end of these very high frequencies ($\gtrsim 30$ day$^{-1}$), and especially for Be stars of later spectral type, it is unlikely that they arise in the Be star. Instead, these may indicate the presence of a stellar companion pulsating at these relatively high frequencies. This may be supported by the higher incidence of very-high-frequency signals in later spectral types (2% (5/218), 3% (2/79), and 6% (6/110) in early-, mid-, and late-types, respectively), which is an expected trend if these signals do indeed indicate a companion since the contrast ratio between the Be star and a hypothetical companion will be less severe for Be stars of later spectral type. However, the number of systems with very high frequencies is far too small to draw meaningful conclusions. Another possibility is that these modes do originate in the Be star despite being outside the theoretical instability strips where modes of these frequencies are expected.

Other studies seem to agree with the rarity of very-high-frequency signals. Only one star (HD 45314 = PZ Gem = TIC 438306275) of 15 studied in Nazé et al. (2020b) has a very-high-frequency signal ($f = 18.080$ day$^{-1}$). However, this is a marginal detection, and although this star is included in the present study, this signal is not considered significant (but may indeed be real). Two of the Be stars in the 57 star sample of Balona & Ozuyar (2020), TIC 139385056 (HD 58978 = FY CMa) and TIC 148316007 (HD 49319 = HR 2507), are listed as having very high frequencies up to $\sim 20$ day$^{-1}$. FY CMa is also analyzed in the present study, but again the pair of signals near 18 and $22$ day$^{-1}$ reported in Balona & Ozuyar (2020) are of very low significance and were not included here. Likely due to both the low amplitudes and rarity of frequencies $> 15$ day$^{-1}$, these sort of signals are not discussed much in the literature for Be stars and are typically thought to reflect relatively high-frequency p-mode pulsation.

All 13 Be star systems (plus three systems rejected from the sample) with very high frequencies are briefly discussed in Appendix B.1, including plots of their frequency spectra. These should be studied further to confirm or reject their classical Be nature, and to attempt to resolve the origin of these very-high-frequency signals. For stars with 2 minutes cadence light curves, frequencies out to $360$ day$^{-1}$ can be probed, as opposed to $24$ day$^{-1}$ for 30 minutes cadence data. This limits our ability to detect the highest-frequency signals to the $\sim 65\%$ of systems with 2 minutes cadence data, and so the percentage of systems in which very high frequencies are detected is likely underestimated.

### 5.6.3. Binarity and Composite Frequency Spectra

Binarity is common in the massive star population, with the majority of B-type stars being members of a multiple star

system (e.g., Kouwenhoven et al. 2007). Oudmaijer & Parr (2010) found essentially the same binary fraction and properties when comparing Be versus B stars. While the binary fraction of Be stars is poorly constrained (due largely to their broad absorption lines and high luminosity), many Be stars are known to have a binary companion (often sdO stars; Wang et al. 2021). The frequency spectra of binary systems will then contain signals from two (or more) sources if they are variable. In fact, the prototype of the $\beta$ Cephei class, $\beta$ Cephei, is a binary with a slowly rotating $\beta$ Cephei primary and a rapidly rotating mid-B-type classical Be star (Wheelwright et al. 2009). Often in Be binaries, the companion to the Be star is an evolved object, such as a hot subdwarf (sdOB) star (Gies et al. 1998; Chojnowski et al. 2018; Wang et al. 2018, and references therein), a neutron star (i.e., a Be X-ray binary system; Ziolkowski 2002), or even a black hole (Casares et al. 2014). Despite the fact that a Be companion does not need to be significantly evolved, a recent study (Bodensteiner et al. 2020a) found no firm evidence for main-sequence companions around early Be stars (note, however, that low-mass main-sequence companions are difficult to detect because the total flux of the system tends to be dominated by the Be star, and that this study does not include the aforementioned $\beta$ Cephei system where the Be star is of a later spectral type and is distant from the B0.5 primary). Some unknown fraction of this TESS sample of Be stars is in binaries. It is therefore possible that some of the high- and very-high-frequency signals detected in the TESS data arise not from the Be star, but from a companion.

### 5.6.4. Potential Pulsating Companions

Some examples of stars with pulsation at greater than $15$ day$^{-1}$ include the main-sequence A-type roAp stars with periods between 6 and 20 minutes ($f \sim 72-240$ day$^{-1}$) and amplitudes of a few millimagnitude (Kurtz 1982), and the (usually, but not always, main-sequence A/F type) $\delta$ Scuti stars with typical periods between 0.02 and 0.25 day ($f \sim 4-50$ day$^{-1}$) and typical amplitudes between 3 and 10 mmag, but sometimes up to and in excess of 0.3 mag (the high-amplitude $\delta$ scuti stars; Breger 2000; Garg et al. 2010).

Perhaps most relevant for Be star systems are those with sdOB companions, where at some earlier time the initially more massive star expanded and donated mass and angular momentum to the present-day Be star, shedding its envelope and leaving behind its core (the sdOB star). Such Be+sdOB systems typically have (near) circular orbits with periods of roughly 100 days (e.g., Bozic et al. 1995; Bjorkman et al. 2002; Peters et al. 2016; Chojnowski et al. 2018). sdOB stars are also known to pulsate in this high-to-very-high-frequency regime. In sdOB stars, p-mode pulsations typically have periods between 2 and 10 minutes ($f \sim 144-720$ day$^{-1}$), and amplitudes of about 1% ($\sim 10$ mmag; Kilkenny et al. 1997), and g-modes oscillate with periods between about 45 minutes and 2 hr ($f \sim 12-32$ day$^{-1}$) and amplitudes typically between 0.1 and 0.5% ($\sim 1-5$ mmag; Green et al. 2003; Kawaler et al. 2010). However, in some cases, the amplitudes can be much larger, and an sdOB star can pulsate in both p- and g-modes. For example, the sdO star PG 1605+072 shows over 55 modes, the strongest of which has a period of 8.03 minutes ($f = 180$ day$^{-1}$) and an amplitude of over 50 mmag (about 5%, Pereira & Lopes 2004). Sahoo et al. (2020) analyzed TESS data for three hot subdwarfs, finding frequencies between about 9
and 300 day⁻¹, with most of the signals being between about 15 and 30 day⁻¹.

In most cases, the Be star dominates the visible flux in Be binary systems. This can severely dilute any signals originating in a companion, to the point of being undetectable even with space photometry. However, there are many possible configurations where such a signal in a companion can rise above the detection threshold despite the contaminating flux from the Be star. For example, in the well-known B0.5Ve+sdO binary φ Per (Poeckert 1981; Gies et al. 1998), the sdO star contributes approximately 3% of the total visible flux (Mourard et al. 2015). If the sdO star were to pulsate with an amplitude of 50 mmag (like PG 1605+072), this signal would be diluted down to an amplitude of ~1.5 mmag, which is easily detectable in space photometry provided the observing cadence is such that the frequency of the signal is not beyond the Nyquist limit, and the exposure time is sufficiently short. The flux ratio between the Be star and its companion can be even lower for later-type Be primaries, such as the B5e primary in the 7 Vul binary system (which is suspected to have an sdOB companion; Vennes et al. 2011). However, it should be noted that most confirmed Be+sdOB binaries have primaries with spectral types between B0 and B3 (Wang et al. 2018). Two minutes cadence TESS light curves are much better suited to this task (with a Nyquist limit of 360 day⁻¹) compared to the 30 minute light curves (Nyquist limit of 24 day⁻¹).

These systems with very high frequencies are therefore good candidates for further study to search for evidence of binarity. This can be done through radial velocity monitoring of the Be star and/or its disk (Bjorkman et al. 2002; Miroshnichenko et al. 2002), direct detection of spectral features (and their radial velocity) of a companion (Bjorkman et al. 2002; Chojnowski et al. 2018; Wang et al. 2018), modeling of the radial structure of the Be star disk to infer truncation from a binary companion (Klement et al. 2019), and to observe the radial and azimuthal structure of the disk to detect density waves that propagate at the orbital period, which are caused by tidal forces from a binary companion (Panoglou et al. 2016; Chojnowski et al. 2018; Panoglou et al. 2018, 2019; Cyr et al. 2020). Direct interferometric detection may also be possible.

This is an important topic because the binary fraction of Be stars is poorly constrained, and one of the proposed evolutionary channels by which a Be star achieves its near-critical rotation is through binary interaction (Pols et al. 1991; de Mink et al. 2013), as opposed to, or in addition to, outward angular momentum transfer from a contracting core in a single stellar evolution scenario (Ekström et al. 2008; Granada et al. 2013).

## 5.7. No Variability Detected

There are 13 Be star systems where no variability is detected, spanning magnitudes between $T_{mag} = 6.4–10.3$ (see Table 2). Eleven (two) of these systems are of late (mid) spectral type. This does not necessarily mean that these systems do not have variability, as signals below the detection threshold remain possible. Very low photometric amplitudes attributed to pulsation are often found in late-type Be stars. For example, a few signals with photometric amplitudes less than 0.05 mmag exist in Kepler data for the Be star ALS10705 (Rivinius et al. 2016). Further, in three of these cases (TIC 143543729, 206840215, and 247589847), there is enough spectroscopic data to confirm the presence of a disk that is variable over time (e.g., in TIC 247589847 = BD+13 976 = ABE-083, APOGEE

spectra show a double-peaked, slightly variable Br11 emission feature over a 393 day baseline; Chojnowski et al. 2015). That is, at least these three systems are in fact variable (at least in regards to their mass ejection rate) despite the nondetections in TESS.

# 6. Discussion of Overall Results

Virtually all (97%) of the Be stars observed by TESS in sectors 1–13 are variable. Their photometric signals carry information about the underlying physical processes causing this variability, and a careful study of these observations can elucidate the nature of Be stars as a population especially in regards to pulsation and mass ejection episodes. While there remains much to be explored with these data, the overview of Be star variability seen by TESS presented in this work provides a solid background of their photometric behavior as a population on short timescales and is a first step toward fully leveraging the unique data set provided by TESS for these and similar objects. A discussion for each characteristic variability pattern is given in Section 5. The remainder of this section highlights correlations between different classes of variability.

## 6.1. Correlations between Variability Classifications

Figure 11 shows correlations between the observed characteristics of the sample (including spectral type). One of the most interesting correlations is between flickers and frequency groups. If there were no correlation, then 91% of the early-type stars with flickers are expected to also have frequency groups (61 out of 67), yet every star with flickers exhibits frequency groups (and, in 83% of systems with flickers, there is an increase in the amplitude of groups coinciding with the flicker event). This is suggestive of a physical link between frequency groups and the short-duration mass-loss events that flickers are interpreted to represent.

This is not to say that frequency groups are necessary for a Be star to build up or sustain a disk. For example, HR 6819 (TIC 118842700) does not contain frequency groups, but rather is dominated by low frequency and stochastic variability and unambiguously supports a disk (Rivinius et al. 2020).<fnref n="11" /> A caveat is that, given the transient nature in many cases of Be star variability patterns, it remains possible that groups sometimes exist in HR 6819 outside of the short observing baseline provided by TESS. Similar examples include HD 84567 (TIC 11972111; Shokry et al. 2018), HD 53667 (TIC 177204351), HD 254647 (TIC 291385725), HD 44637 (TIC 438103655), and omi Pup (TIC 127493611; Koubský et al. 2012).

In stars with flickers that have the typical group configuration (where there are two or more groups that stand out above the noise, and the higher-frequency groups have, very roughly, nearly twice the frequency of the first), there is a difference in the relative strength of $g1$ and $g2$ when compared with the full sample (as shown in Figure 7). In the 71 stars with flickers and typical groups, 12 have $g1$ stronger than $g2$ (17%), 31 (44%) have $g1$ and $g2$ having similar strength, and 28 (39%) have $g2$ stronger than $g1$. In other words, in stars with flickers, there is a tendency for $g2$ to be relatively stronger compared to the relative group strengths of the full sample. Although not quantified at present, there is a qualitative trend of $g2$ seeing a

<footnote n="11">¹¹ HR 6819 is a multiple system that includes a narrow-lined B3 III star, which contributes about 50% of the total visible flux.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Correlations between various characteristics of the sample (blue (red) = positive (negative) correlation). The size of each dot is proportional to the strength of the (anti)correlation. From top to bottom, “early ST,” “mid ST,” “late ST,” and “unknown ST” refer to the spectral type bins. “Groups” include stars with one or more frequency groups, and “typical groups” refer to those with the typical group configuration. “Low freq,” “flickers,” “stochastic,” and “isolated freqs.” refer to the variability classifications introduced in Section 4. “High freqs.” and “very high freqs.” refer to systems with high ($6<f<15\ \text{day}^{-1}$) and very high ($f>15\ \text{day}^{-1}$) frequencies. The remaining categories describe the relative strength of $g1$ and $g2$ in systems with the typical group configuration, where $g1$ is stronger (“$A_{g1}>A_{g2}$”), where both groups are of similar strength (“$A_{g1}=A_{g2}$”), and where $g1$ is weaker than $g2$ (“$A_{g1}<A_{g2}$”).</figcaption></figure>

more dramatic enhancement than $g1$ during flicker events (e.g., Figures 4, 8), but there are exceptions to this.

In systems where low-frequency stochastic variability is a prominent feature, there is a negative correlation with frequency groups, and a positive correlation with high-frequency signals, as well as positive correlations with flickers and being dominated by low frequencies. These (anti) correlations are apparent in the examples shown in Figure 10.

There is an anticorrelation between late-type stars and most aspects of the variability of the sample, except for isolated frequencies (including high and very high frequencies). This is perhaps due to photometric amplitudes in late-type stars being relatively low (in regards to pulsation, low-frequency stochastic excess, and any disk variability due to intrinsically low density), which may be further exacerbated by any stellar companions. Or, signals may be very closely spaced such that they are not resolved in systems observed in only one or a few TESS sectors (e.g., TIC 279430029 = HD 53048, B6Vne,

where the two strongest signals in $g1$ are separated by only $0.0022\ \text{day}^{-1}$ but are resolved because this system was observed in all 13 sectors in cycle 1).

## 7. Conclusions

Following in the footsteps of decades of observations from the ground, and in more recent years also from space, the TESS mission now, for the first time, provides high-precision space photometry for the majority of known Galactic Be stars, allowing for a systematic study of their variability on timescales from minutes up to tens of days and down to amplitudes of approximately 50 ppm. Analysis of TESS data for 430 classical Be stars observed in its first year of operation confirms that virtually all Be stars are variable (97% of this sample is variable at the level of precision available with TESS). The stars in this sample show a variety of characteristic signals, the rates of incidence of which are summarized in
Table 1. Understanding the cause of these features and their incidence rates and patterns is an important step toward elucidating the physical processes in Be stars and other rapid rotators.

In general, Be stars show a higher level of pulsational variability than nonrapidly rotating stars of the same spectral type (Diago et al. 2009), and there is mounting evidence that rapid rotation enhances pulsational amplitudes and the number of excited frequencies (Rieutord 2009; Neiner et al. 2012a, 2020). It is becoming increasingly clear that pulsation is a common, and likely ubiquitous, element of Be stars (Rivinius et al. 2013; Semaan et al. 2018), which is supported by this study. In particular, the commonality of frequency groups may have important consequences in regards to the Be phenomenon, as this is the most common feature of the sample studied in this work. While it is possible that rotation contributes in some degree to the observed variability of the sample (via an inhomogeneous stellar surface and/or corotating clouds), rotation alone is insufficient to explain the majority of the observed signals. Rather, nonradial pulsation (NRP) has properties that can explain much of the observed variability in this sample. NRP in rapid rotators naturally forms frequency groups (Saio 2013; Kurtz et al. 2015; Saio et al. 2018a), NRP-modes can transport angular momentum to the surface layers (Bowman et al. 2019b, 2020), which can assist in the triggering of outbursts (Neiner et al. 2020), and NRP-modes can couple in the resonant cavity in which they oscillate to produce combination frequencies with amplitudes higher than the sum of the parent modes (which can explain the lowest-frequency group, g0, in many Be stars, and higher-order groups; Kurtz et al. 2015; Baade et al. 2018a). It is clear that high-frequency signals (taken to be $f > 6$ day$^{-1}$, but excluding harmonics of lower-frequency signals) are related to pulsation, and not rotation (as such a rotational frequency would exceed the critical rotation rate). Further, stochastically driven internal gravity waves (a type of nonperiodic pulsation) have emerged as a plausible explanation for the low-frequency excess and stochastic variability observed in many OB stars (Bowman et al. 2019b, 2020).

Conversely, it is also important to study the Be stars that do not exhibit groups, but instead often are dominated by low-frequency stochastic variability and/or one or more isolated signals. Understanding the physical properties of the stars belonging to these two categories (with and without groups), especially their rotation rate and evolutionary status, may lead to an improved understanding of how different Be stars eject mass. For example, a hypothesis that can be tested with spectroscopic or polarimetric time-series data is that systems with a high degree of low-frequency stochastic variation, but without groups, may feed a disk in a more continuous fashion, compared to the episodic mass-loss episodes commonly seen in systems with well-defined groups. Semaan et al. (2018) found that both of the Be stars that are evolved slightly past the main sequence, out of their sample of 15, show only low-frequency stochastic variability, in contrast to the majority of the sample, which exhibits frequency groups (12 out of the 13 main-sequence stars). Constraining the evolutionary status of the systems observed by TESS may therefore be of great importance in interpreting the observed variability patterns.

A nonnegligible fraction of the Be stars in this sample (18% of the total, and 31% for early-type stars) exhibit flickers, which are interpreted as short-lived episodes of mass ejection

(see Figure 9 for a proof of concept, and Figures 4 and 8 for additional examples). TESS is uniquely suited for capturing such events, which can last for only a few days and change the net brightness of the system by a few percent. Similar events have been observed from space in many Be stars prior to TESS (Rivinius et al. 2013; Baade et al. 2016, 2017; Semaan et al. 2018), but the large number of Be stars observed with TESS finally allows for a substantial sample to be studied. Future work will aim to better quantify the occurrence rate, shape, timing, amplitude, duty cycle, dependence on inclination angle, and dependence on the strength of any preexisting disk of these events. Hydrodynamic and radiative transfer codes will allow for the modeling of these events, providing estimates of the flux of mass and angular momentum from the star, which will inform estimates of the degree to which these relatively small-scale mass ejection events contribute to the total mass budget of the disk.

It is well known that mass-loss episodes in Be stars can last for decades, years, months, weeks, or days (Labadie-Bartz et al. 2017, 2018; Rímulo et al. 2018). Whatever mechanisms are responsible for opening the mass-loss valve in Be stars must be able to account for this large range in timescales. The increased amplitude of frequency groups coinciding with flicker/outburst events that are commonly seen in this sample (frequency group enhancement coinciding with flickers is seen in 83% of systems displaying these events) is suggestive of pulsation playing a key role, at least in the relatively short-lived events captured by TESS. However, we caution that in some situations, the enhancement of frequency groups can instead be a consequence of mass ejection if the variability is in some part circumstellar (or perhaps a temporary increase in the amplitude of r-modes). TESS, with its short observational baseline, is not sensitive to variability with timescales of tens of days and longer.

The contribution to massive star science with TESS has only just begun, yet it is already producing important results for large samples (e.g., David-Uraz et al. 2019; Balona & Ozuyar 2020; Bowman et al. 2020; Burssens et al. 2020). Where other space photometry missions have paved the way toward studying Be stars at high photometric precision, TESS is continuing these opportunities for by far the largest sample yet of bright OB stars. The results presented in this paper are intended to be an overview of the types of variability that are seen in Be stars with TESS, and more in-depth results regarding the nature of these signals and their significance will be explored further in forthcoming works by our group, and undoubtedly others. Especially exciting, given the brightness of TESS targets, is the relative ease of monitoring these stars with other observational techniques (spectroscopy, polarimetry, multiband photometry, and interferometry) from the ground, which, when combined with TESS data, can greatly improve our understanding of these objects.

J.L.-B. acknowledges support from FAPESP (grant 2017/23731-1). A.C.C. acknowledges support from CNPq (grant 311446/2019-1) and FAPESP (grants 2018/04055-8 and 2019/13354-1). T.A. acknowledges support from FAPESP (grant 2018/26380-8). A.R. acknowledges support from CAPES (grant 88887.464563/2019-00). A.L. acknowledges support from CAPES (grant 88882.332925/2019-01). P.S. acknowledges support from FAPESP (grant 2020/04445-0). The authors are grateful to the amateur spectroscopy
community, whose observations have directly supported this and forthcoming works, in particular Stephane Charbonnel, Christian Buil, Paolo Berardi, Tim Lester, Alain Maetz, Robin Leadbeater, Erik Bryssinck, Olivier Garde, Franck Houpert, and Olivier Thizy whose data was used in this publication. We thank Martin Paegert for assisting in the calculation of the contamination ratios. We thank the anonymous referee for numerous suggestions that greatly improved the clarity and presentation of this work. This work makes use of observations from the LCOGT network. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes. Funding for the TESS mission is provided by NASA's Science Mission directorate. TESS Guest Investigator program G011204 provided the 2 minutes cadence data. This project makes use of data from the KELT survey, including support from The Ohio State University, Vanderbilt University, and Lehigh University, along with the KELT follow-up collaboration. 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. This research has made use of NASA's Astrophysics Data System. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This work has made use of the BeSS database, operated at LESIA, Observatoire de Meudon, France: http://basebe.obspm.fr This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration et al. 2018). This research made use of Astropy,<fnref n="12" /> a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018). This work used the BeSOS Catalogue, operated by the Instituto de Física y Astronomía, Universidad de Valparaíso, Chile: http://besos.ifa.uv.cl and funded by Fondecyt iniciación No. 11130702, using the PUCHEROS instrument (Vanzi et al. 2012). The page is maintained thanks to FONDECYT No. 11190945. This research has made use of the services of the ESO Science Archive Facility.

Facilities: TESS, Gaia, KELT.

Software: astropy (Astropy Collaboration et al. 2013, 2018), Period04 (Lenz & Breger 2005), VARTOOLS (Hartman 2012).

## Appendix

Appendix A describes the clustering algorithm used to characterize frequency groups. The remaining sections provide remarks on individual systems, beginning with those having very high frequencies (Appendix B.1). Systems confirmed to not be classical Be stars are then briefly discussed in Appendix B.2, followed by systems strongly suspected to not be Be stars (and are thus rejected from the sample; Appendix B.3), and systems weakly suspected of not being classical Be stars, but kept in the sample due lacking solid grounds for rejection (Appendix B.4). Appendix C includes a table (Table 2) for all of the stars in the sample and the variability characteristics assigned to them. Appendix D

describes the plots made for each star in the sample and shows an example of such a figure.

## Appendix A
### Frequency Group Clustering Algorithm

In order to determine groups of frequencies and their “center of mass,” a python routine was developed to automatically identify these properties from the discrete iteratively pre-whitened frequency spectrum for each star. This clustering algorithm served to identify groups and provide a number to describe their weighted center and net amplitude.

Broadly speaking, this routine starts with identifying the signals with the highest amplitude (the “seed” signal for that group), and then moves outwards from these frequencies to identify signals that belong to the same group, while retaining information about the strength of each signal and ultimately computing the group center (with weights given in proportion to the amplitude of each signal belonging to a given group). As the algorithm moves outwards from the highest-amplitude signals, each frequency is compared to the seed signal of the group by the equation

$$A(f_a) + A(f_b) \geq \frac{|f_a - f_b|^2}{d}, \quad (1)$$

where $A(f)$ is the amplitude of the frequency $f$, $f_a$ and $f_b$ are the frequencies being compared, and $d$ is a weight parameter, which is calculated by

$$d = \frac{\alpha \times N_f \times I}{A_{\text{total}}}, \quad (2)$$

where $\alpha$ is a free parameter (described in the following paragraph, and with units of day⁻¹), $N_f$ is the total number of signals in the frequency spectrum out to 24 day⁻¹, $I$ is the interval of frequencies considered ($I = 24 \text{ day}^{-1}$ in our analysis), and $A_{\text{total}}$ is the sum of the amplitudes of all $N_f$ frequencies. If Equation (1) is false for the given frequencies, the frequency is not included in the group. Otherwise, the compared frequencies are merged inside the group by replacing both with a new signal having a frequency of

$$f_{\text{new}} = \frac{f_a \times A(f_a) + f_b \times A(f_b)}{A(f_a) + A(f_b)}$$

and an amplitude of

$$A(f_{\text{new}}) = A(f_a) + A(f_b),$$

which is subsequently compared to the remaining signals in the spectrum. After all frequencies are analyzed in this way, the routine returns the resulting clusters of frequencies centered at the weighted average for the group, and with a characteristic strength computed from the sum of the amplitudes of the original signals determined to belong to the group.

Because of the diversity of frequency spectra in the full sample, there is no single value for the $\alpha$ parameter that returns reliable results for all stars. Instead, this algorithm was applied to each star with three empirically determined values of $\alpha$ ($5 \times 10^{-4}$, $9 \times 10^{-4}$, and $13 \times 10^{-4}$), where the authors then choose from these three outputs the results that best reflect the behavior of a given frequency spectrum (considering both the 30 minutes and 2 minutes data whenever available). In practice, often all three values of $\alpha$ return virtually the same group

<footnote n="12">¹² http://www.astropy.org</footnote>
information, but there are cases where a certain value of $\alpha$ is clearly preferable (e.g., in densely populated frequency spectra when groups are close together, a high value of $\alpha$ may erroneously merge two adjacent groups). A consensus among the authors as to the existence, location, and relative strength of $g1$ and $g2$ for each system was finally arrived at after considering the output from the clustering algorithm.

Although this method is adequate as a first step toward quantifying frequency groups, alternative methods may provide more robust results. For example, a simultaneous fit to the frequency spectrum of a red noise profile plus periodic signals seems better suited to the classical Be stars. Then, the red noise profile can be multiplied by some (small) constant and used as a threshold so that the individual frequencies that make up a group can then be determined and used in calculating the group width, center, and overall amplitude.

## Appendix B
### Notes on Individual Systems

#### B.1. Stars With Very High Frequencies ($f > 15 \text{ day}^{-1}$)

All 13 of the Be star systems (plus three rejected from the sample) with very high frequencies are briefly discussed below. Their frequency spectra are shown in Figure 12. All of these were tested to see if the very high frequencies are combinations of lower-frequency signals. Unless otherwise stated, the listed signals in these stars do not appear to be combination frequencies.

**TIC 11411724 = StHA 52:** B1.5V. Likely not a classical Be star (see Section B.3). Embedded in strong reflection nebula. There is a weak pair of frequencies around $20 \text{ day}^{-1}$ ($20.072 \text{ day}^{-1}$ and $20.609 \text{ day}^{-1}$, both at $\sim 0.3$ ppt), and many signals between 8 and $15 \text{ day}^{-1}$, the strongest of which include $8.465 \text{ day}^{-1}$ (0.7 ppt), $8.635 \text{ day}^{-1}$ (0.6 ppt), $8.791 \text{ day}^{-1}$ (0.7 ppt), $8.924 \text{ day}^{-1}$ (0.5 ppt), $9.678 \text{ day}^{-1}$ (0.5 ppt), $9.948 \text{ day}^{-1}$ (0.9 ppt), $10.325 \text{ day}^{-1}$ (0.6 ppt), $11.161 \text{ day}^{-1}$ (1.1 ppt), $11.423 \text{ day}^{-1}$ (0.4 ppt), $13.775 \text{ day}^{-1}$ (0.4 ppt), and $14.156 \text{ day}^{-1}$ (0.4 ppt). These values are from 2 minutes cadence data.

**TIC 75581184 = HD 80156:** B8.5IVe. The most prominent feature is a group centered at $2.70 \text{ day}^{-1}$. Harmonics of this group seemingly extend out to even $50 \text{ day}^{-1}$. It is most likely that these very-high-frequency signals represent harmonics of lower-frequency signals, but this system is unique in the extent of the apparent harmonic series. A possible explanation is that the dominant low-frequency signals are extremely noninusoidal; however, in a wavelet plot, the strengths of these groups are not apparently correlated (nor are they anticorrelated). Instead, these may represent a series of groups of increasing $|m|$ (out to $|m| \gtrsim 8$). A notable feature is that the groups near odd-numbered harmonics of the first group are relatively strong. That is, the group at $\sim 4 \times g1$ (i.e., the third harmonic) is stronger than the group at $\sim 3 \times g1$, and similarly for $\sim 6 \times g1$ versus $\sim 5 \times g1$, and so on, similar to what is seen in panels (7), (8), and (9) of Figure 5. This analysis is based on 2 minutes cadence data.

**TIC 80719034 = HD 67985:** B8Vne. Two typical, but very narrow, groups centered at $2.99$ (the strongest signal has amplitude of 2 ppt) and $5.96$ (the strongest signal has amplitude of $6.5$ ppt) day$^{-1}$, and a low-frequency signal near $0.36 \text{ day}^{-1}$ (amplitude $\sim 6$ ppt; and possibly its first harmonic). There are

many signals between 11 and $22 \text{ day}^{-1}$ with amplitudes between 0.5 and 2 ppt. The strongest of these, at $18.510 \text{ day}^{-1}$ (2.1 ppt), is the third strongest signal overall. Other prominent signals include $11.327 \text{ day}^{-1}$ (1.7 ppt), $12.040 \text{ day}^{-1}$ (1.4 ppt), $12.423 \text{ day}^{-1}$ (0.5 ppt), $12.952 \text{ day}^{-1}$ (0.6 ppt), $13.284 \text{ day}^{-1}$ (0.7 ppt), $17.512 \text{ day}^{-1}$ (0.6 ppt), $18.042 \text{ day}^{-1}$ (0.6 ppt), $18.501 \text{ day}^{-1}$ (0.6 ppt), $20.403 \text{ day}^{-1}$ (0.6 ppt), $21.476 \text{ day}^{-1}$ (0.7 ppt), and $22.300 \text{ day}^{-1}$ ($=2 \times f_{11.327} - f_{0.354}$, 0.4 ppt). These values are from 2 minutes cadence data.

**TIC 123828144 = HD 62894:** B8e. There is a triplet centered around $15.6 \text{ day}^{-1}$, made up of the signals $15.202 \text{ day}^{-1}$ (0.12 ppt), $15.641 \text{ day}^{-1}$ (0.21 ppt), and $16.165 \text{ day}^{-1}$ (0.05 ppt). There is an isolated signal at $6.222 \text{ day}^{-1}$ (0.25 ppt). The dominant signals in the LC are a relatively strong (2 ppt) signal at $0.96 \text{ day}^{-1}$, and its subharmonic (at $0.48 \text{ day}^{-1}$, 0.4 ppt), plus some lower-level stochastic variability. Blending is a concern, as there is a neighbor of equal brightness (CD-42 3473) at a distance of $24''$, which is completely blended in the TESS pixels (contamination ratio = 1.036). This object is included in Labadie-Bartz et al. (2020) as a potential $\beta$ Cephei pulsator, where in KELT data (Labadie-Bartz et al. 2017), the $0.96 \text{ day}^{-1}$ signal is recovered, with the only other detected frequency being at $5.52537 \text{ day}^{-1}$, and it is noted that “There is an A-type star of equal brightness $24''$ distant from this object, which could be a $\delta$ Scuti pulsator.” It should first be determined if these very-high-frequency signals arise from the neighbor (which seems likely considering the spectral types). Only 30 minutes cadence is available for this target.

**TIC 140001327 = HD 72014:** B1.5Vnne. This frequency spectrum is rich and includes many combination frequencies. There are many signals between 5 and $19 \text{ day}^{-1}$ including at $5.303 \text{ day}^{-1}$ ($=2 \times f_{0.401} + 3 \times f_{1.501}$, 0.47 ppt), $7.481 \text{ day}^{-1}$ (0.20 ppt), $8.239 \text{ day}^{-1}$ (0.21 ppt), $8.354 \text{ day}^{-1}$ (0.19 ppt), $9.121 \text{ day}^{-1}$ (0.49 ppt), $9.831 \text{ day}^{-1}$ (0.30 ppt), $11.283 \text{ day}^{-1}$ (0.35 ppt), $11.331 \text{ day}^{-1}$ ($=f_{1.501} + f_{0.831}$ or $=f_{3.093} + f_{8.239}$, 0.31 ppt), $12.784 \text{ day}^{-1}$ ($=f_{1.501} + f_{11.283}$ or $=f_{5.303} + f_{7.481}$, 0.29 pp), $14.285 \text{ day}^{-1}$ ($=2 \times f_{1.501} + f_{11.283}$, 0.14 ppt), $16.586 \text{ day}^{-1}$ ($=f_{5.303} + f_{11.283}$, 0.28 ppt), and $18.087 \text{ day}^{-1}$ ($=2 \times f_{5.303} + f_{7.481}$, 0.23 ppt). Two signals are tentative at $25.568 \text{ day}^{-1}$ ($=2 \times f_{5.303} + 2 \times f_{7.481}$, 0.04 ppt) and $30.873 \text{ day}^{-1}$ ($=3 \times f_{5.303} + 2 \times f_{7.481}$, 0.04 ppt). There are lower-frequency groups, but with a somewhat atypical pattern being centered at $1.52 \text{ day}^{-1}$, $3.09 \text{ day}^{-1}$, $3.61 \text{ day}^{-1}$, $4.52 \text{ day}^{-1}$, and $5.31 \text{ day}^{-1}$. Additionally, there are short flickers and low-frequency stochastic excess. This analysis is based primarily on 2 minutes cadence data.

**TIC 144028101 = mu Lup:** B8Ve. Very-low-amplitude signals in TESS at frequencies of $1.53$ and $2.44 \text{ day}^{-1}$ (amp = 0.07, 0.06 ppt, respectively), and also two pairs of signals at $14.515 \text{ day}^{-1}$ (0.02 ppt) and $14.845 \text{ day}^{-1}$ (0.11 ppt), and $17.673 \text{ day}^{-1}$ (0.05 ppt) and $17.989 \text{ day}^{-1}$ (0.06 ppt). There is a nice symmetric double-peaked $H\alpha$ profile in BeSOS (http://besos.ifa.uv.cl) in two epochs of 2015 (February and July). The high- and very-high-frequency signals are convincing despite their low amplitude. mu Lup itself is a close visual double (separated by about $1''$) where both components are approximately the same brightness and spectral type (Fabricius et al. 2002). Further, an early A-type neighboring star, HD 135748 (approximately 2 mag fainter), lies $23''$ away. These three sources share a common proper motion (Gaia
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Frequency spectra of systems with very high frequencies. Vertical dotted lines are at 6, 12, and $15\ \text{day}^{-1}$ for reference. In the top panel, the x-axis extends to 100 day$^{-1}$, while the remainder share the same axis, out to $40\ \text{day}^{-1}$. The horizontal red line is $4 \times$ the mean amplitude for $f > 15\ \text{day}^{-1}$, shown only for reference. TIC IDs marked with an asterisk are rejected as classical Be stars (see Sections B.2 and B.3).</figcaption></figure>

Collaboration 2018). With TESS data alone, it is impossible to distinguish the origin of the various photometric signals. 30 minutes cadence data were used for this analysis.

TIC 167110617 = HD 43285: B6Ve. Very unusual frequency spectrum. There are three prominent pairs or triplets roughly evenly spaced at 5, 7.5, and $10\ \text{day}^{-1}$, made up of
signals at $4.859\ \text{day}^{-1}$ (0.07 ppt), $4.975\ \text{day}^{-1}$ (0.15 ppt), $5.114\ \text{day}^{-1}$ (0.06 ppt), and $7.446\ \text{day}^{-1}$ (0.10 ppt), $7.602\ \text{day}^{-1}$ ($=2 \times f_{4.859} - f_{2.117}$, 0.21 ppt), and $9.714\ \text{day}^{-1}$ (0.05 ppt), $9.896\ \text{day}^{-1}$ (0.20 ppt). There are weaker signals between 12 and $24\ \text{day}^{-1}$, including at $11.974\ \text{day}^{-1}$ ($=2 \times f_{5.117} + f_{1.741}$, 0.04 ppt), $14.784\ \text{day}^{-1}$ (0.06 ppt), $19.582\ \text{day}^{-1}$ ($=3 \times f_{5.117} + 2 \times f_{2.117}$, 0.02 ppt), and $23.048\ \text{day}^{-1}$ (0.03 ppt), along with clear low-frequency signals consistent with SPB pulsation at $1.741\ \text{day}^{-1}$ (0.1 ppt), $2.117\ \text{day}^{-1}$ (0.06 ppt), and $3.148\ \text{day}^{-1}$ (0.03 ppt). Since there are both high (presumably p-mode pulsation) and low (presumably g- and/or r-mode pulsation) frequencies and combinations between them (e.g., at 11.974 and $19.582\ \text{day}^{-1}$), it seems as though all of these signals arise in the same star. However, the spectral type of B6Ve puts this star comfortably outside of the typical range of $\beta$ Cephei pulsators. A BeSS spectrum from 2020 March 31 shows $H\alpha$ in pure absorption, with no shell features. Three professional BeSS spectra with FEROS and ELODIE show weak $H\alpha$ emission, so the Be classification seems reasonable. This is a strange system given the LC and relatively late spectral type. Two minutes cadence data were used for this analysis.

TIC 190393155 = HD 75925: B4Vnne. The strongest feature in the LC is a group centered at $2.5\ \text{day}^{-1}$ (the strongest of these signals being at $2.535\ \text{day}^{-1}$, 2.46 ppt). Lower-frequency signals are found at $1.419\ \text{day}^{-1}$ (0.71 ppt) and $0.947\ \text{day}^{-1}$ (0.46 ppt). There is a group centered at $9.3\ \text{day}^{-1}$ (with its strongest two signals having an amplitude of $\sim 0.16$ ppt), a group at $13.5\ \text{day}^{-1}$ (strongest signal having an amplitude of $\sim 0.1$ ppt), a third group centered at $18.0\ \text{day}^{-1}$ (its two strongest signals at $17.831\ \text{day}^{-1}$ (0.16 ppt) and $18.162\ \text{day}^{-1}$ (0.15 ppt)), and a weaker fourth group at $\sim 22.7\ \text{day}^{-1}$ (with its strongest signal at 0.05 ppt). It is unclear if these four groups at higher frequencies are related to lower-frequency groups. This target is blended with the star SAO 220579, which is about 2 mag fainter and 13" distant (the contamination ratio is 0.1401). Included in the Renson & Manfroid (2009) catalog of Ap and Am stars, with a note saying that the star is of "doubtful nature." Two minutes cadence data were used for this analysis.

TIC 237651093 = HD 50696: B1.5IIe. There are typical groups centered at $2.16\ \text{day}^{-1}$, $4.31\ \text{day}^{-1}$, and $8.52\ \text{day}^{-1}$. There are apparently isolated signals at $6.403\ \text{day}^{-1}$ (0.41 ppt), $11.100\ \text{day}^{-1}$ (0.28 ppt), and $15.834\ \text{day}^{-1}$ (possibly a combination of $=3 \times f_{2.127} + 2 \times f_{4.722}$ or $=f_{2.361} + 3 \times f_{4.487}$, 0.29 ppt). It is unclear if this highest-frequency signal is a genuine combination frequency or if it is a coincidence. There is typical, symmetric, double-peaked $H\alpha$ emission in one BeSS spectrum from 2001, but more recent spectra do not show evidence of $H\alpha$ emission, consistent with a Be star with a variable disk. Two minutes cadence data were used for this analysis.

TIC 301435200 = HD 307350: B2Ve. At lower frequencies, there are two pairs of signals at $2.974\ \text{day}^{-1}$ (3.62 ppt) and $3.125\ \text{day}^{-1}$ (1.20 ppt), and $3.617\ \text{day}^{-1}$ (0.97 ppt) and $3.467\ \text{day}^{-1}$ (0.47 ppt). In both of these pairs, the two signals are separated by $0.15\ \text{day}^{-1}$. Further signals, not apparently in groups include $7.803\ \text{day}^{-1}$ (0.16 ppt), $9.176\ \text{day}^{-1}$ (0.17 ppt), $10.776\ \text{day}^{-1}$ (0.20 ppt), $12.152\ \text{day}^{-1}$ (0.16 ppt), $13.018\ \text{day}^{-1}$ (0.10 ppt), $15.276\ \text{day}^{-1}$ (0.12 ppt), and $16.537\ \text{day}^{-1}$ (0.09 ppt). Thirty minutes cadence data were used for this analysis since this

target was observed in two sectors (10 and 11), but 2 minutes cadence is only available for sector 10.

TIC 308951795 = HD 306145: B2Vne. High- and very-high-frequency signals exist between about 9 and $40\ \text{day}^{-1}$, including $9.779\ \text{day}^{-1}$ (0.71 ppt), $12.666\ \text{day}^{-1}$ (0.30 ppt), $29.993\ \text{day}^{-1}$ (0.01 ppt), and $40.020\ \text{day}^{-1}$ (0.01 ppt). The S/N of the last two frequencies is marginal ($\sim 4.2$, calculated with a $4\ \text{day}^{-1}$ window with period04). Possibly not a classical Be star; see Appendix B.4.

TIC 384471407 = HD 78328: B9.5IIe. Has an isolated very-high-frequency signal at $32.369\ \text{day}^{-1}$ (0.14 ppt), which is tied for being the strongest periodic signal along with a frequency at $0.558\ \text{day}^{-1}$ (0.14 ppt), which is surrounded by lower-amplitude signals. One low-resolution BeSS spectrum shows $H\alpha$ in weak emission.

TIC 399669624 = 2 Ori: A1Vne. Many spectra from BeSS show a weak, symmetric, double-peaked $H\alpha$ emission profile that is roughly stable over 13 yr. However, this is not a classical Be star, but is rather a $\lambda$ Boo star. Discussed further in Section B.2. There are multiple signals in the traditional g- and p-mode regimes, at $4.151\ \text{day}^{-1}$ (0.22 ppt), $9.143\ \text{day}^{-1}$ (0.16 ppt), $3.011\ \text{day}^{-1}$ (0.12 ppt), and $1.510\ \text{day}^{-1}$ (0.07 ppt). Multiple very-high-frequency signals include $16.559\ \text{day}^{-1}$ (0.05 ppt), $17.152\ \text{day}^{-1}$ (0.05 ppt), $21.303\ \text{day}^{-1}$ (0.04 ppt), $25.452\ \text{day}^{-1}$ (0.04 ppt), $50.414\ \text{day}^{-1}$ (0.01 ppt), and $50.910\ \text{day}^{-1}$ (0.01 ppt). Although low amplitude, the two highest-frequency signals have S/N = 9.2 and 7.3, respectively (using a $4\ \text{day}^{-1}$ window with period04). Two minutes cadence data were used for this analysis.

TIC 427400331 = HD 290662: A0Vpe. Has typical groups centered at 1.05 and $2.12\ \text{day}^{-1}$ (maximum amplitudes of 0.56 and 0.19 ppt, respectively), and a low-frequency signal at $0.248\ \text{day}^{-1}$ (1.42 ppt). The periodogram is largely featureless beyond this except for a group centered near $75\ \text{day}^{-1}$, composed of signals at $73.811\ \text{day}^{-1}$ (0.12 ppt), $74.100\ \text{day}^{-1}$ (0.28 ppt), $74.600\ \text{day}^{-1}$ (0.25 ppt), $74.817\ \text{day}^{-1}$ (0.16 ppt), and $75.050\ \text{day}^{-1}$ (0.63 ppt). Possibly not a classical Be star; see Section B.4. Two minutes cadence data were used for this analysis.

TIC 451280762 = HD 99146: B3Ve. Has a relatively complex frequency spectrum with low-frequency stochastic excess and numerous periodic signals. This complexity makes it difficult to determine if groups are present. High- and very-high-frequency signals are found at $8.892\ \text{day}^{-1}$ (0.19 ppt), $10.018\ \text{day}^{-1}$ (0.46 ppt), $10.609\ \text{day}^{-1}$ (0.22 ppt), $11.771\ \text{day}^{-1}$ (0.09 ppt), $15.494\ \text{day}^{-1}$ (0.07 ppt), and $15.546\ \text{day}^{-1}$ (0.07 ppt). The frequency spectrum is remarkably rich. Discussed further in Appendix B.4. Two minutes cadence data were used for this analysis.

TIC 467027607 = HD 306111: Oe. Has low-frequency signals at $0.171\ \text{day}^{-1}$ (3.87 ppt) and $0.408\ \text{day}^{-1}$ (2.24 ppt), and two typical groups centered at 1.8 and $3.7\ \text{day}^{-1}$. Further high and very high frequencies are found at $8.745\ \text{day}^{-1}$ (0.37 ppt), $9.614\ \text{day}^{-1}$ (0.27 ppt), $12.011\ \text{day}^{-1}$ (0.14 ppt), $18.577\ \text{day}^{-1}$ (0.25 ppt), and $18.809\ \text{day}^{-1}$ (0.35 ppt). Two minutes cadence data were used for this analysis.

## B.2. Firmly Nonclassical Be Stars Rejected from the Sample

TIC 14709809 = RY Gem: A known "moderately interacting Algol-type eclipsing binary" (Plavec & Dobias 1987). TESS data likewise indicate this is an eclipsing binary (EB) with a period of about 9 days, where the primary eclipse reaches a
~60% depth. Many BeSS spectra show $H\alpha$ with a clear double-peaked profile and shell absorption, with variable V/R ratios, likely formed in an accretion disk around the A0V (present-day) primary.

**TIC 30562668 = HD 76838:** Short-period EB (~4 days) with $\beta$ Cephei pulsations and reflection effect. Embedded in a nebula as seen in WISE images, which may have been formed in a past binary interaction that resulted in the current binary configuration.

**TIC 53063082 = HZ CMa:** The TESS light curve looks much more like rotation or binarity with $P$ = 15 or 30 days. Strong single-peaked $H\alpha$ emission in a single high-resolution BeSS spectrum from 2009, with an emission level relative to the continuum (E/C) ~2.3. This system is included in many studies of Be stars (Sterken et al. 1996a; Zhang et al. 2005; Touhami et al. 2011; Labadie-Bartz et al. 2017; Wang et al. 2018; Klement et al. 2019). This is a binary with an accretion disk and ellipsoidal variability with $P_{orb}$ = 28.6 days, and apparently has a Roche-lobe overfilling K giant companion (Sterken et al. 1994), where the authors interestingly note that the high vsini value (300 km s⁻¹; Meisel 1968) of the B star may be a consequence of accretion; i.e., this could be a Be star in the making (but currently has an accretion disk as it is being spun up). Listed as B6IVe+A in Slettebak (1982). Besides the orbital modulation, the TESS data contains only stochastic variability.

**TIC 90296023 = FY Vel:** KELT data show this to be a 33.75 day period EB ellipsoidal variable, as reported in Labadie-Bartz et al. (2017), and is a known Beta Lyrae type EB (Thackeray et al. 1970). Not a classical Be star. Presumably there is mass transfer via Roche-lobe overflow, leading to the strong $H\alpha$ emission seen in the four available BeSS spectra, which show $H\alpha$ E/C ~8–10.

**TIC 123545883 = V743 Mon:** An unclassified B[e] star, or potentially a Herbig Ae/Be star (Varga et al. 2019), where severe low-frequency stochastic variability is its main feature in TESS.

**TIC 141973945 = bet Hya:** The TESS light curve resembles rotation or binarity with a dominant signal at 0.43 day⁻¹ (12.5 ppt), and lower-amplitude harmonics. Four BeSS spectra all show a shell signature with weak emission. Simbad lists this as an alpha2 CVn type variable. A spectral type of kB8hB8HeA0VSi is given in Garrison & Gray (1994), and is therefore not a classical Be star.

**TIC 151300497 = V1075 Sco:** Low-frequency stochastic variability is prominent in TESS, as are longer-term trends. There are narrow ill-defined groups that stand out above the noise. Not a classical Be star, despite $H\alpha$ emission in BeSS spectra. The GOSS survey gives a spectral type of O7.5V((f))z (e) (Sota et al. 2014).

**TIC 151131426 = HV Lup:** There is one dominant signal with a double-wave pattern at $P$ = 5.66 days with a ~10% semi-amplitude, plus low-frequency stochastic variability. The TESS LC suggests that this is an interacting binary with a 5.66 day orbital period. The behavior of the light curve and frequency spectrum most closely resembles other interacting binaries rather than a classical Be star. No BeSS spectra are available. An HARPS spectrum downloaded from the ESO database (http://archive.eso.org/cms/data-portal.html) shows $H\alpha$ with a P Cygni profile, and also similar but weaker emission in $H\beta$. Most of the helium lines seem to have two components as expected for an SB2 with similar effective

temperatures. Considering the TESS and HARPS data together, this appears to be a binary and not a classical Be star.

**TIC 178719204 = HD 70340:** The light curve resembles rotation or binarity with a double-waved morphology and a period of 12 days. Thirty-one BeSS spectra taken between 2009 and 2021 show no sign of emission. A possibly chemically peculiar early A-type star (Ghazaryan et al. 2018).

**TIC 187458882 = HD 57682:** Many spectra in BeSS show the $H\alpha$ emission as being single peaked, very variable in strength, narrow, and seemingly variable in position relative to the main $H\alpha$ absorption. The MiMeS survey find this to be a strongly magnetic O9.5 IV star (Grunhut et al. 2017), and is therefore incompatible with a classical Be (or Oe) designation. Its only feature in TESS is strong low-frequency stochastic variability. There is also a diffuse nebulae surrounding the star in the WISE red band, and a strong SED excess at the longest radio wavelengths.

**TIC 203452834 = HL Lib:** A binary system with spectral types G6: III + B8.5 V and an orbital period of 24.615 days (Dempsey et al. 1990), where the authors state “Photometric observations indicate random changes superimposed on regular ellipsoidal light variations, the latter probably the result of tidal distortion of the giant primaries. Mass transfer and loss is apparent in inverted mass ratios derived from orbital analysis, strong wind features present in the spectra, and the presence of circumsystem shells.” The variability in TESS is mainly a double-waved signal at the orbital period and low-frequency stochastic variability.

**TIC 207176480 = HD 19818:** High-resolution spectroscopy acquired by our group shows this is an SB2, with an A0 main-sequence star and a cooler giant. Rotational modulation is apparent in TESS and KELT (Labadie-Bartz et al. 2017) with the same period of 3.3 days. Hydrogen and other emission features are variable on short timescales and seem incompatible with a Keplerian disk. We suggest that, due to the rapid rotation of the cool giant star, a strong magnetic dynamo exists causing the observed rotational modulation (due to spots), and the highly energetic flares seen in the TESS data (Labadie-Bartz & Carciofi 2020b). Chromospheric activity is likely the cause of the transient emission features (and the strong X-ray flux).

**TIC 220322383 = 15 Mon:** 15 Mon is a binary with an O7V ((f))z primary (Sota et al. 2014) in a ~25 yr eccentric orbit, with a secondary that is probably a very late O or very early B spectral type (Gies et al. 1993). The system is a known X-ray source (Ramírez et al. 2004; Nazé 2009) and is a cluster member (NGC 2264). The signals in TESS are of low amplitude, the strongest of which is 0.5 ppt at 12.5 day⁻¹ (with a nearby signal at 11.9 day⁻¹), which may be p-mode pulsation in the secondary.

**TIC 224244458 = bet Scl:** The TESS data more closely resemble rotational modulation with a period of about 2 days, with a spectral type of B9.5IIpHgMnSi (Abt & Morrell 1995). Rejecting as a Be star, this is a chemically peculiar star exhibiting rotational modulation.

**TIC 234813367 = AX Mon:** Puss & Leedjärv (2002) classified this as an interacting binary system with a K giant and a B(e) star that is accreting matter, noting the P Cygni profile and unusual emission variations over the orbit. This system appears in numerous studies of populations of Be stars, but an interacting binary with an accretion disk around the hot component is incompatible with the decretion disks of classical Be stars. The TESS light curve is dominated by low-frequency
stochastic variation, but there are clear periodic signals, including one that stands out near $3.5 \text{ day}^{-1}$. Many BeSS spectra show a clear P Cygni profile in $H\alpha$.

TIC 234835218 = EM\* GGA 395: A Herbig Ae/Be star (Li et al. 2002). It is very clearly in a strong nebula as seen in WISE images, and the SED indicates a large amount of cool dust. It is characterized as having strong low-frequency stochastic variability and an isolated frequency near $2.6 \text{ day}^{-1}$ in TESS.

TIC 238791674 = CD-49 3441: A Herbig Ae/Be star according to The et al. (1994). High-amplitude low-frequency stochastic variability is the most prominent feature in TESS. A 6.77 day signal is found in the KELT data (Labadie-Bartz et al. 2017). Three low-resolution BeSS spectra show strong $H\alpha$ emission about seven times the continuum level (which would be abnormally high for a B8 classical Be star).

TIC 246189955 = HD 328990: The SED is incompatible with that of a classical Be star, suggesting two stellar components (likely a cool giant, and a late B or early A main-sequence star). There is only stochastic variability present, and the variability is dominated by low frequencies. The amplitude is very high (~10% max–min) relative to other confirmed classical Be stars of its spectral type, A0e.

TIC 253212775 = V495 Cen: An eclipsing interacting binary with a 33.5 day orbital period; a cool evolved star, and a hot mid-B dwarf with an accretion disk (Rosales Guzmán et al. 2018). The TESS light curve is dominated by the orbital modulation, but there are also clear (possibly aperiodic) oscillations with amplitudes of about 1% and frequencies between 0.5 and $3 \text{ day}^{-1}$.

TIC 253380837 = HD 113573: The TESS data unambiguously show this to be an EB with an orbital period of about 1 day, and primary and secondary eclipses that are slightly different in depth, but both around 3%. Although the field is somewhat crowded, the use of differently sized apertures suggests the target star is the source of the EB signal. Besides the orbital period and its many harmonics, no other signals are present in the data.

TIC 289877581 = d Lup: The TESS variability looks much more like rotation than pulsation. Given the “p” in the spectral type (B3IVpe), this may have surface spots (chemical or magnetic) leading to modulation at the $0.48 \text{ day}^{-1}$ signal—the strongest signal in the LC (assumed to be rotation). This star is included in catalogs of chemically peculiar stars (Romanyuk & Kudryavtsev 2008), supporting the plausibility of rotationally modulated brightness. The single BeSS spectrum is narrow-lined in $H\alpha$, with no sign of emission. Arcos et al. (2018) analyzed two epochs of spectra from 2013 and 2015, saw no variability or emission in $H\alpha$, and fit blue photospheric absorption lines to arrive at a projected rotational velocity of $v\text{sini} = 30 \text{ km s}^{-1}$, which would be unusually slow for a rapidly rotating classical Be star.

TIC 305090822 = HD 157273: Just one long sinusoidal signal exists, with $P = 15.72 \text{d}$ (also very clear in KELT data, unpublished). Clearly comprising emission in $H\alpha$ in the three low-resolution BeSS spectra. The SED appears to be formed from two components. This is probably an interacting binary with an accretion disk and possibly ellipsoidal variation.

TIC 320228013 = HD 308829: The TESS data show this to be a short-period EB ($P_{\text{orb}} \sim 6.77 \text{ days}$) with significant asymmetric out-of-eclipse variability possibly caused by some combination of a reflection effect and rotation of an

inhomogeneous stellar surface. This system is a cluster member (Cl\* IC 2944 THA 51). The relatively dense region of the sky makes blending in TESS problematic, so it is not yet certain that the EB signal can be attributed to HD 308829. The SED has peculiarities. This is a known X-ray source (NRS2013), included in the ROSAT all-sky bright source catalog (Voges et al. 1999) and also observed with XMM-Newton, as discussed in Nazé et al. (2013).

TIC 323612875 = HD 102369: TESS data show a very simple frequency spectrum with only one low-frequency ($\sim 0.45 \text{ day}^{-1}$) signal and a weak first harmonic, which resembles rotation as opposed to pulsation. No spectra are available on BeSS or BeSOS (http://besos.ifa.uv.cl).

TIC 333670665 = V863 Cen: A known magnetic He-strong star (Shultz et al. 2019). The TESS variability looks more like rotation or binarity than pulsation with a period of about 1.3 days, and a weak harmonic. Three BeSS spectra show narrow lines with no sign of emission from 2012, 2017, and 2018.

TIC 380117288 = AI Cru: A short-period EB ($P_{\text{orb}} \sim 1.4 \text{ days}$) with $\sim 50\%$ primary depth. Previously known to be an eclipsing binary (Kreiner 2004).

TIC 381641106 = CSI-62-12087: Listed as a WR star, MR 41 (Roberts 1962). The TESS light curve is purely stochastic at a high amplitude and low frequency.

TIC 394728064 = DR Cha: A clear EB with primary and secondary eclipses separated by about 20 days. No BeSS spectra exist. None of the six entries in the catalog of Skiff (2009) indicate emission. It is unclear what the reason for the historical Be star classification of this system, but perhaps its binary nature has lead to confusion.

TIC 399669624 = 2 Ori: A1Vne. This is not a classical Be star, but is rather a $\lambda$ Boo star (Abt & Morrell 1995; Murphy & Paunzen 2017); they have peculiar abundance patterns and often pulsate in $\gamma$ Dor and/or $\delta$ Scuti modes. Many spectra from BeSS show a weak, symmetric, double-peaked $H\alpha$ emission profile that is roughly stable over 13 yr, which, given the $\lambda$ Boo classification is probably a circumstellar accretion or debris disk or similar. This system has a high rotation rate, with $v\text{sini} = 211 \text{ km s}^{-1}$ from Glebocki & Gnacinski (2005), or $v\text{sini} = 261 \text{ km s}^{-1}$ from Ammler-von Eiff & Reiners (2012).

TIC 408757239 = V716 Cen: Short-period ellipsoidal variable and EB, with $P_{\text{orb}} \sim 1.5 \text{ days}$. A known EB (Kreiner 2004).

TIC 442240473 = HD 84511: An EB with period $= 32.996 \text{ days}$ (Labadie-Bartz et al. 2017) and strong ellipsoidal variation (also seen in TESS).

TIC 443616529 = *phi Leo*: This system may have a debris disk, as there is some evidence for exocomets in the system (Eiroa et al. 2016). With a spectral type of A7IVne, this is too late to belong to the class of classical Be stars. The strongest signal by far is an isolated frequency at $6.48 \text{ day}^{-1}$, plus some lower-level variability that may be stochastic, being strongest at the lowest frequencies. It is unclear why this is in lists of Be stars, although a debris disk may cause peculiarities (shell features) that can be confused with a gaseous decretion disk. No BeSS spectra from 2007 to 2020 show any sign of emission, nor do the three spectra on BeSOS (http://besos.ifa.uv.cl). This is a high proper-motion star.

TIC 450276053 = V338 Car: Very clearly an EB or ellipsoidal variable in KELT with $P_{\text{orb}} = 73.8 \text{ days}$ (unpublished).
**TIC 455463415 = HD 135160:** A very strange EB with a strong reflection effect, a primary eclipse depth of about 2%–2.5%, and an orbital period of about 6 days. There are odd “bumps” in the LC not synchronized with the orbit, repeating roughly every 8 days. This is a known SB2, confirmed in Wang et al. (2018). None of the references or spectral types in the eight entries in the Skiff (2009) catalog include mention of emission, so the reason for this being included in the Jaschek & Egret (1982) catalog of Be stars is unclear.

**TIC 467065657 = HD 97253:** Prominent low-frequency stochastic variability with high amplitude. Classified as having a spectral type O5III(f) from the GOSS survey (Sota et al. 2014), which is incompatible with being a classical Be star.

### B.3. Possible Nonclassical Be Stars Rejected from the Sample

**TIC 11411724 = StHA 52:** Embedded in a strong reflection nebula, NGC 2023 (Rozhkovskij & Kurchakov 1968). Unusual frequency spectrum. Four frequencies of similar strength spread out over 0.5 day⁻¹, centered at 8.75 day⁻¹, and other high frequencies between 10 and 11.5 day⁻¹. Not listed in BeSS, and has no emission features in APOGEE (ABE-A42) instead showing only broad absorption. Given the strength and size of the reflection nebula in the visible, it is unclear exactly what is contributing to the TESS photometry.

**TIC 22825907 = HD 148877:** Dominated by a 10 day period with an amplitude of 4%. Looks like rotation and/or binarity + stochastic variability. KELT data (unpublished) show a strong sinusoidal signal at 9.95 days, which is possibly double-waved at 19.91 days. The SED appears to have two components. This is more likely a mass-transfer binary than a classical Be star. The stochastic variability peaks at a frequency of about 3 day⁻¹.

**TIC 23091719 = NW Pup:** B2IVne. Lower-frequency variability, including stochastic features, dominates the TESS LC, but no obvious groups exist. No sign of emission, and narrow-lined, including a narrow He 6678 line in one high-resolution BeSS spectrum from 2020-04-15. Unclear if this is truly a classical Be star given the narrow lines and lack of emission in 2019–2020. No sign of emission in 2014 in two epochs of BeSOS spectra (http://besos.ifa.uv.cl), and their best-fit model gives $v\text{sini} = 50 \text{ km s}^{-1}$. Far-radio excess in the SED. Included in the catalog of Egret & Jaschek (1981), where it is listed as an He-abnormal chemically peculiar star, and in the Renson & Manfroid (2009) catalog of Ap and Am stars, being listed as B3 He var. BeSS spectral type and emission-line designation seems to come from only Hiltner et al. (1969). Very narrow-lined in the blue in Chauville et al. (2001).

**TIC 26175330 = 17 Sex:** A1Ve. No signals in TESS. Extremely deep and narrow $H\alpha$, dropping to about 0.175 times the continuum level. No sign of any emission in many BeSS spectra, including a professional spectrum from ELODIE. Listed as a Herbig Ae/Be star on Simbad, and classified as an A-shell star in Montesinos et al. (2009). Likely falsely classified as a Be star because of shell features in this A-type star.

**TIC 140132301 = HD 72126:** Extremely strong, stable, single frequency at 3.04 day⁻¹, with a semi-amplitude of about 10%. Very unusual LC and frequency spectrum for a Be star. No other Be stars in the sample have such a simple frequency spectrum with only one extremely strong frequency. There is some diffuse nebula in the vicinity, so perhaps this is misclassified as a Be star. There are no BeSS spectra. SED

does not appear to show any fluctuation in the IR/radio. It is in the Be star catalog of Jaschek & Egret (1982), but is not included as a Be star here.

**TIC 213153401 = HD 154538:** The TESS LC looks similar to that of Sigma Ori E (Landstreet & Borra 1978). Two BeSS spectra are too low-resolution to provide a detailed view, but they are both clearly in absorption. This system could be an EB, but without any ellipsoidal variation and a very short period of $P = 0.5$ day and equal primary/secondary depths (or no secondary at all), or eclipses could be caused by magnetically confined corotating clouds (like in Sigma Ori E). Very close visual double, so it is hard to say from which star the signals are coming from. It is possible that the fainter of the pair is the one with emission (see Skiff 2009). At any rate, the LC is too contaminated to analyze in the context of this work.

**TIC 269087549 = 19 Mon:** Possibly a hybrid SPB/$\beta$ Cephei pulsator, with a pair of $\beta$ Cephei-like pulsations centered around 5 day⁻¹, with a combination frequency near 10 day⁻¹. Known $\beta$ Cephei star in Simbad. H and He line profiles on BeSS look significantly deformed, perhaps from line-profile variations from pulsation, and/or binarity. No sign of any emission from 2003–2020 in hundreds of spectra in BeSS. It is in the Jaschek & Egret (1982) catalog of Be stars. In the Skiff (2009) catalog, there are 11 entries, and only one of them hints at emission, with the note “em?” from Irvine (1975). Found to be an SB2 based on nine spectra (Chini et al. 2012). Assuming this is not a classical Be star, and that the potential “em?” classification from 1975 was the result of being an SB2.

**TIC 284230347 = HD 55806:** High-amplitude low-frequency variation, 3% semi-amplitude, plus low-frequency stochastic variation. Apparently located in the seismology fields of CoRoT (Gutiérrez-Soto et al. 2007), although to date, no CoRoT data have been reported on for this object. Frémat et al. (2006) noted for HD 55806 that they could not find any set of fundamental parameters allowing a simultaneous fit of the observed He and Mg spectral lines. All of the He lines show unusual line shapes, probably related to the presence of a close companion. It is the only star in their sample of 64 to show this behavior. The low-frequency signal is also very apparent in the KELT data, being single-waved at 6.996 days or maybe double-waved at 13.98 days (Labadie-Bartz et al. 2017). This suggests a close binary, and perhaps an accretion disk scenario. The TESS data seem more similar to other interacting binaries compared to “normal” classical Be stars.

**TIC 319854805 = HD 47359:** A suspicious and interesting case with features that look like simple flickers with precursor phases and no apparent change in the frequency spectrum, and a pair of high-frequency signals at 11.8 and 13.3 day⁻¹. Just one dominant frequency and a very-low-amplitude harmonic of it, plus obvious low-frequency signals are present. Three BeSS spectra show weak $H\alpha$ emission, sometimes double-peaked, sometimes single-peaked, with $E/C \sim 1.2$. Both the 12 day “flicker” signal and the $f = 1.545$ day⁻¹ are very clear in KELT (Labadie-Bartz et al. 2017). The $P = 12.13$ days signal is interesting with how regular it is in KELT and how much it resembles the TESS signal, and maybe is related to rotation or binarity. This probably is not a normal Be star flicker. Listed as B0.5Vp in (Yudin 2001). With $v\text{sini} = 443 \pm 40 \text{ km s}^{-1}$ (Frémat et al. 2006), this would have to be very close to edge-on, yet the $H\alpha$ profile is clearly not that of a shell star. The
nature of this object is unclear, and it is not included in the statistics of this sample.

TIC 322104948 = HD 306989 = V644 Cen: Suspicious LC. Clear 25 day period in KELT (Labadie-Bartz et al. 2017), and also apparent in TESS, but also with low-frequency stochastic variability. Many high-frequency signals in TESS are evident. No BeSS spectra. SED is structured in a way that suggests two components. Young stellar object (YSO) candidate from Marton et al. (2016). Possibly a very long-period (200 yr) eclipsing binary with an eclipse duration of 17+ yr according to O'Connell (1951) and Davies et al. (1987). It is likely not a classical Be star.

TIC 342257745 = HD 322422: Unusual LC with an isolated low-frequency signal dominating ($f = 0.39 \text{ day}^{-1}$), probably a weak group at $1.61 \text{ day}^{-1}$, and some high-frequency signals at $6.91, 8.3$, and $13.4 \text{ day}^{-1}$. The low-frequency signal is not apparent in the lower-precision KELT data (unpublished). Embedded in a nebula. Listed as having emission in seven references in Skiff (2009). $v\text{sini} = 170 \pm 11 \text{ km s}^{-1}$ from Zorec et al. (2016). Possibly a Herbig Ae/Be star, listed in the extreme emission-line objects table (The et al. 1994); these stars are usually most likely LBVs, B[e] stars, HAEBEs, PNs, or Symbiotics.

TIC 376077639 = V862 Ara: B7IIIe. A clear case of sum and difference frequencies and harmonics of the main group. There is no emission in the two BeSS spectra with high enough resolution. The Be designation does not seem to have much evidence. The only reference for being an emission-line star at all is from Hipparcos photometry, where it earned the VSX designation “BE:” in which the colon indicates uncertainty (Samus et al. 2017). There does not seem to be reasonable evidence to claim this is a Be star.

TIC 455809360 = CD-61 4751: Suspicious LC dominated by low-frequency stochastic variation, with an SED that looks like two components. Behavior is more reminiscent of an interacting binary compared to that of Be stars, but it is unclear what the nature of this system is. No BeSS spectra are available.

TIC 457546452 = HD 126986: Low-frequency periodic double-waved LC with 4% amplitude is the main feature seen in TESS—almost definitely binary or rotation related (with a 7.2 or 14.4 day period). Unlikely to be a classical Be star just based on TESS data. No BESS spectra are available. The SED seems to have two components.

## B.4. Possibly Nonclassical Be Stars Kept in the Sample

TIC 21249978 = V868 Ara: Unusual LC. Strong sinusoidal signal with period of 22.51 days (possibly double-waved at twice this period). Therefore unlikely to be a classical Be star. One high-resolution BeSS spectrum shows $H\alpha E/C \sim 3$, with a deformed triple-peaked profile, and $H\beta$ is mostly filled in. The SED has considerable structure and may indicate two components. This may be either an ellipsoidal variable (but not eclipsing) with mass transfer and an accretion disk, or a classical Be star with a decretion disk in a 45 day (or 22.5 day) binary. There seem to be two frequency groups centered around $4.5$ and $9.0 \text{ day}^{-1}$, plus the low-frequency orbital modulation and low-frequency stochastic variation. It has a somewhat high proper motion, of $12.772$ and $-6.280 \text{ mas yr}^{-1}$ (Gaia Collaboration 2018).

TIC 51288359 = HD 151083: EB with a $0.901$ day (or $1.8$ days) period, plus higher-frequency signals near $6.25 \text{ day}^{-1}$

and $11.5 \text{ day}^{-1}$, and some low-frequency stochastic variability. Eclipse depths are approximately 2%–3%. While such a short-period EB would be incompatible with a classical Be star, a situation where such an EB is blended with a Be star cannot be ruled out; there is little information on this target available in the literature, and it lacks BeSS spectra.

TIC 140031673 = HD 71510: B2Ve. There appear to be groups, but their pattern and delineation is unclear. None of the spectra in BeSOS (2014, 2015; http://besos.ifa.uv.cl) or BeSS (2014/15/17/19) show $H\alpha$ in emission. $v\text{sini} = 150 \text{ km s}^{-1}$ from BeSOS. It is in the Jaschek & Egret (1982) catalog of Be stars and is apparently a known visual binary (e.g., Wackerling 1970). None of the six references in Skiff (2009) list emission. The visual binary pair is listed as B3V + G3V (Pallavicini et al. 1992), but that does not exclude the B-type star being a Be star. There is an IR nebula centered on the target (Bodensteiner et al. 2018). This star is of unclear nature. Perhaps the IR excess from the surrounding nebula is the reason for its Be designation, but without further evidence, this cannot be excluded as a classical Be star. This is included in the statistics for the sample. Variability in TESS is largely stochastic, but there do seem to be coherent signals in the traditional g-mode regime.

TIC 147244857 = HD 70234: Unusual LC for a Be star. Has only two isolated frequencies with constant power at $1.55$ and $2.27 \text{ day}^{-1}$ causing a beating pattern with a very short envelope. First reported as an emission-line star in Henize (1976), where it is listed as B9III(e). One low-resolution BeSS spectrum clearly shows $H\alpha$ in emission at $E/C = 2.5$.

TIC 215573117 = CD-27 5181: This looks more like an EB or a Sigma Ori E type case, with two strong signals being harmonically related (at $0.6$ and $1.2 \text{ day}^{-1}$), which resemble eclipses of about 1% depth, plus low-frequency stochastic variability. No BeSS spectra are available. Included in the Be star catalog of Jaschek & Egret (1982), and is apparently of a relatively early spectral type, being classified as OB+e in Stephenson & Sanduleak (1971). The same photometric frequency is reported in Labadie-Bartz et al. (2017). Without a spectrum, we cannot determine whether or not this is a classical Be star, and is thus not rejected from the sample.

TIC 191312952 = HD 129772: One dominant signal at $3.15 \text{ day}^{-1}$, which is considered to be a frequency group given the apparent beating pattern in the LC. There are multiple harmonics of this signal/group. It is in the Renson & Manfroid (2009) catalog of Ap, HgMn, and AM stars, where it is listed as having a spectral type of B8 Ca, suggesting an abundance of calcium, but its classification is also of “doubtful nature.” No BeSS spectra are available. There is insufficient evidence to rule this out as a classical Be star.

TIC 207580161 = HD 119835: This is listed in the catalog of YSO candidates (Marton et al. 2016), and there are no signals in the TESS data (apart from some mild systematics), which is not too unusual for its late spectral type (B8.5IIIne). However there does not seem to be any evidence pointing to this not being a classical Be star.

TIC 215983126 = HD 144970: The star is very reddened, but is apparently B0V or B0Ve (Feast et al. 1961; Cannon & Pickering 1993). No BeSS spectra are available. TESS data show only low-frequency stochastic variability. It could be some other type of object, but this is inconclusive and so is included as a Be star in this work.
TIC 256994805 = V715 Mon = HD 49567 = HR 2517: Only a weak disk is present in one professional BeSS spectrum. This system is apparently a high-mass X-ray binary (Khalak et al. 1998) and has also been observed to show “flares” (Sterken et al. 1996b), which are strange since they are brighter in the bluer bands (Be outbursts are generally of higher amplitude in redder bands). TESS sees two fairly typical groups (but at relatively low frequencies, $0.50 \text{ day}^{-1}$ and $1.0 \text{ day}^{-1}$), and an isolated signal at $\sim 5.6 \text{ day}^{-1}$.

TIC 282808223 = HD 50820 = HR 2577: A binary B3IVe +K2II system with a period of 58 yr where the Balmer emission lines of the B star are variable in a fashion unrelated to the orbit (Hendry 1982). The TESS LC shows mostly slow, gradual, low-amplitude variation plus mild low-frequency stochastic variation and a group near $1 \text{ day}^{-1}$. Included as a Be star in this work, although it is possible that variability from the K giant contributes to the TESS data.

TIC 308951795 = HD 306145: Unusual LC, with many short, small “bumps” that are not obviously related to pulsation or flickers and are seemingly oscillatory although do not appear to be strictly periodic. Although these bumps may be morphologically similar to flickers, their very short timescales of $\sim 1$ day (from start to finish) and high occurrence rate compared to the majority of flicker events may indicate that they represent something other than discrete mass ejection episodes. There are also definitely high and maybe very high frequencies, and low-frequency stochastic variation out to high frequencies. No emission in the single low-resolution BeSS spectrum. Very close visual double (clearly seen in the Two Micron All Sky Survey). The strongest signal is at $0.25 \text{ day}^{-1}$, or $P = 3.9$ days (also apparent in KELT, unpublished). The variability in TESS is very unusual, but cannot be ruled out as a classical Be star based on existing data.

TIC 315679257 = HD 146596: Primarily low-frequency signals ($\sim 1\%$ in amplitude) and stochastic variability. No BeSS spectra are available. It has $H\alpha$ emission from HARPS according to Rainer et al. (2016), and is therefore classified as a Be star, but without estimates of stellar parameters. The SED shows an excess in the mid-IR with a clear down-turn. The $H\alpha$ line (from http://sisma.brera.inaf.it/index.jsp) looks typical for a Be star, with $E/C \sim 2.3$, and $V/R \sim 1$. Possibly some weak and very narrow emission in other non-Balmer lines (near 4922 Å, 5020 Å, and some others), which appear to have P Cygni profiles. The nature of this object is unclear, but cannot be ruled out as a classical Be star.

TIC 316792722 = HD 99771: The two lowest-frequency groups appear to be in the typical configuration, at 0.61 and $1.23 \text{ day}^{-1}$. There is a third group near $2.19 \text{ day}^{-1}$, and then many isolated frequencies with exact harmonics. It is possible that the signal near $4.96 \text{ day}^{-1}$ is split, and there may be a slightly unusual SED excess at long wavelengths. This system is embedded in a large cloud. Included in the Jaschek & Egret (1982) catalog of Be stars, but there do not seem to be any more recent works or spectra that are relevant. There are no signals in KELT (unpublished), as expected, due to their low TESS amplitudes (maximum of 0.4 ppt).

TIC 322233181 = HD 306962: Dominated by low-frequency stochastic variability, and an isolated signal at $5.97 \text{ day}^{-1}$. This system was selected as a YSO candidate based on broadband WISE photometry (Marton et al. 2016). There is not enough evidence to discount this being a classical Be star, but it is suspect.

TIC 381747495 = HD 105753: Odd LC. Dominated by low-frequency stochastic variation, and a somewhat strong group near $8.5 \text{ day}^{-1}$. Clear long-term trends in KELT data (unpublished). There is no information to rule this out as a classical Be star, but the light curve and frequency spectrum are suspicious. No similar cases are seen in the sample. No BeSS spectra.

TIC 405520863 = 39 Cru: Looks more like rotation or binarity. Only a single signal at $1.30 \text{ day}^{-1}$, and a few weak harmonics. Two BeSS spectra have $H\alpha$ emission with $E/C \sim 3.5$, and a symmetric double peak, and two spectra on BeSOS (http://besos.ifa.uv.cl) from 2014 and 2015 have the same profile. Unusual in its simplicity and constant $H\alpha$ emission. There is no obvious evidence against this being a classical Be star, but its TESS light curve is remarkably simple.

TIC 427400331 = HD 290662: Has a very-high-frequency group at $75 \text{ day}^{-1}$ with 0.6 ppt amplitude, plus typical lower-frequency groups centered at $1.05$ and $2.12 \text{ day}^{-1}$ and with amplitudes of around $0.25-1.5$ ppt. This is listed in the Renson & Manfroid (2009) catalog of Ap and Am stars, where this high-frequency group could then plausibly be roAp or $\delta$ Scuti pulsation. Renson & Manfroid (2009) gives a spectral type of B9 Fe. One low-resolution BeSS spectrum does not obviously show any emission. Skiff (2009) also notes an A0Vp designation from 1971 saying that Fe II lines prominent. It is unclear what the nature of this star is, and there is insufficient evidence to reject this as a Be star.

TIC 440399815 = HD 113605: Low-level low-frequency stochastic variability is the main feature in TESS. In KELT there is a fairly strong periodic signal at $P = 32.087$ days, which could be related to binarity (unpublished). Without any spectra, the nature of this object is unclear, but there is no convincing evidence against this being a classical Be star.

TIC 451280762 = HD 99146: Remarkably rich frequency spectrum in TESS, with high levels of low-frequency stochastic variation, and many clearly periodic signals out to $\sim 16 \text{ day}^{-1}$ that are not organized in a way that resembles the typical groups. No BeSS spectra are available. There is clear long-term variation and possibly outbursts in KELT (unpublished). The frequency spectrum is unusual, but there is no evidence to this not being a classical Be star.

TIC 466715331 = HD 308217: The TESS LC is dominated by low-frequency stochastic variation and low frequencies, but there is a frequency “group” near $6.2 \text{ day}^{-1}$, which looks like p-modes, and also an isolated signal at $8.44 \text{ day}^{-1}$. There are events that vaguely resemble flickers (with a quasi-period of about 10 days), but could be related to rotation. The single BeSS spectrum is in absorption. We cannot rule this out as a classical Be star.

## B.5. Be Star Systems of Interest

TIC 14088298 = HD 33453: B8Vne. Two very narrow groups with strong beating patterns. Groups are at $3.28$ and $6.53 \text{ day}^{-1}$, which is one of the widest separations in the whole sample. Normal looking symmetric $H\alpha$ spectrum in BeSS. The extreme group separation may make this interesting to study (perhaps it has a relatively high rotational frequency).

TIC 65803653 = 27 CMa: Looks like a composite spectrum with an isolated signal at $0.76 \text{ day}^{-1}$ and its first harmonic, plus two groups at $1.33$ and $2.69 \text{ day}^{-1}$. There are also isolated signals around $10.9$ and $13.5 \text{ day}^{-1}$, plus a pair of signals near $5.16$ and $5.90 \text{ day}^{-1}$. It is a shell star with a strongly varying
disk, including strong asymmetries from seven BeSS spectra between 2006 and 2019. It is a good candidate for binarity and a known β Cephei pulsator (Stankov & Handler 2005). Balona & Krisciunas (1994) found the same β Cephei frequency—they report 10.893 day⁻¹, and a low-frequency signal of 0.796 day⁻¹. In fact, this is a very close visual double where the two components are separated by ~0.11, and differ in magnitude by about 0.5 (Mason et al. 2009), so it is plausible that both of these stars contribute to the observed variability.

**TIC 118842700 = QV Tel:** B3IIIpe. One strong low-frequency signal (around 0.1 day⁻¹), and strong low-frequency stochastic variability. There are no obvious patterns in the LC. From BeSOS (http://besos.ifa.uv.cl), $Hα$ is symmetric and double-peaked with $E/C \sim 1.5$. There are two narrow C II lines to the right of $Hα$. Arcos et al. (2018) estimated vsini = 50 km s⁻¹, which even seems too high comparing the fit to the data. $Hα$ is definitely variable in the 17 BeSS spectra from 2011–2020. This is HR 6819, the potential triple (BH + B) + Be system (Rivinius et al. 2020), which can explain the narrow lines since the giant B star (with narrow lines) contributes about 50% of the total visible flux. Alternatively, this system may be a binary where the narrow-lined B star is in a short-lived evolutionary phase where it is currently contracting toward being a typical sdO star, having recently donated mass to the Be star (Bodensteiner et al. 2020b).

**TIC 127493611 = omi Pup:** BIVnne. Strong, clear, coherent sinusoidal signal at $f = 0.495$ day⁻¹ (amplitude ~6 ppt), and a few high-frequency signals near 8.5, 10.2, and 14.5 day⁻¹ (amplitudes between 0.1 and 0.3 ppt). $Hα$ is almost flat-topped in BeSS, but slanted down to the blue. $Hα$ does not change much over 2 yr in BeSOS spectra (http://besos.ifa.uv.cl). It is likely a Be + sdO binary with a 28.9 day period (Koubský et al. 2012). Helium emission can change very rapidly, with clear variability seen from night to night in some BeSS spectra. There is no sdO spectrum detected in Wang et al. (2018). Despite having no frequency groups, this system seems to be able to support a rather strong disk at all observed times.

**TIC 139385056 = FY CMa:** B0.5 IVe. A well-known Be +sdO binary (Peters et al. 2008). The TESS light curve is unusual compared to the rest of the sample. There is low-frequency stochastic variability. The most prominent feature is a pair of signals centered at 8.50 and 8.99 day⁻¹, of roughly equal strength (~2.2 ppt). There are many frequencies between 3 and 9 day⁻¹, which may form some doublets or triplets, and an apparently isolated signal at ~13.4 day⁻¹ (amplitude ~0.5 ppt).

**TIC 207020262 = HD 45828 = BD+08 1366:** There are seemingly five groups in a harmonic series (with the first at 0.75 day⁻¹), but all with comparable amplitudes and all very narrow. Three APOGEE spectra show weak Br11 emission, possibly with shell absorption.

**TIC 21249978 = V868 Ara:** See Section B.4. If this is a classical Be star, then it is possible that the 22.5 day photometric period is related to binarity where the gravitational influence of a companion causes phase-locked structure in the disk at the orbital period (i.e., $m = 2$ density waves; Panoglou et al. 2016, 2018). The triple-peaked $Hα$ profile available on BeSS is qualitatively similar to those of other Be binaries (e.g., HR 2142 and HD 55606; Peters et al. 2016; Chojnowski et al. 2018), and the short photometric period would place this among the shortest-period Be binaries known. This is

speculative, and a more detailed spectroscopic analysis of the system should be undertaken to determine its nature.

**TIC 216158265 = HD 155436:** This is apparently a classical Be star, yet shows only variability consistent with multiple p-mode pulsation and lacks any signals short of 6.5 day⁻¹, but many signals are present between 6.5 and 11 day⁻¹. This object is in the Jaschek & Egret (1982) catalog of Be stars, and is listed as having a spectral type of B0.5IIIn(e) in Garrison et al. (1977). One low-resolution BeSS spectrum from 2012 shows $Hα$ in emission with $E/C \sim 1.8$. While the TESS frequency spectrum is unusual in not having any power at low frequencies, there is nothing to suggest this is not a Be star.

**TIC 216875138 = HD 156172:** The TESS variability of this star is dominated by a pair of signals near 7 day⁻¹, with a clear harmonic of the stronger of the two, and likely multiple combination frequencies. Two low-resolution BeSS spectra show $Hα$ with double-peaked emission, and a high-resolution echelle spectrum shows $Hα$ and $Hβ$ profiles consistent with a typical Be star disk (Appendix B.4 in Labadie-Bartz et al. 2021). There may be weak frequency groups in the TESS data (near 1 and 2 day⁻¹), but since these cannot be reliably distinguished from the red noise floor, we do not report frequency groups for this star.

**TIC 234752466 = HD 150533:** Mainly shows low-frequency stochastic variability and signals consistent with p-modes. There are no lower-frequency periodic signals.

**TIC 234933597 = HD 46484:** B0.5IVe. Very rich frequency spectrum with signals up to 9 day⁻¹. However, there do not appear to be any typical frequency groups. BeSS spectra show a weak and variable disk.

**TIC 279430029 = HD 53048:** B6Vne. Observed in all 13 sectors in cycle 1. The features in the frequency spectrum resemble very narrow groups centered at 1.77 day⁻¹, 3.56 day⁻¹, 5.35 day⁻¹, and 7.14 day⁻¹. These seem to form a typical harmonic series. However, these groups are extremely narrow, and we are only able to resolve them into multiple frequencies thanks to the year-long observing baseline. If this were observed in only one sector, only isolated frequencies would be detected. The wavelet plot does show the groups to be variable in strength over months-long timescales, consistent with beating of very closely spaced modes. Alternatively, these groups could represent isolated frequencies with slowly varying amplitudes, but they are interpreted here as groups.

**TIC 281047621 = ome Ori:** Variable disk, double-peaked, and usually roughly symmetric in BeSS spectra and sometimes shows He wing emission (suggesting a variable but sometimes dense and hot inner disk). It is an unusually simple LC, with a very strong peak at 1.07 day⁻¹, and a harmonic at 2.14 day⁻¹, and lower-amplitude frequency groups at around 0.5 and 1.54 day⁻¹. It looks like it is ionizing some nearby gas. It is difficult to determine if there are groups or just isolated frequencies, but since the amplitude for the main 1.07 day⁻¹ signal clearly (but slowly and mildly) varies, it is assumed to represent a group. It perhaps hosts a weak magnetic field (Neiner et al. 2012c). Interestingly, the rotation frequency derived in Neiner et al. (2003) is 0.765 day⁻¹ ($P = 1.307$ days). This is obviously different from the main TESS signal. However, there is some weak power in TESS at $2 \times 0.765$ day⁻¹.

**TIC 307225534 = V767 Cen:** Variability at the 15% level with a timescale of around 30 days, plus some much lower-amplitude frequency groups at around 0.8, 1.7, 3.5, and
5.2 day⁻¹. Absorption lines seem narrow in BeSS, and Hα emission is single-peaked and about E/C ~2.5 usually, but up to 4. B2IIIPep is from Slettebak (1982), and vsini = 70 km s⁻¹ is from Frémat et al. (2005). The brightness drops by 15% in about 10 days, which is a bit unusual but could be possible for a pole-on system where the inner disk is rapidly dissipating. If this is a classical Be star, and there is no convincing evidence that it is not, then it is at a very low inclination angle due to having a low vsini, single-peaked emission, and narrow absorption lines. Its X-ray flux being hard, thermal, and variable on short timescales is consistent with this being a γ Cas analog (Nazé 2009; Schöller et al. 2017; Nazé & Motch 2018).

TIC 324940394 = HD 148907: B6Ve. There are two frequency groups that are roughly harmonically related, but the first of these is centered at 4.91 day⁻¹, and the second is at ~10 day⁻¹ (i.e., the groups are at unusually high frequencies). There are also isolated frequencies at 6.94 day⁻¹ and 8.95 day⁻¹. There may also be a weaker group centered near 2.5 day⁻¹. This star is located in a somewhat crowded field, and blending is a concern. A more careful analysis is needed to determine if these higher-frequency signals and groups originate in the target star, which would be unusual given the spectral type.

TIC 330435560 = HD 123131: B6.5Ve. The most prominent feature in the frequency spectrum is a narrow group centered at 1.20 day⁻¹, but there are also multiple isolated frequencies at 2.15 day⁻¹, 3.22 day⁻¹, and 4.29 day⁻¹ (2× 2.15 day⁻¹). This system appears to be a very close visual double in high-resolution sky images, and it is possible that the frequency spectrum has contributions from multiple sources.

TIC 338738989 = HD 152060: The amplitude spectrum is dominated by two signals at 9.57 day⁻¹ (amplitude 6.0 ppt) and ~12.19 day⁻¹ (amplitude 5.2 ppt), and two lower-amplitude signals near 10.46 day⁻¹ (1.31 ppt) and 9.25 day⁻¹ (0.65 ppt), and there appears to be a small degree of low-frequency stochastic excess. This is a known β Cephei pulsator (Pigulski & Pojmański 2008). The lack of any apparent lower-frequency pulsation is unusual for a classical Be star.

TIC 400136687 = HD 76568: B1Vnne. Unusual frequency spectrum made up of many narrow groups (or perhaps isolated frequencies) located at 1.28 day⁻¹, 2.28 day⁻¹, 2.91 day⁻¹, 3.56 day⁻¹, 4.07 day⁻¹, and 4.55 day⁻¹.

TIC 401635731 = V1012 Cen: B3Vne. The light curve displays unusual low-frequency variability compared to the majority of the sample. This low-frequency activity may represent short flickers (timescales of ~3 days) that seem to repeat at regular intervals (roughly every 5 days). However, this low-frequency variability appears to be modulated by a beat envelope with a period of about 30 days. If these low-frequency signals do represent flickers, then this system is unusual in that the frequency groups (g1 located at 2.36 day⁻¹, and g2 at 4.80 day⁻¹) are strongest when the amplitude of the low-frequency variability is lowest. KELT data (unpublished) for this system show a fairly strong periodic signal at a period of 5.5 days, which is approximately the timescale where these flicker-like signals repeat in TESS. A single low-resolution BeSS spectrum from 2013 shows Hα in emission with E/C ~3.

TIC 437371531 = V480 Car: B2.5Ve. Both Hipparcos and KELT data show flickers that occur semi-regularly with a period of ~39.22 days (Labadie-Bartz et al. 2017). The TESS data (sectors 9 and 10) seem to confirm this and also show typical frequency groups. This regularity may make this system convenient to study.

TIC 439164152 = CPD-63 2495: O9.5Ve. Has one eclipse near the end of the data set, reaching a depth of about 0.5%. No optical BeSS spectra, but three IUE spectra exist. It is a high-mass X-ray binary (HMXRB; Lutovinov et al. 2013). Chernyakova et al. (2014) gives an interesting study of the object from its 2010–2011 periastron passage. It is known to have an orbit of P<sub>orb</sub> = 1236.9 day, e = 0.87 (Johnston et al. 1992, 1994). It seems like a classical Be star, and is not a supergiant. The eclipse is probably unrelated to the HMXRB and may be a blend or some other source. With an orbital period of about 3.4 yr, an eclipse is extremely unlikely.

TIC 468095832 = 2E 1118.7-6138 = WRAY 15-793: HMXRB and a gamma-ray source, with a spectral type of O9.5III/Ve (Janot-Pacheco et al. 1981). Dramatic variability is seen in KELT (unpublished), with events having amplitudes up to 2 mag, and duration of ~100 days. TESS data show low-frequency trends (including low-frequency stochastic variability) and isolated signals (including at high frequencies).

**Appendix C**

**Table of the Full Original Sample**