title: Short Term Variability of Evolved Massive Stars with TESS II: A New Class of Cool, Pulsating Supergiants
authors: Trevor Z. Dorn-Wallenstein Emily M. Levesque Kathryn F. Neugent James R. A. Davenport Brett M. Morris Keyan Gootkin
text: 
  
  
  
  Trevor Z. Dorn-Wallenstein
  tzdw@uw.edu
  
  0000-0003-3601-3180]Trevor Z. Dorn-Wallenstein
  University of Washington Astronomy Department 
  
  Physics and Astronomy Building, 3910 15th Ave NE  
  
  Seattle, WA 98105, USA 
  
  0000-0003-2184-1581]Emily M. Levesque
  University of Washington Astronomy Department 
  
  Physics and Astronomy Building, 3910 15th Ave NE  
  
  Seattle, WA 98105, USA
  
  0000-0002-5787-138X]Kathryn F. Neugent
  University of Washington Astronomy Department 
  
  Physics and Astronomy Building, 3910 15th Ave NE  
  
  Seattle, WA 98105, USA
  Lowell Observatory, 1400 W Mars Hill Road, Flagstaff, AZ 86001
  
  0000-0002-0637-835X]James R. A. Davenport
  University of Washington Astronomy Department 
  
  Physics and Astronomy Building, 3910 15th Ave NE  
  
  Seattle, WA 98105, USA 
  
  0000-0003-2528-3409]Brett M. Morris
  Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland
  
  
  0000-0003-0922-138X]Keyan Gootkin
  University of Washington Astronomy Department 
  
  Physics and Astronomy Building, 3910 15th Ave NE  
  
  Seattle, WA 98105, USA
  
  
  
  
  Massive stars briefly pass through the yellow supergiant (YSG) phase as they evolve redward across the HR diagram and expand into red supergiants (RSGs). Higher-mass stars pass through the YSG phase again as they evolve blueward after experiencing significant RSG mass loss. These post-RSG objects offer us a tantalizing glimpse into which stars end their lives as RSGs, and why. One telltale sign of a post-RSG object may be an instability to pulsations, depending on the star's interior structure. Here we report the discovery of five YSGs with pulsation periods faster than 1 day, found in a sample of 76 cool supergiants observed by at two-minute cadence. These pulsating YSGs are concentrated in a HR diagram region not previously associated with pulsations; we conclude that this is a genuine new class of pulsating star, Fast Yellow Pulsating Supergiants (FYPS). For each FYPS, we extract frequencies via iterative prewhitening and conduct a time-frequency analysis. One FYPS has an extracted frequency that is split into a triplet, and the amplitude of that peak is modulated on the same timescale as the frequency spacing of the triplet; neither rotation nor binary effects are likely culprits. We discuss the evolutionary status of FYPS and conclude that they are candidate post-RSGs. All stars in our sample also show the same stochastic low-frequency variability (SLFV) found in hot OB stars and attributed to internal gravity waves. Finally, we find four α Cygni variables in our sample, of which three are newly discovered.
  
  
  
  
  
  
  § INTRODUCTION
   
  
  The exact evolutionary pathway a star takes in its late phases and its final fate as a function of initial mass and composition are both incredibly sensitive to the physics of mass loss. Recent advancements have explored the important roles that interior mixing, pulsations, rotation, binary interactions, and magnetic fields can play on the evolution of the most massive stars <cit.>. However, a discrepancy still exists between the highest mass observed red supergiants (RSGs) and the highest mass observed supernova II-P progenitors (known as the red supergiant problem, see ). Many solutions have been proposed by observers and theorists alike, including considering extinction effects to account for underestimated progenitor luminosities (e.g. , ); reexamining bolometric corrections that are extremely sensitive to effective temperature to accurately estimate supernova progenitor luminosities(e.g. ); attemping to quantify biases in progenitor mass estimates (e.g. , ); intricately mapping the landscape of explodability in stellar models <cit.>; and incorporating well-tested prescriptions for RSG mass loss showing that the highest-mass (≥20M_⊙) RSGs may simply evolve back to the blue side of the H-R diagram before explosion <cit.>.
  
  A more direct solution is to find yellow or blue stars that have likely already experienced a RSG phase <cit.>, explicitly determining the mass at which stars no longer end their lives as RSGs. Such post-RSGs allow us to place critical observational constraints on which stars do not simply evolve redward from the main sequence and then explode. Many methods of finding post-RSGs have been attempted. Surface abundance enhancements of CNO-cycle elements are indicative of both envelope loss and convective mixing that extends from the envelope to the core during the RSG phase. Alternately, stars with evidence of past strong mass loss are likely candidate post-RSGs, detected either via infrared excesses caused by warm circumstellar dust, or by direct detection of ejected mass <cit.>. This, however, requires that the CSM is detectable. Finally, though stars' first crossing of the HR diagram proceeds relatively unimpeded, a small number of very luminous yellow supergiants are observed to encounter the “yellow void” where their envelopes become dynamically unstable, resulting in outbursts <cit.>. 
  
  A possibility that has only recently been explored is to search for stable pulsations in evolved massive stars. Many main sequence massive stars are known pulsators <cit.>. Hot stars also exhibit stochastic low-frequency variability (SLFV, manifested as red-noise in the periodogram) that may be attributable to internal gravity waves (IGW; see , with caveats in ). However, beyond lower-mass massive stars found in the upper Cepheid instability strip, κ-mechanism pulsations are not expected as massive stars evolve redward. α Cygni variables are B and A supergiants (blue supergiants, BSGs) that exhibit microvariability and line profile variations, identified as strange-mode pulsations by <cit.>. These strange modes arise in stars with high ratios of luminosity to mass (L/M). <cit.> proposed that α Cygni variables achieve such high L/M values after significant mass loss in the RSG phase. However, while the predicted frequencies roughly correspond to those that are observed, the predicted and observed surface abundances do not match up (though the discrepancy may be resolved by adopting different criteria for convective instability; see ).
  
  While α Cygni variables are promising candidate post-RSG objects, stars evolving from RSGs to BSGs can experience a variety of evolutionary scenarios, including binary interactions or the aforementioned dynamical instabilities, as they evolve blueward. Therefore, a better route may be to search for pulsating yellow supergiants (YSGs) that are closer to the RSG phase, and which may pulsate for the same reason as α Cygni variables. The discovery of a group of pulsators that are cooler than α Cygni variables, and separate from both the “yellow void” and the brightest Cepheids, would thus be of great use. Unfortunately, theoretical modeling of envelope stability in this regime of the HR diagram still encounters convergence difficulties <cit.>.
  
  Thankfully, where theory falls short, observations may provide a path forward. The Transiting Exoplanet Survey Satellite (, ) is collecting lightcurves of the brightest stars across 85% of the sky. While its primary mission has been to search for exoplanets, many of the brightest massive stars have already been observed at two-minute cadence for approximately 27 days at minimum, with observations of stars in the northern and southern continuous viewing zones (CVZs) lasting an entire year. has thus allowed us to measure microvariability in an unprecedentedly large sample of cool supergiants. In <cit.> (Paper I hereafter), we examined a small sample of evolved massive stars that had been observed in Sectors 1 and 2, and found evidence for fast pulsations in three YSGs. Here, we utilize the first 22 Sectors of data, and report the discovery of a group of yellow supergiants that exhibit rapid (<1 d) multiperiodic variability. These stars are more luminous and warmer than the classical Cepheid instability strip, fainter than outbursting yellow “hypergiants,” and notably cooler than the coolest α Cygni variables. We describe our sample and methodology in <ref>, then characterize and discuss the SLFV that is ubiquitous in the sample. We present a new class of fast yellow pulsating supergiants in <ref>, discuss the importance of this new class in <ref> before concluding in <ref>.
  
  
  
  
  § METHODOLOGY
  
  
  
  
   §.§ Sample Selection
  
  
  We first created a sample of cool supergiants with well-measured effective temperature (T_eff) and luminosity (L). <cit.> used spectra obtained with the Hydra multi-object spectrograph on the Cerro Tololo 4-meter telescope to confirm the membership of a large sample of YSGs and RSGs in the Large Magellanic Cloud (LMC), along with updated formulae derived from Kurucz <cit.> and MARCS <cit.> to obtain T_eff and log L/L_⊙ from J-K photometry. Because none of the RSGs published in <cit.> have been observed by (as described below), we also include the Galactic RSGs from <cit.>, <cit.>, and the unique RSG WOH G64 <cit.>. Finally, we discard stars with log L/L_⊙≤ 4 to avoid contamination by lower-mass evolved stars (see, for example, ). 
  
  
  
   §.§ Observations
  
  
  We crossmatched our sample of cool supergiants to the latest version of the Input Catalog (TIC, ) available on the Mikulski Archive for Space Telescopes (MAST). From this sample, we selected all stars with a magnitude in the bandpass fainter than T=4 (where begins to saturate), and brighter than T=12 to obtain sufficient signal-to-noise ratios (SNR) in the lightcurves to detect sub-ppt-level variability (see Paper I). We also omitted any stars with calculated  values above 0.1 to mitigate contamination by nearby stars[See <cit.> for the exact definition of ]. We then downloaded target lists for Sectors 1-22[target lists are available online at <https://tess.mit.edu/observations/target-lists/>], and select the cool supergiants that have been observed at two-minute cadence. This results in a total of 28 YSGs and 48 RSGs. The positions of the YSGs in the HR diagram have a typical error of 0.015 dex in log T_eff and 0.10 dex in log L/L_⊙ respectively.
  
  None of the observed RSGs are in the catalog from <cit.>. For RSGs observed by <cit.>, we used their published estimates of T_eff and log L/L_⊙ (derived from M_bol) where available. The T_eff measurements have typical uncertainties of ±25 K for M stars (T_eff≲ 3810 K, log T_eff≲ 3.581), and ±100 K for K stars. The log L/L_⊙ measurements have typical uncertainties of ∼0.1 dex. One RSG, V772 Cen (= HD 101712), is a known RSG+B star binary; due to this, the derived T_eff and log L/L_⊙ from the TIC are both significantly higher than expected for an RSG. A spectrum of V772 Cen is published in <cit.>, and available on the Vizier online service <cit.>. We obtained this spectrum, and used it to estimate T_eff and log L/L_⊙ for the RSG member of the binary following <cit.>, with comparable uncertainties. For the fourteen remaining RSGs, no suitable archival spectrophotometry exists from which we could estimate T_eff and log L/L_⊙, and so we used the T_eff estimate published in the TIC, as well as the radius measurement to estimate log L/L_⊙. For stars in both the TIC and , the parameters from the TIC show generally good agreement with the results in  to within the errors. However, we do expect the errors on both parameters (especially luminosity) to be significant as the relations used to compute stellar parameters are only validated on dwarfs and giants <cit.>. Table <ref> shows the name, TIC number, coordinates, proper motions, magnitude, log T_eff, and log L/L_⊙ of each star, as well as the source used to determine their position in the HR diagram, whether the star is a RSG or YSG, and if the star is an α Cygni variable or belongs to the newly identified class of pulsating yellow supergiants (see below). 
  
  lccccccccccc
  
  Names, TIC numbers, coordinates, proper motions, magnitudes and positions in the HR diagram of the cool supergiants observed by TESS, ordered by effective temperature from coolest to warmest. The source of the T_eff and log L/L_⊙ measurements is indicated, where N corresponds to <cit.>, L to <cit.>, I to <cit.>, and T to the TIC <cit.>. Typical uncertainties in log T_eff and log L/L_⊙ are 0.015 dex and 0.10 dex respectively in <cit.>. M stars from  have uncertainties of 25 K and 0.1 dex respectively, while the uncertainties in T_eff in K stars are somewhat larger (100 K). Quantities for RSGs derived from the TIC show good agreement with the values published by  where overlap exists. We also indicate whether the star is a RSG or YSG (indicated with “R” or “Y” respectively), and whether the star is a candidate α Cygni variable or belongs to the newly identified class of pulsators.
  Common Name     TIC Number     R.A.     Dec     μ_α     μ_δ      T     log T_eff     log L/L_⊙     Source      RSG/YSG?     Var. Type
  
            [deg]     [deg]     [mas/yr]     [mas/yr]     [mag]     [K]     L_⊙                 
  V1092 Cen     290678703     174.10924243     -61.31944611     -6.709     0.744     5.290     3.534     4.448     T     R     - 
   
   HS Cas     52782147     17.08300080     63.58652909     -2.450     -0.357     5.887     3.535     4.560     T     R     - 
   
  HD 143183     423407817     240.40092730     -54.14322405     -2.301     -3.620     4.245     3.537     5.222     T     R     - 
   
  BD+35  4077     136034302     305.30862034     35.62126593     -2.846     -4.499     5.684     3.556     4.768     L     R     - 
   
   AD Per     348314378     35.12084468     56.99312317     -0.066     -1.423     5.357     3.543     4.587     T     R     - 
   
   KY Cyg     15065085     306.49184826     38.35213201     -3.574     -6.279     4.898     3.544     5.432     L     R     - 
   
  TYC 8626-2180-1     459005094     161.46107152     -59.48870180     -7.080     1.750     4.571     3.547     4.936     L     R     - 
   
   V589 Cas     399355842     26.52283837     60.99352149     -0.952     -0.488     5.843     3.547     4.716     L     R     - 
   
   RS Per     348607532     35.60122973     57.10947226     -0.371     -0.931     5.084     3.550     5.156     L     R     - 
   
   V602 Car     467450857     168.37488668     -60.09134769     -5.425     2.183     4.945     3.550     5.020     L     R     - 
   
   W Per     251118305     42.65788594     56.98341594     0.243     -1.991     5.625     3.550     4.732     L     R     - 
   
   V396 Cen     443405175     199.35433424     -61.58398415     -4.770     -1.758     4.580     3.550     5.212     L     R     - 
   
   BI Cyg     13249363     305.34119647     36.93214587     -2.929     -5.223     4.738     3.553     5.352     L     R     - 
   
   BC Cyg     13325866     305.41061705     37.53303272     -3.856     -5.835     5.094     3.553     5.280     L     R     - 
   
   SU Per     348528265     35.52872734     56.60413801     -0.617     -1.490     4.650     3.553     4.952     L     R     - 
   
   PZ Cas     272324954     356.01366443     61.78949643     -3.110     -1.808     4.972     3.556     5.324     L     R     - 
   
   ST Cep     63963820     337.54474090     57.00085201     -3.517     -2.837     5.150     3.556     4.088     L     R     - 
   
   RW Cyg     15888421     307.21079278     39.98178278     -3.255     -5.511     4.596     3.556     5.156     L     R     - 
   
   TZ Cas     378292562     358.23432055     61.00233067     -3.220     -2.075     5.562     3.556     4.988     L     R     - 
   
   BU Per     264731552     34.72204574     57.42132329     -0.526     -1.106     5.898     3.556     4.764     L     R     - 
   
   V349 Car     457427613     157.39738942     -57.96638247     -7.191     3.632     5.250     3.559     4.808     L     R     - 
   
   V774 Cas     399433806     26.75004525     60.37232574     -1.068     -0.601     5.853     3.559     4.616     L     R     - 
   
  HD  95687     466289471     165.39899669     -61.04883831     -6.746     1.084     4.647     3.559     4.948     L     R     - 
   
   V441 Per     445664243     36.34108308     57.43726049     -0.254     -1.559     5.283     3.559     4.820     L     R     - 
   
  HD 303250     458834083     161.08350153     -58.06484800     -6.875     2.935     5.585     3.559     4.936     L     R     - 
   
   RT Car     458861722     161.19645089     -59.41336782     -7.450     2.914     6.417     3.559     5.260     L     R     - 
   
   V772 Cen     321656644     175.45585098     -63.41457099     -5.508     1.089     5.258     3.560     4.630     I     R     - 
   
  HD 101007     319508664     174.23716722     -61.18277794     -6.647     0.928     4.885     3.562     4.368     T     R     - 
   
   V648 Cas     450147792     42.76645187     57.85553435     -0.184     -1.252     5.912     3.562     4.900     L     R     - 
   
   IX Car     465185147     162.60957843     -59.98238045     -6.054     2.311     4.902     3.562     5.128     L     R     - 
   
   W Cep     65034243     339.11484739     58.42609816     -3.329     -2.132     5.207     3.566     5.466     T     R     - 
   
   V910 Cen     290681168     173.93730736     -61.57806090     -6.699     0.937     5.304     3.568     4.516     L     R     - 
   
   V528 Car     466325776     165.77563786     -60.91072867     -7.130     1.875     4.335     3.568     4.912     L     R     - 
   
   YZ Per     245588987     39.60591607     57.04616613     -0.119     -1.391     5.160     3.568     4.684     L     R     - 
   
   V362 Aur     285640583     81.79257440     29.92105466     -0.678     -2.892     4.886     3.568     4.620     L     R     - 
   
   PR Per     348442493     35.42670692     57.86281915     -0.788     -1.328     5.425     3.570     4.440     T     R     - 
   
   FZ Per     348314886     35.24852231     57.15832387     -0.696     -1.223     5.649     3.571     4.468     T     R     - 
   
   V809 Cas     265186608     349.84905043     62.73977569     -2.257     -2.004     4.117     3.574     4.472     L     R     - 
   
   V439 Per     348671468     35.79610521     57.19943969     -0.308     -0.920     5.803     3.580     4.420     L     R     - 
   
   V605 Cas     348436054     35.09359712     59.67136417     -0.711     -0.953     5.710     3.585     4.920     T     R     - 
   
    41 Gem     337334476     105.06593015     16.07900049     -2.088     -4.853     4.230     3.597     4.341     T     R     - 
   
   RW Cep     422108142     335.77923003     55.96322672     -3.616     -2.349     4.370     3.597     5.470     T     R     - 
   
  HD 155603     188405014     258.61523030     -39.76665102     -0.900     -1.087     4.138     3.601     4.870     T     R     - 
   
   NR Vul     435670188     297.54969991     24.92338263     -2.320     -5.807     5.421     3.602     5.348     L     R     - 
   
   QY Pup     334352580     116.91052662     -15.99068889     -2.162     3.511     4.898     3.608     4.756     T     R     - 
   
  HD  17958     390806332     44.10270614     64.33244354     -3.739     0.017     4.124     3.623     4.548     L     R     - 
   
  HD  33299     367172191     77.64572922     30.79754031     -0.015     -3.142     5.274     3.633     4.044     L     R     - 
   
   AZ Cas     444831689     25.56880634     61.42120644     -2.198     -0.263     7.053     3.656     4.550     T     R     - 
   
  SK -67   57     40603917     77.96699736     -67.16603943     1.552     0.179     11.736     3.656     4.519     N     Y     - 
   
   HV   883     30526897     75.03151958     -68.45001791     1.785     -0.043     11.196     3.680     4.841     N     Y     - 
   
  HD 269953     404850274     85.05069622     -69.66801469     1.718     0.692     9.267     3.692     5.437     N     Y     FYPS 
   
  HD 269110     40404470     77.29420213     -69.60339017     2.081     0.252     10.038     3.750     5.251     N     Y     FYPS 
   
  HD 268687     29984014     72.73273606     -69.43125133     1.833     -0.114     10.465     3.784     5.169     N     Y     FYPS 
   
  HD 269840     277108449     84.04200662     -68.92812902     1.487     0.682     10.132     3.791     5.335     N     Y     FYPS 
   
  HD 269902     277300045     84.53992899     -69.10592146     1.707     0.628     9.790     3.793     5.352     N     Y     FYPS 
   
  HD 269331     179206253     79.50763757     -69.56049032     1.772     0.291     10.114     3.810     5.307     N     Y     - 
   
  RMC 137     404768745     84.65400862     -69.08552151     1.886     0.908     11.878     3.847     4.543     N     Y     - 
   
  CPD-69   430     277172433     84.23672142     -69.27176831     1.778     0.620     11.922     3.857     4.581     N     Y     - 
   
  W61 27-27     277025859     84.01580306     -69.02503035     1.526     0.556     10.770     3.861     4.493     N     Y     - 
   
  HD 269392     179376451     79.96588265     -69.88570140     1.979     0.257     11.961     3.865     4.605     N     Y     - 
   
  HD 269128     40518041     77.59495623     -68.77328288     1.862     0.284     9.189     3.872     5.134     N     Y     - 
   
  HD 269700     425081475     82.96784116     -68.54412683     1.602     0.401     8.808     3.882     5.069     N     Y     - 
   
  HD 270151     389749856     87.26183616     -70.04170277     1.756     0.847     10.561     3.897     4.635     N     Y     - 
   
  CPD-69   491     404852071     85.20335494     -69.28089004     1.801     0.533     10.298     3.914     4.665     N     Y     - 
   
  HD 270754     294872353     71.76854552     -67.11475533     0.754     0.830     11.191     3.915     4.927     N     Y     - 
   
  HD 269655     391810734     82.65769726     -68.41088876     1.495     0.972     11.123     3.924     4.497     N     Y     - 
   
  HD 269997     404933493     85.33497338     -69.08535421     1.658     0.890     9.313     3.927     4.970     N     Y     - 
   
  W61  6-77     389363675     85.56405367     -69.22241953     1.658     0.424     11.069     3.969     4.502     N     Y     - 
   
  HD 269777     276864600     83.57692409     -67.30380846     1.373     0.795     11.153     3.976     5.067     N     Y     - 
   
  CPD-69   394     276936320     83.65036911     -69.76013942     1.642     0.341     10.745     3.984     4.571     N     Y     - 
   
  HD 269992     404967301     85.36531401     -69.80104224     1.969     0.735     9.256     3.990     5.096     N     Y     - 
   
  HD 269786     277022505     83.76501015     -69.75056435     1.872     0.508     9.655     4.000     5.116     N     Y     - 
   
  HD 269101     40343782     77.43830593     -68.76940998     1.760     -0.125     10.577     4.027     4.799     N     Y     α Cyg 
   
  SK -69   68     40515514     77.49501756     -69.11716271     2.311     0.180     11.526     4.029     4.611     N     Y     α Cyg 
   
  HD 268798     30317301     74.28356450     -68.42008185     1.942     -0.074     10.103     4.033     5.071     N     Y     α Cyg 
   
  HD 269769     276936458     83.62856915     -69.78104781     1.783     0.526     10.700     4.037     4.714     N     Y     α Cyg 
   
  
  
  
  Using the Python package , we queried MAST and downloaded all available two-minute cadence lightcurves for each target. The data are provided by the Science Processing Operations Center (SPOC), and include two flux measurements as a function of time: a simple aperture photometry measurement () and flux measurements that have been corrected for systematic trends in the data (). The time at each cadence is the photon arrival time at the solar system barycenter, correcting for the position and movement of the spacecraft. For the following, we used the  lightcurves. To stitch together lightcurves from different sectors, we divided each sector's lightcurve by its median flux. 
  
  
  
  
  
  
   §.§ Stochastic Low Frequency Variability Across the Upper HR Diagram
  
  
  The top panel of Figure <ref> shows the lightcurve of the YSG HD 269953, a star that we concluded is a likely post-RSG in Paper I based on its pulsations and infrared excess measured by Spitzer <cit.>. Individual  measurements are shown as black points[A rapid dimming/brightening event is visible in the Sector 11 lightcurve that arises due to a combination of systematics in the detrending and a discontinuity after the mid-sector downlink. However, as revealed in the wavelet analysis below as well as a by-eye inspection of the periodogram computed only on the light curve before this event, this discontinuity only manifests itself as a low frequency transient and has no effect on the recovered frequencies.]. Dash-dotted vertical lines show the boundaries between sectors, and sector numbers are indicated. The middle panel shows a zoom-in to three different thirty-day portions of the lightcurve, now plotting the data after smoothing with a 128-cadence rolling median. Oscillations can clearly be seen. The bottom panel shows the Lomb-Scargle periodogram of the unsmoothed data <cit.>, calculated with the  package. We use the  normalization option, and divide the power by the number of points in the lightcurve. The resulting quantity is equivalent to the absolute value of the power spectral density, |PSD|, in units of normalized flux squared. We use the default  heuristics for choosing the frequency grid; the maximum frequency is set by the pseudo-Nyquist frequency, f_Ny = 1/(2⟨Δ t⟩) (where ⟨Δ t⟩ is the mean difference in time between two consecutive observations) and the frequency spacing is five times smaller than the Rayleigh frequency, f_R = 1/T where T is the time baseline of the entire lightcurve. The black periodogram is plotted with a linear scaling (corresponding to the black axes labels), and the grey periodogram is plotted with logarithmic scaling on both the frequency and power axes (corresponding to the grey axes labels). Both are scaled to be in units of %/d^-1. As found previously in Paper I, the periodogram displays prominent peaks, superimposed upon a frequency-dependent background. The background shows rising power at low frequency that levels off at the lowest frequencies — i.e., red noise or SLFV – that is clearly visible in log scaling. 
  
  Examining the periodograms of the entire sample shows that SLFV is ubiquitous throughout this region of the HR diagram. SLFV has been identified in hot O and B stars <cit.>, and its presence throughout this sample of A-M supergiants suggests that it is in fact a ubiquitous feature of massive stars. Figure <ref> shows the periodograms of four stars. The power is normalized to have a maximum value of 1, and an arbitrary offset constant is added for clarity. The top two periodograms are calculated for two “normal” supergiants that are representative of the overall sample: the red supergiant BD+35 4077, and the yellow supergiant HD 270754. Their periodograms are dominated by SLFV, and they display no strong peaks. The bottom two periodograms belong to two yellow supergiants: HD 269101 and HD 268687. While both stars' variability are dominated by SLFV, they also show visible peaks in their periodograms.
  
  
  
  Searching for significant peaks in a periodogram where the background is independent of frequency (i.e., white noise) is a well-established problem, and peaks can be selected based on various signal-to-noise estimates, or by determining the false alarm probability (FAP), e.g., the probability that a peak with a given power could randomly arise given the null hypothesis of white noise[Note that this is not the probability that the peak is not real. In symbols, P( data| noise) ≠ P( noise| data)] <cit.>. Methods for calculating FAPs in the case of red noise are often incredibly computationally expensive. Fast alternatives have been proposed <cit.>, however these methods are highly model-dependent, and choosing the wrong noise model can dramatically affect the FAP estimate. Currently, no physical theory for SLFV has been uniformly agreed upon by the community, and so we refrain from utilizing these methods to calculate the FAP.
  
  Instead, we follow <cit.> and <cit.>, and use the  routine within the  package to fit the amplitude spectrum (α(f), obtained by taking the square root of the PSD) with a phenomenological model. We adopt the function
  
      α(f) = α_0/1+(2πτ f)^γ+α_w
  
  from <cit.>, where f is the frequency, α_0 is the amplitude as f→0 in units of normalized flux, τ is a characteristic timescale in days on which the noise is correlated, γ sets the slope of the red noise, and α_w is an additional parameter we add in to model the white noise floor at the highest frequencies, also in units of normalized flux. We note that Equation (<ref>) is equivalent to the function adopted by <cit.>, and the characteristic frequency in those works is equivalent to ν_char=(2πτ)^-1. We fit the base-10 logarithm of the amplitude spectrum, calculated as the square root of the PSD, to avoid artificial weighting of real peaks at high frequencies. Note that we do not first prewhiten the coherent variability discussed below from the lightcurves. However, compared to the periodograms calculated from, e.g., the CoRoT lightcurves of hot stars studied by <cit.>, the power of the low frequency excess seen in this sample is far stronger than the power of the observed peaks, with the exception of HD 269953, and so we don't expect the fit parameters to be significantly affected. In the case of HD 269953, fitting the logarithm of the amplitude spectrum mitigates any significant effect. To test this, we performed the iterative prewhitening procedure described below, and recorded the values of α_0, τ, and γ at each step. We found that, after prewhitening the four significant frequencies found below, α_0 changed by a factor of 1.4,  τ by a factor of 1.7, and γ by a fact of 0.95. As HD 269953 is the worst case scenario, we decided that these changes do not change the results shown in Figure <ref>.
  
  lccccc
  
  Names, TIC numbers, fit parameters from Equation (<ref>), and corresponding errors for all stars in our sample.
  Common Name     TIC Number     α_0     τ     γ     α_w 
  
             ppt     d          ppt  
  V1092 Cen     290678703     194.9581 ± 17.3575     0.3464 ± 0.0318     1.755 ± 0.049     7.3583 ± 0.0077 
   
   HS Cas     52782147     414.3034 ± 33.9049     0.1181 ± 0.0095     1.784 ± 0.042     11.6105 ± 0.0196 
   
  HD 143183     423407817     114.8980 ± 12.2582     0.2139 ± 0.0267     1.601 ± 0.061     7.9003 ± 0.0111 
   
  BD+35  4077     136034302     198.9679 ± 25.9046     0.6764 ± 0.0985     1.544 ± 0.056     8.5699 ± 0.0107 
   
   AD Per     348314378     173.5579 ± 10.9315     0.0597 ± 0.0039     2.268 ± 0.092     18.0802 ± 0.0291 
   
   KY Cyg     15065085     286.2191 ± 37.0462     0.6387 ± 0.0998     1.359 ± 0.042     13.1549 ± 0.0172 
   
  TYC 8626-2180-1     459005094     546.2201 ± 69.9013     0.3636 ± 0.0479     1.456 ± 0.032     9.0229 ± 0.0141 
   
   V589 Cas     399355842     92.4473 ± 7.3775     0.0915 ± 0.0080     2.142 ± 0.113     12.1910 ± 0.0188 
   
   RS Per     348607532     383.0987 ± 28.7313     0.0859 ± 0.0069     1.654 ± 0.037     14.4608 ± 0.0264 
   
   V602 Car     467450857     163.9396 ± 10.5384     0.1408 ± 0.0098     2.165 ± 0.093     21.4533 ± 0.0220 
   
  This table is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.
  
  
  
  
  Figure <ref> shows the square of Equation (<ref>), |α(f)|^2, with α_0=1, and α_w=10^-3. The top (bottom) panel shows the effects of varying τ (γ) at constant γ (τ). Figure <ref> shows a summary of our fits to all stars in our sample. The top row shows a histogram of the best-fit values of each of the four fit parameters. We also plot these parameters as a function of log T_eff in the bottom row. Errorbars are calculated from the covariance matrix returned by , and the color corresponds to increasing log L/L_⊙ from darker to lighter colors.
  
  
  
  Examining the bottom panels of Figure <ref> from left to right: α_0, the amplitude of the SLFV, is suppressed at temperatures around log T_eff∼3.7-3.8. τ slightly increases with increasing temperature. While we don't place any bounds on τ when performing the fits, there is an implicit upper limit to τ that can be seen in the data, as the SPOC processing pipeline and our sector-combining procedure effectively erases correlations in the lightcurves on timescales longer than a few days. The slope of the noise, γ, clusters around γ = 2 for the RSGs, while favoring slightly smaller values for YSGs. Finally, the white noise, α_w is systematically higher for most RSGs. This is not surprising, as all the RSGs in our sample are Galactic, and most were only observed for a single sector. All YSGs are in the LMC, and thus were observed nearly continuously during the first year of observations, with the exception of small gaps when e.g. individual stars passed through the gaps between CCDs in camera 4. Therefore, we expect the signal in the periodogram (relative to the intrumental noise) to be much higher for these stars.
  
  
  
    §.§.§ Convection or Something More?
  
  
  As demonstrated above, the characteristic timescale of the background noise, τ, appears to increase with increasing temperatures. Given that τ roughly corresponds with the characteristic timescale over which the stochastic variability is correlated, this gives us some clues as to the origin of the noise.  Much lower mass yellow stars like our Sun show similar low-frequency power excesses due to granulation in their outer layers. <cit.> demonstrate that, from first principles, the characteristic convective timescale, τ_conv (denoted τ_eff in ) scales with the surface gravity, g, and effective temperature as g^-0.85T_eff^-0.4. Because g∝ M R^-2 and R^-2∝ T_eff^4 L^-1, this implies τ_conv∝ T_eff^-3.4(L/M)^0.85. Thus as a massive star evolves rightward at essentially constant L in the HR diagram, the convective timescale increases strongly as a function of decreasing temperature (with a small boost as the star's L/M increases as the star loses mass). This is the exact opposite of the trend we observe, implying that the low-frequency variability that we see in our sample is not (solely) a result of surface convection, at least in the warmer stars. 
  
  In much hotter massive stars, correlated stochastic variability has been linked with sub-surface convection zones <cit.>, that may interact with pulsations <cit.>. More recently, it was suggested by <cit.> to be a sign of internal gravity waves (IGWs) arising from the boundary of internal convective and radiative layers <cit.>. Perhaps the stochastic variability in this sample is connected with that seen in hot stars? Recently, <cit.> characterized stochastic variability in a sample of 70 OBA stars spanning a range of temperatures above ∼10^4 K and masses between ∼5 and 80 M_⊙ (see Figure 2 in that work). Unfortunately, their sample has very few post-main sequence stars, especially in the mass range of our sample, and so we are unable to construct a complete sequence of α_0, γ, or τ≈ 1/ν_char as stars evolve across the HR diagram. If stochastic variability is attributable to IGWs, such an evolutionary sequence would be an incredibly powerful means of applying asteroseismology to massive stars as they near the ends of their lives.
  
  It is also possible that sub-surface processes that aren't IGWs are causing the low frequency stochastic variability. Therefore, we cannot uniquely identify the stochastic variability with any particular source, and reserve such identification for a study of low frequency stochastic variability in massive stars across the entire HR diagram. That said, the RSGs display power law slopes clustered around γ=2, consistent with what was found by <cit.> in AAVSO data — though in a significantly higher frequency range, and with lower amplitude than the lightcurves studied by  — and attributed by those authors to convective processes. Indeed, though the observed scaling of τ with T_ eff is inconsistent with surface convection, the timescales of simulated turbulent convection in stellar interiors (albeit in lower mass stars) do not show this simple scaling <cit.>.
  
  However, <cit.> found that there is no point at which the power spectra of RSGs turn over; as observing time increases, more power at low frequencies is recovered. Because the observations of RSGs only span a few sectors at most, and correlations on long timescales are smoothed out by the detrending performed by the SPOC, we might expect to see values of τ clustered at the maximum value possible given the detrending (a few tens of days). However, RSGs display the smallest observed values of τ seen in our sample, suggesting that, at the precision of , RSG microvariability is not entirely consumed by convective noise. One possible way to distinguish between subsurface convection and IGWs from core convection may be to measure the macroturbulent velocity of these stars, and compare these measurements with the stars' locations in the HR diagram and observed values of α_0, τ, and γ (see both  and ). It is also entirely possible that both scenarios are at play, and contributing to the observed stochastic variability in this sample.  Regardless of the physical origin, stochastic variability is ubiquitous across the upper HR diagram, from hot stars <cit.> to cool supergiants (this work and ). Characterizing this variability and determining its origin has the potential to offer critical insight into the evolution of massive star interiors from birth to death.
  
  
  
  
  
  
  § FAST YELLOW PULSATING SUPERGIANTS
  
  
  After fitting the amplitude spectrum of each star to obtain the best-fit model, α̂(f), we divide its signal out of the power spectrum by computing |PSD| / α̂(f)^2. The resulting quantity has no formal definition, but is incredibly useful at showing the power of peaks relative to the background. Hereafter, we refer to the background-normalized power spectrum as the residual power spectrum (RPS). Figure <ref> shows the HR diagram, where each star in the sample is replaced by its RPS between 0 and 5 d^-1, normalized by its maximum value and scaled to fit in the plot. Note that because of this scaling, the relative heights of peaks in two different RPS have no relation, but the relative heights of two peaks in the same RPS are meaningful. In particular, plotting the RPS in this way allows us to simultaneously assess the approximate signal to noise in the periodogram as a function of each star's location in the HR diagram. A subset of the nonrotating, Z=0.006 evolutionary tracks calculated with MESA and described below are plotted as thin black lines with the initial masses.
  
  
  
  For the majority of stars in this sample, their lightcurves appear to be entirely composed of SLFV, and their RPS show either no peaks, or small peaks with low signal to noise. There are, however, two groupings of stars with high signal-to-noise peaks in their RPS. The first group, comprised of a vertical strip of four stars with log T_eff≈ 4.0 lies in the region of the HR diagram where α Cyg variables are expected to be found[The frequencies found in the lightcurves of these stars are somewhat higher than in other α Cyg variables <cit.>. However, most observations of α Cyg variables have been taken from the ground, where detecting frequencies around ∼1 d^-1 is difficult. While not the focus of our work, we note these candidate α Cyg variables for use by other authors.], which we highlight in blue. We searched for any past work that has studied their variability, as listed in SIMBAD <cit.>, and find the following:
  
  
      
    * HD 268798 was previously identified as an eclipsing binary by <cit.>, and as a rotational variable with ellipsoidal variations by <cit.>. It has not previously been identified as an α Cyg variable.
      
    * HD 269101 was identified as a candidate α Cyg variable by <cit.>, and incorrectly identified as a Slowly Pulsating B-star by <cit.> (likely due to its entry in SIMBAD erroneously listing it as an early-B supergiant).
      
    * HD 269769 has not previously been studied in the time domain, and thus has not previously been identified as an α Cyg variable.
      
    * Sk-69^∘ 68 has not previously been studied in the time domain, and thus has not previously been identified as an α Cyg variable.
  
  
  Apart from these candidate α Cyg variables, we also find a cluster of five stars, all with 5.1 ≤log L/L⊙≤ 5.5 and 3.69 ≤log T_eff≤ 3.8 with high signal to noise peaks in their RPS. We highlight these stars in green in Figure 6. This region of the HR diagram contains no other stars[One star, HD 269331, has a similar luminosity and a log T_eff that is 0.02 dex higher than the warmest identified pulsator. Its RPS is low signal-to-noise and shows no significant peaks. In Paper I we identified two prominent bumps in the first two sectors of data. Examination of the remaining sectors shows that these bumps are present throughout the lightcurve, and that the star is displays variability with an amplitude of ∼1 ppt. With no significant RPS peaks, we exclude HD 269331 from our subsequent analyses; however it is possible that HD 269331 is a genuine member of this novel class of supergiant pulsator.]. We list the common names (found on SIMBAD), TIC numbers, coordinates, temperatures, and luminosities for the five stars in Table <ref>. They are well-separated in the HR diagram from the lower-right edge of the yellow void — a region of the HR diagram occupied by a very small number of stars that exhibit extreme variability and mass loss due to dynamical instabilities in their atmospheres <cit.> — which is shown in goldenrod, and the upper-left edge of the Cepheid instability strip (shown in cross-hatched orange, derived from Z=0.006, nonrotating stellar models on their first crossing of the HR diagram; see ). This group includes the three pulsating YSGs previously identified in Paper I, as well as two newly identified stars. Many of the frequencies found in their lightcurves (see below) are on timescales shorter than 1 day, and, as discussed in Paper I, are hard to explain with rotational or orbital effects given the large radii of YSGs. Furthermore, while it is possible that these frequencies may arise in the winds of these stars, we deem it unlikely that only YSGs in this region of the HR diagram would show coherent modulations in their winds.[Regardless, with the notable exception of a small number of incredibly luminous YSGs that have undergone outbursts, YSG winds remain poorly understood.] Finally, spectra of all five stars from <cit.> indicate that they are all fairly typical YSGs, though HD 269902 has a slightly weaker Ca2 triplet. Therefore, we adopt the name “Fast Yellow Pulsating Supergiants” (FYPS)[The authors acknowledge the poor adjective ordering in this acronym. However, we believe that FYPS is easier to pronounce than FPYS.] for these stars, and discuss them below.
  
  llllcc
  
  Names, TIC numbers, coordinates, temperatures, and luminosities of the five pulsating YSGs. 
  Common Name     TIC #     R.A.     Dec.     log T_eff     log L/L_⊙ 
  
            ^∘     ^∘     K     L_⊙ 
  HD 269953    	404850274    	85.050696    	-69.668015     3.692    	5.437 
  
  HD 269110    	40404470    	77.294202    	-69.603390     3.750    	5.251 
  
  HD 268687    	29984014    	72.732736    	-69.431251     3.784    	5.169 
  
  HD 269840    	277108449    	84.042007    	-68.928129     3.791    	5.335 
  
  HD 269902    	277300045    	84.539929    	-69.105921     3.793    	5.352 
  
  
  
  
  
  
   §.§ Chance or New Class?
  
  
  All five of the FYPS are located in the LMC. Each pixel is 21” on a side (∼17 ly at the distance of the LMC). Furthermore, YSGs are found in crowded regions with many hot young stars, making it highly unlikely that the starlight in the optimal aperture defined by the SPOC is coming only from these stars. This effect is somewhat mitigated by 's relatively red passband (centered at 7865 Å); while in bluer passbands, the flux in the aperture may contain significant flux from nearby O and B stars, cool, evolved evolutionary phases of massive stars that dominate the flux in the aperture are significantly rarer due to their shorter lifetimes. Additionally, the binary fraction of massive stars is high <cit.>; even the most evolved red supergiants that are the most likely to have interacted and merged with a companion have a binary fraction of ∼20% <cit.>. Therefore it is possible that we are recovering five stars with pulsating companions that happen to be located in a small region of the HR diagram by chance. Assuming the initial massive star binary fraction (as well as the initial period and mass ratio distributions) is roughly constant for the masses of stars in our sample, it is effectively equally likely for any star in our sample to have a pulsating companion (with decreasing likelihood for the largest stars in our sample). Finally, detrending of time series photometry can generate spurious low frequency peaks, which propagate to higher frequencies via harmonics and combinations with real peaks <cit.>. Potentially, it could appear as if these five stars, which may be ordinary YSGs, are pulsating when they are in fact not.
  
  
  Therefore, we need to assess the likelihood that the five stars identified above are otherwise normal YSGs, and their lightcurves are all contaminated by starlight from actual pulsating stars (whether from nearby stars in the aperture or a binary companion), or contain spurious periodic signals introduced by detrending. If five stars were randomly selected as “pulsators” from our sample due to contamination, we would not expect them to be found in such a small region of the HR diagram. To determine the extent to which crowding may have influenced this detection, we can ask the question: assuming any of the stars in our sample could have randomly been contaminated by pulsators, how likely is it to draw five stars from our sample and have them form a grouping in the HR diagram with equal or lesser size. We can answer this question directly with a simple bootstrap analysis. For each of the 18,474,840 unique subsamples of 5 stars, we calculate the dimensions, Δlog T_eff and Δlog L/L_⊙, of the smallest box in the HR diagram that contains the subsample. Note that this analysis does not account for the decreasing likelihood of finding a binary companion around a large supergiant, and so in just the case of contamination by a binary companion, the derived likelihood is an uppper limit. Figure <ref> shows the two-dimensional histogram of Δlog T_eff and Δlog L/L_⊙. The actual values of Δlog T_eff and Δlog L/L_⊙ calculated from the five FYPS is indicated by the red star. We find that only 8.8% (1.2%) of all possible subsamples have equal or lesser ranges of Δlog T_eff (Δlog L/L_⊙). All told, only 0.07% of all possible subsamples are bounded by a smaller region in the HR diagram. If we repeat this calculation with only the LMC YSGs, this number decreases to 0.03%. Thus, we deem it exceedingly unlikely that the lightcurves of these five stars (and only these five stars) happened to have randomly been contaminated by a nearby pulsator or pulsating companion, and conclude that we have discovered a genuine new class of pulsating star.
  
  
  
  
  
  
   §.§ Variability and Pulsation Frequencies of FYPS
  
  
  
  
    §.§.§ Prewhitening
  
  
  The extraction and measurement of individual pulsation frequencies from lightcurves with frequency-independent noise is a fairly well-defined procedure: the strongest peak can be identified based on its amplitude or FAP, and a prewhitening procedure can be applied to iteratively fit and subtract sinusoids from the lightcurve corresponding to the extracted frequencies until some noise threshold is reached <cit.>. The resulting lists of frequencies, amplitudes, and phases in principle, completely describe the coherent variability found in the lightcurve. However, in the case of frequency dependent noise (in this case, stochastic low frequency variability of astrophysical origin), spurious peaks that are random fluctuations superimposed on the noise are extracted, while true peaks that lie on top of the much lower-amplitude white noise at higher frequencies can be ignored. To account for this effect, we adopt the procedure used by <cit.>, with the following modifications and stopping criterion:
  
  
      
    * At each stage of prewhitening, we fit the (amplitude) spectrum as described above, and obtain the RPS before selecting a frequency to prewhiten. We note that this is the opposite procedure adopted by <cit.>, who prewhiten coherent frequencies before fitting the amplitude spectra to characterize the stochastic background. This is because the power of the stochastic variability is much higher than that of the peaks in the periodogram, especially at low frequency. However, because we are not removing the stochastic variability from the lightcurve itself, this only helps us locate the peaks in the RPS. The frequencies, amplitudes, and phases we obtain are otherwise identical to what we would obtain following <cit.>.
      
    * Paper I identified multiple harmonics of some recovered frequencies. To properly treat potential harmonics, at each stage of prewhitening, we fit both the selected frequency and the amplitude and phase of its first two harmonics. There are some instances where we select a frequency that is itself a harmonic of another lower amplitude frequency, and so the fundamental is not removed by the harmonic fit. We note all instances when this occurs below.
      
    * In addition to saving the best-fit parameters of each sinusoid, we calculate the associated errors on each parameter, using the formulae given in <cit.> and <cit.>:
      
      ϵ(f_j)    = √(6/N)1/π Tσ_j/A_j
      ϵ(A_j)    = √(2/N)σ_j 
      ϵ(ϕ_j)    = √(2/N)σ_j/A_j
  
      where N is the number of points in the lightcurve, T is the time baseline of the lightcurve, f_j, A_j, and ϕ_j are the frequency, amplitude, and phase extracted at the j^ th prewhitening stage, and σ_j is the standard deviation of the flux at the same stage. We also record the value of the RPS at the selected frequency, and the SNR, calculated as the peak height divided by the standard deviation of the RPS in a narrow window between f_max±2 f_R and f_max±7 f_R.
      
    * As a stopping criterion, we proceed until we reach a minimum in the Bayesian Information Content (BIC, ) of the fit:
      
      BIC    = -2ln(ℒ) + mln(N) 
      
              -2ln(ℒ)    = ∑_i=1^N(y_i - F(t_i,Θ_m))^2/σ_i^2
  
      where ℒ is the likelihood (to within a constant), m is the number of free parameters in the fit (7j+7 at the j^ th stage of prewhitening, beginning with j=0), y_i are the original fluxes observed at times t_i, F(t_i,Θ_m) is the sum of all of the fit sinusoids evaluated at fit parameters Θ_m, and σ_i are the normalized errors in the original lightcurve <cit.>.
  
  
  As a postprocessing step, we discard frequencies that are quite similar to each other (i.e., the difference in frequencies is within 1.5 f_R, keeping the earliest frequency found). These similar and spurious frequencies can arise due to the short length of the observing baseline <cit.>. The unique frequencies, amplitudes, and phases, corresponding formal errors, RPS peak heights, and RPS SNRs found for each star are listed in Appendix <ref>. The frequencies extracted from the FYPS are all between ∼0.5 and ∼4.6 d^-1 with semi-amplitudes ranging between ∼40 and ∼280 ppm.
  
  From the final list of extracted frequencies, we search for harmonics of the form f_j/f_i = n, where n is an integer greater than 1, that satisfy
  
      nf_i - f_j ≤√((nϵ(f_i))^2 + (ϵ(f_j))^2)
  
   i.e., f_j is an exact integer multiple of f_i to within the errors, and the k^ th harmonic corresponds to n = k+1 (e.g., the first harmonic is n=2). We also search for frequency combinations in the form f_i+f_j=f_k, such that 
  
      f_1 + f_j - f_k ≤√((ϵ(f_i))^2+(ϵ(f_j))^2+(ϵ(f_k))^2).
  
  
  
  
  
    §.§.§ Wavelet Analysis
  
  
  In addition to the frequencies extracted from the entire lightcurve, we also attempt to determine whether the frequencies and amplitudes are stable. To that end, we employ a time-frequency analysis to search for variability in the dominant frequencies. We calculate the Weighted Wavelet Z-transform <cit.>, an extension of wavelet analysis with the Morlet wavelet:
  
      ψ(t,τ,ω) = e^iω(t-τ)-cω^2(t-τ)^2
  
  where here τ is the time of the center of the wavelet, ω = 2π f is the angular frequency, and c sets the width of the Gaussian envelope and is chosen to be sufficiently small so that the wavelet decays appreciably over the course of one cycle. Here we adopt c=0.0125, following <cit.>; in principle, smaller or larger values of c can be chosen to alter the time and frequency resolution, which is frequency dependent. The discrete wavelet transform can be converted into a projection onto the continuous basis functions
  
      Φ_1(t,ω,τ)    = 1 
      Φ_2(t,ω,τ)    = cosω(t-τ)
      Φ_3(t,ω,τ)    = sinω(t-τ)
  
  and we now fold the Gaussian envelope into a weighting function for each data point at time t_i that depends on the frequency and time center of the wavelet:
  
      w_i(τ,ω) = e^-cω^2(t_i-τ)^2
  
  This change ensures that a wavelet centered on a gap in the data won't pick up small amplitude random fluctuations on either side of the gap, thereby suppressing false power that can often arise in wavelet transformations of unevenly sampled data <cit.>. At each τ and ω, we calculate the number of effective data points,
  
      N_eff(τ,ω) = (∑_i w_i(τ,ω))^2/(∑_i w_i^2(τ,ω))
  
  and the weighted variance of the flux (which we denote x):
  
      V_x(τ,ω) = ⟨ x|x⟩ - ⟨Φ_1 | x⟩^2
  
  where the inner product of two functions, ⟨ f|g⟩ is defined as
  
      ⟨ f|g ⟩ = ∑_iw_i(τ,ω) f(t_i)g(t_i)/∑_jw_j(τ,ω)
  
  We also calculate the weighted variance of a sinusoidal fit to the model:
  
      V_y(τ,ω) = ⟨ y|y⟩ - ⟨Φ_1 | y⟩^2
  
  defining
  
      y = y⃗_Φ·Φ⃗
  
  where Φ⃗ is a vector containing the basis functions, Φ. y⃗_Φ is a vector containing the coefficients of the projection onto the basis functions,
  
      y⃗_Φ = S^-1y⃗_b
  
  with the entries of the matrix S equal to S_ab=⟨Φ_a | Φ_b⟩ and entries in the vector y⃗_b are y_b = ⟨Φ_b | x ⟩. With these ingredients, we can calculate the WWZ at each τ and ω as
  
      WWZ = (N_eff - 3) V_y/2(V_x-V_y)
  
  Finally, we set WWZ=0 if min(t_i-τ) > 2π/ω (i.e., the nearest data point is more than one cycle away from the center of the wavelet), to reduce computation cost.
  
  
  
    §.§.§ HD 269953
  
  
  
  
  HD 269953 is the brightest FYPS discovered with log L/L_⊙ = 5.437. We initially analyzed it in Paper I, and suggested that it was a post-RSG, largely due to its infrared excess hinting at past mass loss. It was previously studied by <cit.>, who noted its variability. Figure <ref> shows the RPS of HD 269953, with the four unique frequencies recovered by prewhitening in grey. After searching for harmonics, we find one instance where the fundamental at 1.335 d^-1 has a lower amplitude than the first harmonic at 2.671 d^-1 and so both frequencies are recovered. No combination frequencies are recovered, implying the presence of three independent frequencies in the lightcurve of HD 269953. The presence of these non-aliased frequencies may indicate the presence of multiple oscillation modes in HD 269953. Alternately, the aperture may contain two different pulsating stars. Without a better model of YSG pulsations, we cannot conclusively determine which scenario is more likely; however, if the pulsation mechanism of FYPS is identical to that of α Cygni variables, the presence of multiple modes may not be surprising <cit.>.
  
  Figure <ref> shows the lightcuve (top, black), rolling median of the flux (top, green), RPS (right), and WWZ (center) calculated on 500 linearly spaced time points, and 1000 frequency points on a log_2 grid between log_2(2π f) = -1 (f≈0.08 d^-1) and 5.5  (f≈7.2 d^-1). Frequencies extracted via prewhitening are shown as horizontal white (WWZ) or grey (RPS) lines. At frequencies below 1 d^-1, the WWZ shows transient events that are associated with times where systematics in the data appear to be present — e.g. the discontinuity after the mid-sector downlink at Time≈1620 days. However, the wavelet map demonstrates that these transient events have no effect on the highest amplitude frequencies in the RPS. At higher frequencies, the frequency of maximum power in the WWZ rapidly changes as a function of time, appearing to oscillate between the three lowest extracted frequencies. There even appear to be times when the peak in the WWZ almost disappears (e.g., around Time=1400; note the gap at Time=1500 is due to a gap in the data). This indicates that the detected pulsations are not stable, with amplitudes changing on timescales of days. Perhaps the modes are stochastically excited and damped on these timescales. YSG interiors are complex, with multiple boundaries between convective and radiative zones, which may be responsible for driving the pulsations. However, with no reliable models of YSG pulsations, we can only speculate at this time. Of the seven higher frequency peaks detected, five are harmonics of lower-frequency signals as discussed. Some brief, low amplitude transient events are associated with the peaks in the RPS corresponding to harmonics of lower frequency signals. Unfortunately, a drawback of the WWZ (and most time-frequency analyses in general) is that potential interesting high-frequency features are smeared out in exchange for increased time resolution.
  
  
  
  
  
    §.§.§ HD 269110
  
  
  HD 269110 was also discussed in Paper I, and had the lowest frequency signals detected there. It has not been analyzed by any other modern variability studies. With our updated prewhitening scheme, we extract four unique frequencies, including three incredibly close to each other: a main peak with the highest SNR at f = 0.553 d^-1, and two small peaks, each separated from the main peak by Δ f = ± 0.011 d^-1 (∼10 times the resolution of the periodogram). The peak in the RPS corresponding to the first harmonic of this frequency shows similar structure, though the higher frequency subpeak is very low signal to noise. The lightcurve, RPS, and WWZ are shown in Figure <ref>. Only the three closely spaced peaks can be seen in the WWZ, and those peaks appear to fade and reappear semi-regularly. While the modulation in the WWZ is not sinusoidal, extracting the WWZ in a 0.2 d^-1 band around the closest frequency to the main peak and calculating the power spectrum reveals a strong peak at f = 0.010 d^-1, quite similar to the frequency difference between the three peaks recovered by prewhitening. One possibility is that HD 269110 is a binary system. Pulsating stars in binaries exhibit frequency modulation similar to what we see <cit.>, and eccentric close binaries can induce modulations in pulsation amplitudes on the orbital timescale <cit.>, also in line with what we have detected. An alternate hypothesis is that the f = 0.553 d^-1 is split by rotational effects <cit.>, and 0.011 d^-1 is the rotational frequency. 
  
  If HD 269110 is a binary system with a 20 M_⊙ YSG primary and a 1/0.011 ≈ 91 d orbital period, a companion star would have to be ∼100 M_⊙ in order for the semimajor axis of the orbit to be larger than the stellar radius derived from the parameters listed in Table <ref>. Such an object would have to be a black hole in order to not be significantly brighter than the primary, in which case it would still be more massive than any known stellar mass black hole, a scenario we deem to be incredibly unlikely. Alternately, rotational modulation of the WWZ on a ∼90 d timescale requires invoking a misalignment of the pulsational and rotational axes. Furthermore, a 91 day rotation period is either quite close to, or exceeds the critical rotation periods of YSGs in the Geneva models, depending on the mass. Regardless, the similarity between Δ f and the characteristic timescale extracted from the WWZ is intriguing, and warrants follow-up observations. We discuss the implications on the evolutionary status of FYPS below. 
  
  
  
  
  
  
    §.§.§ HD 268687, HD 269840, & HD 269902
  
  
  
  
  HD 268687 was the third candidate pulsating YSG found in Paper I. However, the most prominent frequencies we had found previously were superimposed upon the part of the power spectrum dominated by low frequency stochastic variability. Apart from Paper I, it was identified as a variable star in OGLE photometry <cit.>. With prewhitening, we extract six unique frequencies. The RPS and extracted frequencies are shown in the left panel of Figure <ref>. The RPS displays a broad comb of peaks, all but one of which correspond to periods faster than one day. None of the low-frequency peaks identified in Paper I are recovered; however, the peaks previously identified as harmonics of the (now nondetected) dominant peak are recovered. From this, we conclude that the dominant peak identified in early data without correcting for the SLFV may have belonged to the overall pattern of peaks, even though it is no longer detected. Searching for harmonics in the recovered frequencies reveals one frequency at 1.844 d^-1 with a detected first harmonics at 3.687 d^-1. A search for combination frequencies shows that the sums of the lowest frequency peak and both the fundamental and first harmonic are also recovered frequencies. Interestingly, the four recovered frequencies in ascending order starting with 1.844 d^-1 are roughly equally spaced, with a frequency difference of Δ f = 0.922 d^-1.
  
  Pulsations in both HD 269840 and HD 269902 are newly discovered; the former was identified as a variable in OGLE photometry <cit.>. Similar to HD 268687, the RPS of both stars show a broad frequency comb (center and right panels of Figure <ref>) with a total of seven and three unique frequencies extracted from HD 269840 and HD 269902 respectively. We do not find any remaining harmonics or combination frequencies amongst the frequencies extracted from both stars. As in HD 268687, the four frequencies recovered in HD 269840 above 2.270 d^-1 are regularly spaced with a spacing Δ f ≈ 0.567 d^-1 (though the spacing between the highest two frequencies is 0.001 d^-1 higher). No such regular spacing is found in HD 269902.
  
  Motivated by the apparent regular spacing of peaks in the RPS of all three stars, we searched for evidence that each of the three FYPS exhibits a harmonic chain of peaks. In all three panels of Figure <ref>, we assume that one of the extracted frequencies is a fundamental frequency and plot its first few harmonics. In the case of HD 269840, we instead assume that the frequency of the tallest peak at f=2.270 d^-1 is four times the frequency of the fundamental, which lines up with a peak in the RPS that is not selected by our prewhitening procedure. In all three cases, we find an exceptional match between the assumed harmonic chain and most of the peaks in the RPS. However, the actual observed frequencies are inconsistent with regular harmonic patterns to within the errors
  
  This fact could be due to one of two causes. The offsets from an even spacing pattern could be caused by, e.g., structural glitches <cit.>. Such glitches can be used to assess sharp features in the stellar structure that would otherwise be inaccessible by other means. Alternately, we may have underestimated the errors for the extracted frequency. Instrumental correlations exist in data taken by the Kepler and CoRoT missions, and can add to the uncertainty in extracted frequencies <cit.>. As the theory of pulsations in YSGs is still nascent, we have no asteroseismic model for these stars, and thus cannot determine whether the observed offset from an even spacing pattern is astrophysical. However, we can determine the extent to which the uncertainty in the extracted frequencies may be underestimated, assuming that the frequencies should be evenly spaced.
  
  For each star, we calculate the difference between the observed frequencies and the closest frequency in an evenly-spaced grid extrapolated from the fundamental frequencies assumed above, Δ f, ignoring the observed frequency if the difference between it and the closest predicted frequency is more than 0.005 d^-1. We also calculated the associated uncertainties, σ_Δ f by adding the errors of the observed frequencies and the assumed fundamental frequency in quadrature. Finally, for the entire collection of Δ f measurements for all three stars, we calculate the reduced chi-squared, 
  
      χ^2_red = 1/N∑Δ f/D σ_Δ f
  
  where D is a correction factor to account for underestimated errors, and N is the total number of frequencies in the lightcurves of all three stars. If any scatter in the values of Δ f around 0 is driven by measurement error, we can find the value of D for which χ^2_red≈ 1 — i.e., the extent to which our errors are underestimated. Doing so yields a value of D≈3, consistent with typical values of D for similar space-based photometric observations <cit.>; of course, this requires that we have correctly identified the right fundamental frequency, that the frequency spacing is indeed a regular pattern, and that these offsets are not actually due to astrophysical effects.
  
  Ultimately, until systematic correlations in data are better quantified and we are able to generate accurate asteroseismic models of FYPS, we will not be able to make a concrete determination of the cause of the offset between the observed frequencies and a regularly spaced harmonic series. The wavelet analysis used above may be able to help diagnose the behavior of the observed frequencies. Unfortunately, HD 268687 has a measured SLFV amplitude of α = 0.064±0.002, the highest of all of the discovered FYPS, and the WWZ of the lightcurve is entirely dominated by low-frequency transient features associated with this stochastic variability. The values of α are smaller for HD 269840 and HD 269902. We plot the WWZs of both stars in Figure <ref>.
  
  
  
  
  In HD 269902, most of the power in the WWZ is contained in the lower frequencies where the stochastic low frequency background dominates — unsurprisingly given that the star has the second-largest value of α out of the five FYPS (0.024±0.001). Only by altering the scaling of the image of the WWZ are we able to see the ridge associated with the lowest frequency peak, and it is only detected with low signal-to-noise. In both stars, the the lowest frequency peak in the RPS is associated with a broad band of power in the WWZ, similar to the WWZ of HD 269953, with a secondary band appearing at ∼2-3 d^-1 in HD 269840, at the approximate location of the highest peak in the RPS.
  
  
  
  § DISCUSSION
  
  
  
  
  
   §.§ Asteroseismic Modeling
  
  
  The discovery of fast pulsations in YSGs has very interesting implications for the study of the interiors of evolved massive stars. The 25 M_⊙ solar-metallicity Geneva model <cit.> has a main sequence lifetime of approximately 7 Myr. It then crosses the HR diagram in under 1 Myr, and half a Myr later, has evolved bluewards once more to become a Wolf-Rayet star <cit.>. Due to the incredibly short lifetime of YSGs on both crossings of the HR diagram, theoretical uncertainty has long stymied our understanding of massive star evolution <cit.>. Pulsation frequencies extracted from long-baseline lightcurves assembled from space-based observations are perhaps the most precise measurements we can make; typical values of f/ϵ(f) for frequencies listed in Tables <ref>-<ref> are ∼10^-4-10^-5. In better-studied stars, such precision allows for the diagnosis of incredibly complicated physics, and is truly the benchmark of testing stellar evolution theory <cit.>. No asteroseismic models for YSGs that reliably converge exist <cit.>, and so the era of precision YSG asteroseismology is not yet upon us. However, we can compare the frequencies observed in the FYPS with the characteristic Lamb and Brunt-Väisälä angular frequencies in a model YSG:
  
      S_ℓ^2    = ℓ(ℓ + 1)c_s^2/r^2
      
          N^2    = g/H_p[δ(∇_ad - ∇) + ϕ∇_μ]
  
  where ℓ is the harmonic degree, c_s is the sound speed, r is the radius of the local shell, g is the gravity, and the remaining variables are defined as in <cit.>. Modes with angular frequencies |ω| > |N| and |ω| > |S_ℓ| are mostly restored by pressure (p-modes), and modes with |ω| < |N| and |ω| < |S_ℓ| are mostly restored by buoyancy (g-modes). While rudimentary, such a comparison would illustrate the approximate regions of the stellar structure the observed pulsations will allow us to probe.
  
  
  
  To do this, we use version 12778 of Modules for Experiments in Stellar Astrophysics (MESA, ) to evolve a grid of nonrotating stellar models with initial masses between 15 and 30 M_⊙, with a spacing of 1 M_⊙, initial metallicity set to Z=0.006, and initial helium mass fraction set to Y=0.25+1.5Z=0.259. We first evolve the models through the entire pre-main sequence stage until the entirety of the star's luminosity comes from nuclear burning then introduce more elaborate physics. Convective mixing follows the mlt++ prescription from <cit.>, with overshoot following <cit.>. For mass loss, we use the “Dutch” <cit.> and “Vink” <cit.> cool/hot wind schemes respectively, and adopt an efficiency of η=0.8 for the former <cit.>. Because the luminosity in the outer layers of red supergiants can exceed the Eddington luminosity due to high opacity in these layers, we follow <cit.>, and increase the mass-loss rate in by a factor of 3 when the luminosity exceeds five times the Eddington luminosity. Finally, we use Type 2 opacities for when the star has extra C/O during and after He burning. With these controls, all but the 28 M_⊙ model successfully ran; as we focus the remainder of our analysis on a single model, we chose not to introduce additional controls for this one model. Our inlist files, including our timestep and spatial resolution controls, are available online at <https://github.com/tzdwi/TESS>.
  
  The post-ZAMS evolutionary tracks are plotted in the left panel of Figure <ref>, with the 15, 20, 25, and 30 M_⊙ tracks labeled. We note that evolutionary modeling of post main-sequence massive stars is fraught with uncertainty; different prescriptions for overshooting, or mass loss <cit.>, or different choices of input physics like binary interactions or rotation <cit.> can radically alter the evolutionary pathway a given stellar model might take. Nonetheless, we can use these evolutionary tracks to find a model that may approximate the structure and evolution of the FYPS. Interestingly, the stars with initial mass M ≥ 18 M_⊙ lose enough mass to begin turning around on the HR diagram. Models more massive than 19 M_⊙ become luminous and warm enough to encounter the yellow void, at which point the models begin to exhibit rapid changes in their luminosities and effective temperatures on ∼month to year timescales — for clarity, we do not show the post-RSG portion of the tracks after they reach an effective temperature hotter than 7000 K. The stars in our sample are plotted as grey points, and the FYPS are plotted as green stars. The dark line shows the model that passes the closest to the position of HD 269953 in the HR diagram, which has an initial mass of 23 M_⊙. We show the evolution of the interior structure of this model throughout its lifetime in the right panel of Figure <ref>. Each colored region shows the part of the star (in mass coordinates) dominated by H (purple), He (light blue), and metals (green) respectively as a function of the age of the star.
  
  The navy point and navy line in both panels of Figure <ref> indicate the pre-RSG timestep whose temperature and luminosity best match the observed values for HD 269953. At this point, the model is ∼9.0 Myr old, has a current mass of 21.6 M_⊙, and has log T_eff = 3.694, and log L/L_⊙ = 5.398. By this time, the star has begun core He fusion, has created ∼0.8 M_⊙ of C, and is losing significant mass from its envelope. Similarly, the goldenrod point and vertical line in Figure <ref> correspond to the post-RSG timestep in the same model that is closest to HD 269953. This model is 9.3 Myr old, has undergone extensive mass loss as an RSG, and has log T_eff = 3.683 and log L/L_⊙ = 5.446. It is still fusing He in its core, but has built up significant C and O mass. It's current mass is only 16.6 M_⊙
  
  
  
  We show the interior structures of both models as a function of mass coordinate in the top panels of Figure <ref>; the pre-RSG model is on the left, the post-RSG model is on the right. The density profiles, normalized by the central density ρ_c, are shown in dark yellow, while the composition is shown as the profiles of X, Y, and Z with identical colors as the right panel of Figure <ref>. We calculate the Lamb and Brunt-Väisälä frequencies in each model. The bottom panels of Figure <ref> show the logarithm of N^2 (blue) and S_1^2 (orange). The frequencies detected in the lightcurve of HD 269953 are shown as horizontal black dashed lines. Two regions exist within the pre-RSG model where g-modes are able to propagate (shaded in blue), while p-modes are able to propagate throughout the envelope of the star (shaded in orange). The inner g-mode cavity corresponds to the region outside of the core that experienced previous H burning, while the outer cavity cavity corresponds to a chemically-stratified outer envelope. In the post-RSG model, the innermost g-mode cavity has moved outward, and combined with the outer cavity. Sharp features in the stellar structure can be seen in both characteristic frequencies.[We note that the sharp steps in the composition profiles of the envelopes of both models (and the corresponding spiky behavior of N^2 in this region) are due to small discontinuities in the resolution of the MESA model, and are not real features.] 
  
  We stress that this is not an asteroseismic analysis of the star; we have made no attempt at predicting the excited frequencies in either model, and are by no means identifying the observed frequencies with p- or g-modes, let alone more complicated phenomena such as strange modes. Furthermore, the exact treatment of mixing (including semiconvection and thermohaline mixing, which we do not include in our simple models) can have an incredibly strong infuence on the pre-supernova structure of the star <cit.>. Indeed, with no reliable interior models of YSGs, we cannot identify modes in order to conduct a full asteroseismic analysis, which may allow us to constrain these physics as well as the evolutionary status of FYPS. However, this rudimentary comparison does illuminate the regions of the stellar structure that the observed pulsations might probe, as well as the drastically different interior structures and pulsational properties seen in pre- and post-RSG models that reside in quite similar regions of the HR diagram. One possible step towards mode identification is to see whether the highest amplitude pulsation frequencies scale with the observed parameters of the stars to ascertain whether the pulsations may be in the acoustic, gravity, or gravito-intertial regime. We searched for correlations between the strongest observed frequencies and log T_eff, log L/L_⊙, and log R^-2 (as a proxy for log g), but a sample of only five FYPS is insufficient to find any obvious trends. Futhermore, the available spectra are of insufficient resolution to measure v sin i, which are typically less than 10 km s ^-1 in YSGs <cit.>. With more FYPS, and higher resolution spectroscopy we may be able to conclusively determine the origin of these pulsations. Again, we emphasize that this demonstration serves only to motivate future work on these stars, and show the potential of FYPS for asteroseismology of YSGs.
  
  
  
  
  
  
   §.§ Evolutionary Status: Leftward vs. Rightward Evolution, and the Red Supergiant Problem
  
  
  Examining Figure <ref>, the FYPS appear to be descended from stars with initial masses of  M_FYPS≳20 M_⊙. As discussed above, the maximum observed luminosity of Supernovae II-P progenitors is significantly lower than the maximum luminosity of field RSGs (the red supergiant problem), implying that stars with an initial mass above M_max≈20 M_⊙ don't end their lives as red supergiants that explode <cit.>. The coincidence between M_max and M_FYPS is consistent with perhaps the most natural solution to the red supergiant problem: high mass RSGs aren't found as supernova progenitors because they do not explode as RSGs. They instead evolve bluewards on the HR diagram after losing significant mass as RSGs, attaining high enough values of L/M to excite the rich spectrum of observed pulsation modes as predicted in the post-RSG models of <cit.>; FYPS are thus post-RSG objects.[Note that the MESA models more massive than 18 M_⊙ lose enough mass as RSGs to evolve bluewards in the HR diagram. This is not seen in the low-metallicity Geneva models <cit.>. If FYPS are genuine post-RSGs, this would require an increase in the mass loss rates used in standard stellar models. RSG mass loss rates have only been directly measured in very small numbers of stars in the Galaxy <cit.>, and in even fewer LMC RSGs <cit.>. RSG mass loss is also episodic, and while constant mass loss prescriptions can reproduce the makeup of RSG populations <cit.>, they are still only an approximation to the actual mass loss histories of individual objects.]
  
  While this hypothesis offers a tantalizing solution, putting it on more solid footing would require additional evidence that FYPS are indeed in a post-RSG phase. One route is to search for evidence of strong past mass loss. The TIC contains photometry from the Wide Field Infrared Exlorer (WISE, ), in four mid-infrared bands: W1 (3.4 μm), W2 (4.6 μm), W3 (12 μm), and W4 (22 μm). For most stars, the WISE bands are in the Rayleigh-Jeans tail of their spectral energy distribution (SED), and they therefore have colors ∼0 in Vega passbands <cit.>. Deviations from 0 may be attributable to infrared molecular bands, especially in the spectra of red supergiants (Dicenzo & Levesque, in prep). However, they may also be indicative of circumstellar dust.
  
  
  We inspected the WISE colors of all stars in our sample. Only the brightest and coolest FYPS, HD 269953, has significantly larger values in each WISE color than nearby stars on the HR diagram. The remaining FYPS are mostly indistinguishable from their fellow stars, though HD 269840 also has a slightly larger W3-W4 value than other nearby stars. The three FYPS discussed in <ref> have somewhat smaller values of W3-W4, more consistent with the RSGs than the YSGs. While suggestive that the two most luminous FYPS may have some warm dust contributing to their SEDs, this is not concrete evidence that all five FYPS are in a post-RSG phase. Of course, these mid-infrared measurements are only sensitive to warm dust located relatively close to the star. Perhaps HD 269953 is a less-evolved FYPS with a warm and close CSM, while the CSM of the remaining FYPS have cooled and are undetectable. Probing cooler dust further from the stellar surface would require higher resolution imaging at longer wavelengths <cit.>. If the mechanism powering the pulsations in FYPS also operates at higher metallicity, it would be possible to find a more nearby FYPS in our Galaxy in order to perform such observations. 
  
  It is also important to note that the circumstellar dust produced by RSGs is known to be larger grain than that found in the interstellar medium, and that incorrectly accounting for “gray" extinction from this dust can lead to underestimates of L in RSGs (e.g. ). If the same is true for the circumstellar dust found around post-RSG YSGs (produced during the RSG phase), it is possible that this could lead to underestimating the L of these stars, which would in turn lead to underestimating stellar mass when comparing the stars to evolutionary tracks on the HR diagram. Quantifying this requires careful observations of these stars' circumstellar environments and dust properties.
  
  One concrete counterargument to the post-RSG hypothesis is found in HD 269110. If HD 269110 is a binary, and the pulsation amplitudes are orbitally modulated, such a binary system would have to be relatively close and eccentric, implying that the two stars have not previously exchanged significant mass. If the primary had previously experienced strong mass-loss as an RSG, the orbit would most likely have circularized instead. However, as discussed above, a companion star with a 91 day orbit would have to be ∼100 M_⊙ in order for the semimajor axis of the orbit to be larger than the stellar radius. As rotation is an equally unlikely culprit due to the low surface gravity/critical rotation speeds typical of YSGs, a more thorough characterization of HD 269110 is required.
  
  Ultimately, it is hard to draw any conclusions from such a small sample of pulsators. A total of 341 YSGs identified from <cit.>, <cit.>, <cit.>, <cit.> — all in the Milky Way — <cit.> — SMC — and <cit.> have entries in the TIC. Assuming our sample of YSGs observed by is representative of the entire population, ∼5/27=19% of YSGs are FYPS, implying that can detect ∼63 FYPS. Even if observes half as many FYPS, we would be able to better characterize the boundaries of the region in the HR diagram in which they reside, and whether the presence/behavior of FYPS is dependent on metallicity. 
  
  
  
  § CONCLUSIONS
  
  
  Our main results are summarized as follows:
  
  
    
    * We study the lightcurves of 76 cool supergiants with accurate temperatures and luminosities. For YSGs located in the Large Magellanic Cloud, these lightcurves span a time baseline of a year. We discover that low-frequency stochastic variability is ubiquitous in these stars, and rule out surface convection as the underlying cause for all but the red supergiants. This implies that this variability, also observed in main sequence O stars, is a constant feature of massive stars throughout their lifetimes.
    
    * After removing the contribution of this background variability from the periodograms of the stars in our sample, we find two regions in the HR diagram with pulsating stars. Four of these stars are candidate α Cygni variables, of which three are newly identified as such. The remaining five pulsating stars are clustered in a region of the HR diagram not previously identified as a region of instability. 
    
    * We rule out binarity, spurious signals, and chance alignment with nearby stars causing us to mistakenly find these pulsators in the same part of the HR diagram by chance, and conclude that these five stars comprise a real class, that we dub Fast Yellow Pulsating Supergiants (FYPS).
    
    * We extract pulsation frequencies from the FYPS lightcurve using a new procedure to account for the stochastic background, search for harmonics and frequency combinations from the extracted frequencies, and calculate the Weighted Wavelet Z-transform (WWZ) of the lightcurves to study the time-dependent behavior of these frequencies 
    
    * In HD 269953 and HD 269840 (and perhaps HD 269902), the lowest frequency strong periodogram peak extracted is coincident with a broad ridge of power in the WWZ that rapidly switches frequencies and amplitudes on ∼day timescales.
    
    * HD 268687, HD 269840, and HD 269902 show a broad comb of frequencies, with most extracted peaks being close to but not exactly in a harmonic chain.
    
    * One of the frequencies found in HD 269110 is split into a triplet, and the WWZ of the lightcurve at that frequency is modulated on the same timescale as the difference between peaks in the triplet. We are unable to determine whether or not this is due to rotation, binary effects, or some other cause.
    
    * We introduce the possibility that FYPS are post-RSG objects that have lost enough mass as RSGs to attain high luminosity-to-mass ratios, and excite pulsations <cit.>. This possibility is bolstered by the coincidence between the lowest estimated mass of the FYPS and the highest mass RSG progenitors. While we are unable to concretely determine the exact evolutionary status of the FYPS, future work to determine whether or not they are genuine post-RSG objects is of extreme importance, both theoretically by modeling their pulsation frequencies, and observationally by finding more FYPS. Regardless of the evolutionary status of FYPS, their pulsational properties will be of great use in the study of the interiors of evolved massive stars.
  
  
  
  The authors acknowledge that the work presented was largely conducted on the traditional land of the first people of Seattle, the Duwamish People past and present and honor with gratitude the land itself and the Duwamish Tribe. This research  was  supported  by  NSF  grant AST 1714285 and a Cottrell Scholar Award from the Research Corporation for Scientific Advancement granted to EML. This work has been carried out in the framework of the PlanetS National Centre of Competence in Research (NCCR) supported by the Swiss National Science Foundation (SNSF). 
  
  The authors thank C. Georgy for sharing the Z=0.006 Geneva evolutionary tracks.
  
  This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
  
  The data described here may be obtained from the MAST archive at
  [doi:10.17909/T9RP4V]https://dx.doi.org/10.17909/T9RP4V.
  
  This work made use of the following software and facilities:
  
  
  
  
  The Transiting Exoplanet Survey Satellite (, )
  
  Astropy v3.1.2 <cit.>, Matplotlib v3.1.1 <cit.>, MESA r12778 <cit.>, NumPy v1.17.2 <cit.>, Pandas v0.25.1 <cit.>,  Python 3.7.4, Scipy v1.3.1 <cit.>, Astroquery <cit.>
  
  
  aasjournal
  
  
  
  
  § FREQUENCIES FOUND VIA PREWHITENING
  
  
  The following tables contain the list of unique frequencies (separated by 1.5/T) found in the lightcurves of the FYPS by the prewhitening procedure described in <ref>. Frequencies that are exact harmonics of other frequencies are indicated. Combination frequencies recovered are noted in the comments of each table. However, these are only the “exact” harmonics and combination frequencies recovered to within the precision of the observed frequencies. As discussed in <ref>, frequencies with spacings that are close to but not exact harmonics and combinations are recovered in HD 268687, HD 269840, and HD 269902; in these stars, only f_4, f_4 and f_6, and f_2 respectively don't belong to these sequences of near-harmonics.
  
  lcccccccc
  
  Unique frequencies, amplitudes, phases, and formal errors for HD 269953 found via prewhitening. For each frequency, we specify the SNR as defined in text, and the height of the associated peak in the RPS at that stage of prewhitening.
  Frequency     f_j     ϵ(f_j)     A_j     ϵ(A_j)     ϕ_j      ϵ(ϕ_j)     SNR     RPS Peak Height
  
       [day^-1]     [day^-1]     [ppt]     [ppt]     [radians]     [radians]            
  f_0    1.59347960     0.00002335     0.17835463     0.00269747     2.2998     0.0151     53.8321     1187.4096 
   
  f_1^*    2.67052158     0.00009247     0.04452312     0.00266624     0.2701     0.0599     32.9642     150.6893 
   
  f_2    1.33523329     0.00006790     0.06058007     0.00266395     -0.5273     0.0440     20.1564     133.2001 
   
  f_3    1.17424693     0.00006774     0.06062141     0.00265962     0.2279     0.0439     29.5097     118.7172 
   
  *: harmonics of f_2.  
   
  
  lcccccccc
  
  Unique frequencies, amplitudes, phases, and formal errors for HD 269110 found via prewhitening. For each frequency, we specify the SNR as defined in text, and the height of the associated peak in the RPS at that stage of prewhitening.
  Frequency     f_j     ϵ(f_j)     A_j     ϵ(A_j)     ϕ_j      ϵ(ϕ_j)     SNR     RPS Peak Height
  
       [day^-1]     [day^-1]     [ppt]     [ppt]     [radians]     [radians]            
  f_0    0.55280981     0.00003493     0.16642430     0.00376526     -1.7283     0.0226     19.2635     173.6190 
   
  f_1    1.76353979     0.00013526     0.04277167     0.00374665     1.0248     0.0876     28.0123     59.3401 
   
  f_2    0.54185885     0.00007463     0.07747331     0.00374458     -1.2620     0.0483     14.7085     40.2324 
   
  f_3    0.56377424     0.00009201     0.06277635     0.00374066     -2.2209     0.0596     7.5121     29.0397 
   
   
  
  lcccccccc
  
  Unique frequencies, amplitudes, phases, and formal errors for HD 268687 found via prewhitening. For each frequency, we specify the SNR as defined in text, and the height of the associated peak in the RPS at that stage of prewhitening.
  Frequency     f_j     ϵ(f_j)     A_j     ϵ(A_j)     ϕ_j      ϵ(ϕ_j)     SNR     RPS Peak Height
  
       [day^-1]     [day^-1]     [ppt]     [ppt]     [radians]     [radians]            
  f_0    2.76530509     0.00006476     0.18020143     0.00755737     -2.9889     0.0419     63.1039     268.6123 
   
  f_1^*    3.68693149     0.00010893     0.10696509     0.00754615     2.7915     0.0705     86.6110     171.2117 
   
  f_2    1.84352713     0.00005631     0.20680066     0.00754233     -2.5874     0.0365     44.4964     148.7754 
   
  f_3    4.60860763     0.00018083     0.06428032     0.00752813     2.2392     0.1171     57.2225     88.0887 
   
  f_4    1.34234726     0.00004934     0.23553062     0.00752674     3.0867     0.0320     50.1173     87.3329 
   
  f_5    0.92187126     0.00004733     0.24496247     0.00750823     -2.2787     0.0307     8.2893     33.9434 
   
  *: harmonics of f_2. f_3 = f_1 + f_5. f_0 = f_2 + f_5.  
   
  
  lcccccccc
  
  Unique frequencies, amplitudes, phases, and formal errors for HD 269840 found via prewhitening. For each frequency, we specify the SNR as defined in text, and the height of the associated peak in the RPS at that stage of prewhitening.
  Frequency     f_j     ϵ(f_j)     A_j     ϵ(A_j)     ϕ_j      ϵ(ϕ_j)     SNR     RPS Peak Height
  
       [day^-1]     [day^-1]     [ppt]     [ppt]     [radians]     [radians]            
  f_0    2.26965893     0.00005279     0.19821747     0.00463827     -1.0031     0.0234     60.4289     490.7758 
   
  f_1    3.40445242     0.00007777     0.13350616     0.00460197     2.2648     0.0345     45.3612     460.6608 
   
  f_2    1.13468853     0.00003708     0.27915368     0.00458731     1.8206     0.0164     41.2142     332.8405 
   
  f_3    2.83703781     0.00020544     0.04966845     0.00452268     2.2047     0.0911     33.1261     61.2621 
   
  f_4    0.71676916     0.00007334     0.13905708     0.00451996     0.5576     0.0325     18.0708     57.1411 
   
  f_5    3.97206611     0.00024415     0.04160807     0.00450260     -0.9063     0.1082     32.0854     55.6234 
   
  f_6    1.38512974     0.00013130     0.07734552     0.00450113     2.4916     0.0582     24.1347     52.1294 
   
   
  
  lcccccccc
  
  Unique frequencies, amplitudes, phases, and formal errors for HD 269902 found via prewhitening. For each frequency, we specify the SNR as defined in text, and the height of the associated peak in the RPS at that stage of prewhitening.
  Frequency     f_j     ϵ(f_j)     A_j     ϵ(A_j)     ϕ_j      ϵ(ϕ_j)     SNR     RPS Peak Height
  
       [day^-1]     [day^-1]     [ppt]     [ppt]     [radians]     [radians]            
  f_0    2.90337860     0.00006011     0.12717671     0.00455148     -1.2740     0.0358     62.4969     333.8202 
   
  f_1    1.45175166     0.00003744     0.20368400     0.00454088     -1.4444     0.0223     53.9124     307.3930 
   
  f_2    0.08163766     0.00002411     0.31421154     0.00451126     -2.6455     0.0144     5.5540     21.4397 
   
   
  
  
  
  
  
  