Recovery of TESS Stellar Rotation Periods Using Deep Learning

We used a convolutional neural network to infer stellar rotation periods from a set of synthetic light curves simulated with realistic spot evolution patterns. We convolved these simulated light curves with real TESS light curves containing minimal intrinsic astrophysical variability to allow the network to learn TESS systematics and estimate rotation periods despite them. In addition to periods, we predict uncertainties via heteroskedastic regression to estimate the credibility of the period predictions. In the most credible half of the test data, we recover 10%-accurate periods for 46% of the targets, and 20%-accurate periods for 69% of the targets. Using our trained network, we successfully recover periods of real stars with literature rotation measurements, even past the 13.7-day limit generally encountered by TESS rotation searches using conventional period-finding techniques. Our method also demonstrates resistance to half-period aliases. We present the neural network and simulated training data, and introduce the software butterpy used to synthesize the light curves using realistic star spot evolution.


INTRODUCTION
Stellar rotation is fundamentally linked to the structure and evolution of stars.In the decade since Kepler, much has been learned about rotation, feeding into asteroseismology, empowering gyrochronology, and changing the way we think about stellar evolution codes.Rotation period estimates are made possible through a variety of methods.Historically, spectroscopy enabled estimates of rotation velocity due to Doppler red/blue shift from the receding/approaching halves of the stellar disk.The projected rotation velocity could then be used to compute an upper limit on the period if the stellar radius was known.Missions like CoRoT (Baglin et al. 2006) and Kepler (Borucki et al. 2010) have shifted the paradigm: the majority of period estimates now employ photometry instead of spectroscopy.This works particularly for stars which, like the Sun, exhibit magnetic dark and bright spots that induce periodic variations to Corresponding author: Zachary R. Claytor zclaytor@hawaii.edu the light curves as the stars rotate.Several techniques have been developed in recent years to extract rotation information from spot-modulated stellar light curves.Namely, Lomb-Scargle periodograms (Marilli et al. 2007;Feiden et al. 2011), autocorrelation analysis (McQuillan et al. 2013(McQuillan et al. , 2014)), wavelet transforms (Mathur et al. 2010;García et al. 2014), Gaussian processes (Angus et al. 2018), and combinations of these (Ceillier et al. 2017;Santos et al. 2019;Reinhold & Hekker 2020) have all been used to infer rotation periods from light curves.
Period-finding methods have paved the way for large studies of stellar rotation.Applied to CoRoT and Kepler, these techniques have delivered tens of thousands of rotation period estimates, which in turn have been used to advance our understanding of stellar and Galactic evolution (e.g., McQuillan et al. 2014;van Saders et al. 2016van Saders et al. , 2019;;Davenport 2017;Claytor et al. 2020;Amard et al. 2020).The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) stands to increase the number of inferred periods by an order of magnitude in its ongoing all-sky survey.
Rotation studies have also brought to light the limitations of period detection methods.For example, con-ventional methods are subject to aliases, and they still struggle to detect rotation in quiet, Sun-like stars (Mc-Quillan et al. 2014;van Saders et al. 2019;Reinhold & Hekker 2020).Furthermore, the traditional methods do not necessarily reveal a star's true period.Rather, they reveal the period(s) of latitudes at which star spots form, which may rotate faster or slower than the star's equator due to surface differential rotation.
Finally, the systematics of TESS have made traditional period searches difficult (Oelkers et al. 2018;Canto Martins et al. 2020;Holcomb 2020;Avallone et al. 2021, in prep.).The lunar-synchronous orbit of TESS has a 13.7-day period, and the telescope is subject to background variations from reflected sunlight causing periodic contamination that is difficult to remove.As a result, dedicated rotation studies struggle to obtain reliable periods longer than about 13 days (e.g., Canto Martins et al. 2020;Holcomb 2020;Avallone et al. 2021, in prep.).New, data-driven methods are needed to overcome these systematics and recover periods.
Deep Learning is relatively new to astronomy, but in a short time deep learning methods have proven to be valuable at mining information from large data sets.Neural networks, and in particular Convolutional Neural Networks (CNNs), efficiently extract information from time series, spectra, and image data.The strength of CNNs comes from their local connectivity, which incorporates the knowledge that neighboring input points are highly correlated.Examples of success using CNNs in astronomy include Hezaveh et al. (2017), who used CNNs to characterize gravitational lenses from image data.Within the realm of stellar astrophysics, Guiglion et al. (2020) used the same techniques to obtain stellar parameters from spectra.Moreover, Feinstein et al. (2020) and Blancato et al. (2020) used CNNs to infer stellar parameters and flare statistics from light curves.
Using convolutional neural networks, we predict stellar rotation periods from wavelet transforms of light curves.The use of supervised machine learning requires the existence of a training data set for which the target, in this case the rotation period, is known.This is not yet possible with TESS due to the difficulties of obtaining reliable periods using traditional techniques.Furthermore, while there is some overlap between TESS and Kepler, the TESS observations of the Kepler field are short, spanning only 27 days.This is enough to recover and validate short rotation periods (Blancato et al. 2020), and possibly some subset of longer periods (Lu et al. 2020), but not enough to be useful for the broader population of stars in our Galaxy.Moreover, even large rotation samples from Kepler are likely contaminated with mismeasured periods (Aigrain et al. 2015).Us-ing a training set of periods obtained with conventional techniques risks imprinting this contamination onto the neural network.To avoid this, we followed the approach of Aigrain et al. (2015) and used a set of synthetic light curves generated from physically motivated star spot emergence models.This is an example of simulationbased inference (e.g., Cranmer et al. 2020), wherein we simulate a physical process and use machine learning to address the inverse problem of inferring the rotation.
We introduce butterpy1 , an open-source Python package designed to simulate realistic star spot emergence and synthesize light curves, followed by a description of the input physics of the simulations.After describing our training set, we outline our convolutional neural network and the methods we use to train, validate, and test the network.We evaluate our trained neural network on synthesized data sets spanning different period ranges to identify for what periods the network is most predictive.Next, we discuss the network's performance on a small set of real light curves for which rotation periods are known.We also compare our network predictions to periods recovered using conventional methods before finally concluding with thoughts on the feasibility of our methods to recover stellar rotation periods from real TESS light curves.

SYNTHETIC LIGHT CURVES: BUTTERPY
Synthesized light curves have several advantages over observed light curves: (1) the true, equatorial period of the simulated star is known, rather than an estimate of the period (which may be wrong), (2) data of any length and cadence can be synthesized, and (3) other physical properties like spot characteristics, differential rotation, and surface activity are known and can be independently probed.
To simulate light curves, we developed butterpy, a Python package designed to generate realistic, physically motivated spot emergence patterns faster than conventional surface flux transport codes.The name butterpy comes from the butterfly-shaped pattern of spot emergence with time exhibited by the Sun (e.g., Hathaway 2015).We built upon the software developed by Aigrain et al. (2015), which relied on the flux transport models of Mackay et al. (2004) and Llama et al. (2012) to generate spot emergence distributions which were in turn used to compute light curves.The original model of Mackay et al. (2004) was designed to reproduce spot emergence patterns of the Sun as well as Zeeman Doppler images of the pre-main-sequence star AB Do-radus.Later, Llama et al. (2012) used this model in tandem with exoplanet transit observations to trace the migration of active latitude bands across the surfaces of stars.Aigrain et al. (2015) used the light curves generated from these spot distributions to test the recovery rates of various period detection techniques, which we seek to emulate.We discuss the method and assumptions here for clarity.

Spot Emergence and Light Curve Computation
There are several variables to consider regarding the emergence of star spots and their effect on the star's light curve, including the latitudes and rates of emergence, the spot lifetimes, and the rotation speed at the latitude of emergence if the star rotates differentially.While observations of these on stars other than the Sun are limited, they are very well characterized for the Sun (Hathaway 2015, and references therein).For our model, we therefore start with the characteristics that are known for the Sun and allow them to vary.

Location and rate of spot emergence
The latitudes of spot emergence on the Sun vary with the Sun's 11-year activity cycle.At the beginning of a cycle, spots emerge within active regions at high latitudes (λ ≈ ±30 • , Hathaway 2015), and the latitude of emergence migrates toward the equator throughout the rest of the cycle.Before the cycle ends, new spot groups begin forming again at high latitudes, indicating some amount of overlap between consecutive cycles.The repeating decay in spot latitude with time gives rise to a butterfly-like pattern known as a "butterfly diagram".Butterfly patterns have been observed in other stars as well (Bazot et al. 2018;Nielsen et al. 2019;Netto & Valio 2020), but some stars show a random distribution of spot emergence latitudes with time (e.g., Mackay et al. 2004).The width of active latitude bands has also been shown to differ even for Sun-like stars (Thomas et al. 2019).In our model, we allow for either a random or butterflylike spot emergence between a minimum and maximum latitude that are unique for each star.For the butterfly pattern, spots begin the cycle emerging at a latitude λ max , decaying exponentially with time to latitude λ min at the end of the cycle (Hathaway 2011).
As for longitude, spots on the Sun tend to emerge in groups, either next to existing spots, or in some cases antipolar to existing spots (Hathaway 2015).Less often, spots will emerge at random longitudes, not necessarily associated with any existing spot groups.We respectively refer to these two cases as correlated and uncorrelated active regions.In our simulations, we follow the approach of Aigrain et al. (2015), dividing the stellar surface into 16 latitude and 36 longitude bins; the probability of spot emergence is distributed across these bins.To account for the relative likelihood of correlated and uncorrelated emergence, bins already containing active regions are assigned a higher probability of emergence.
The rates of sunspot emergence change with spot area and with time throughout the activity cycle.Schrijver & Harvey (1994) expressed the number of spots emerging in area interval (a, a + da) and time interval (t, t + dt) as r(t)a −2 da dt, where r(t) represents the time-varying emergence rate amplitude, the active region area a is in square degrees, and t is the time elapsed in the activity cycle, ranging from zero to one.For the time dependence, spots emerge very slowly at the beginning, more rapidly in the middle, and slowly again at the end (Hathaway et al. 1994).Mackay et al. (2004) modeled this using a squared sine function: r(t) = A sin 2 (πt).The activity level A is an adjustable scale factor controlling both the average rate of spot emergence and the amplitude of light curve modulation of a single spot.It is defined such that A = 1 for the Sun.

Latitudinal Differential Rotation
We define heliographic longitude such that φ = 0 always faces the Earth.As a consequence, spots move in longitude as the star rotates.The Sun rotates more rapidly near the equator than at the poles, a phenomenon known as "latitudinal differential rotation" (henceforth just "differential rotation").While differential rotation is more difficult to observe on other stars, some stars may exhibit "anti-solar" differential rotation, wherein the equator rotates more slowly than the poles.This has been observed particularly in slowly-rotating stars (e.g., Rüdiger et al. 2019).In our model, we allow for solar-like, anti-solar, and solid-body rotation.Following Aigrain et al. (2015), we model the differential rotation profile as (1) Here, λ k and φ k denote the heliographic latitude and longitude of spot k. t max,k is the time at which the spot achieves maximum flux, Ω is the angular velocity at latitude λ k , and α is the differential rotation shear parameter.To include anti-solar, solid-body, and solarlike profiles, we allow α to range from -1 to 1.

Spot-Induced Flux Modulation
Once spot emergence is determined, we simulate spot evolution and flux modulation based on the simplified model of Aigrain et al. (2012Aigrain et al. ( , 2015)).They take the photometric signature δF k (t) of a single spot k to be where β k (t) is angle between the spot normal and the line of sight, accounting for projection on the stellar surface.The inclination i is the angle between the rotation axis and the line of sight, and λ k and φ k (t) are again the latitude and longitude of the spot.The factor f k is the amount of luminous flux removed if spot k is observed at the center of the stellar disk.Aigrain et al. (2015) used an exponential rise and decay to model the rapid rise and slow decay of single spots, but we employ a twosided Gaussian to exhibit smoothness while preserving the same emergence and decay behavior: where f max k is the flux removed by spot k at the time of maximum emergence t max,k , τ is the relevant emergence or decay timescale, and τ spot is a dimensionless parameter used to relate the emergence and decay timescales to the equatorial rotation period P eq .Like Aigrain et al. (2015), we parametrized the emergence and decay timescales as multiples of the equatorial rotation period.The form of the emergence timescale was chosen so that, in general, emergence is five times faster than decay, with a minimum possible emergence timescale of two days.In the simple model of Aigrain et al. (2012), f max k takes into account the spot area and contrast, but the model of Aigrain et al. (2015) relates this factor to the strength of the magnetic field: where A is the activity level, and the constant is chosen such that A = 1 reproduces approximately Sun-like behavior.
With this expression, the single-spot luminous flux modulation is proportional to the strength of the radial magnetic field at that spot.The magnetic field strength or magnetic flux is proportional to the area of the active region, which van Ballegooijen et al. (1998) derive using the angular width of magnetic bipoles emerging from the active region: where B ± r is the radial component of the magnetic field near either the positive or negative bipole, B max is the initial peak magnetic field strength in the active region, β init is the angular width of a single bipole, β 0 is the angular width (in degrees) of the bipole at the time the active region is inserted into the model, accounting for diffusion, and β ± is the heliocentric angle between a field point and one of the bipoles.
van Ballegooijen et al. (1998) assumed that the bipole width β init is proportional to the angular separation between the positive and negative poles, which they call ∆β, with a proportionality factor of 0.4.Assuming spots form within ten degrees of the active region bipoles, the value of the exponential factor differs from unity by less than one percent.For this reason, we approximate the exponential factor as unity.Thus, at the location of a star spot, Equation ( 5) simplifies to Combining this with Equation ( 4), van Ballegooijen et al. ( 1998) and Mackay et al. ( 2004) consider a range of bipole widths from about 3.5 • to 10 • .The distribution of ∆β is the same for every star, so ∆β k k is effectively constant.Putting it all together, we have a final system of equations to describe the change in luminous flux from a single spot: (8) This depends on ∆β k , λ k , φ k (t), and t max,k , which are unique to each spot; and A, τ spot , P eq , and i, which are unique to each star.

Training Set
Using the model in butterpy, we generated one million light curves at thirty-minute cadence and one-year Note-We adopted the distributions used by Aigrain et al. (2015) with minor modifications: (1) we sampled a broader range of periods and activity levels, (2) we used a uniform distribution of periods so as not to impart unwanted bias on the neural network prediction, (3) we include anti-solar differential rotation by allowing the shear parameter to be negative.duration to match the TESS Full-Frame Images (FFIs) in the continuous viewing zones (CVZs).The simulation input parameters are listed in Table 1.We sampled periods uniformly from the range [0.1, 180] days.The period range was chosen to be as wide as possible to simulate the fastest-rotating stars (P eq ≈ 0.1 day) while also capturing anything that would go through at least two rotations under observation in the TESS CVZs (the total baseline is ∼360 days, so an object with P eq = 180 will go through exactly two rotations in that time).We chose the remaining distributions and ranges to reflect those of Aigrain et al. (2015), with minor adjustments in the ranges of activity level and differential rotation shear to search a broader parameter space.Our input distributions assumed no relation between rotation period and activity level.Figure 1 illustrates an example simulation, showing the distribution of spots on the surface as well as their impact on the observed light curve.

TESS Noise Model
To ensure the training light curves properly emulate real TESS light curves, the training set must exhibit TESS -like noise.Aigrain et al. (2015) used light curves from quiescent Kepler stars to achieve this.In their study, a sample of stars from McQuillan et al. (2014) with no significant period detection served as the quiescent data set.Because there are no existing bulk period measurements for stars in the CVZs, we must find another means of simulating TESS noise.
While TESS is a planet-finding mission, the southern CVZ contains thousands of galaxies which should have roughly constant brightness with time.Any changes in the light curves of these galaxies would be due solely to TESS instrument systematics.Thus, the galaxy light curves should reasonably resemble light curves of quiescent stars in TESS.
We selected roughly 2,000 galaxies in the southern CVZ with T mag ≤ 15 as our quiescent sample, removing a handful of galaxies known to be active and in the Half-Million Quasars catalog (Flesch 2015).We queried FFI cutouts from the Mikulski Archive for Space Telescopes (MAST) using Lightkurve and TESScut (Lightkurve Collaboration et al. 2018;Brasseur et al. 2019).Then, we performed background subtraction and aperture photometry on each source using Lightkurve regression correctors, following Lightkurve Collaboration (2020).To summarize, aperture masks were chosen using the create threshold mask function in Lightkurve.This method selects pixels with fluxes brighter than a specified threshold number of standard deviations above the image's median flux value.We specified thresholds based on the target's brightness to exclude background pixels from the aperture.Once the raw light curve was computed, the regression correctors fit principle components of the time-series images and subtracted the strongest components from the raw light curve.All sector light curves for a source were then median-normalized and stitched together to form the final "pure noise" light curve.
The galaxy light curves were linearly interpolated to each TESS cadence to fill gaps, whether for missing observations or entire missing sectors.Cadences missing at the beginning or end of the light curve were filled with the light curve's mean flux value.Finally, a galaxy light curve was chosen at random to be convolved with each of the synthetic light curves, yielding our final set of simulated TESS -like light curves.We divided the set of 2,000 galaxies into two sets of 1,000: one set to be convolved with light curves from the training partition, and one for the validation and test partitions (see Section 4 for more about data partitioning).

DATA PROCESSING/WAVELET TRANSFORM
There are several options for input to a neural network to predict rotation periods.One could use the light curve directly; Blancato et al. (2020) suggest this as the best way to obtain periods using neural networks without loss of information.However, using the light curve as input means that the information conveying periodic- Note-We use three 2D convolution layers, each with ReLU activation and maxpooling.Our implementation uses 2D max-pooling with a 1-dimensional kernel to achieve pooling in the time dimension but not the frequency dimension.This choice preserves frequency resolution but achieves a small amount of translational invariance in the time dimension.The output of the convolution block is flattened to a 1-dimensional array and passed through three fully connected (dense) layers, with ReLU and finally softplus output, to yield two numbers: the rotation period and its uncertainty.
ity is temporally spread out.While neural networks can certainly learn to predict periods this way, a frequency representation concentrates the period information to one location in input space.Lomb-Scargle periodograms (Lomb 1976;Scargle 1982;Feiden et al. 2011) and autocorrelation functions (McQuillan et al. 2013(McQuillan et al. , 2014) ) are two tried-and-true methods of period estimation that have some promise as input to neural networks.While these methods are effective at concentrating periodicity information to one location, real stars' observed periodicity can change with time due to differential rotation.Lomb-Scargle and autocorrelation methods average over these changes, potentially blurring out interesting evolution.The continuous wavelet transform (Torrence & Compo 1998) has also been used to identify rotation periods from stellar light curves (Mathur et al. 2010;García et al. 2014;Santos et al. 2019) and it has the bonus of elucidating changes in periodicity with time.
We chose this method to localize periodic information while allowing the tracing of spot evolution.We used the continuous wavelet transform implemented in SciPy (Virtanen et al. 2020) with the power spectral density correction of Liu et al. (2007).Using a Morlet wavelet, we computed wavelet power spectra for both the noiseless and noise-injected light curves in our training set.Examples of both noiseless and noisy power spectra are shown for the same simulated star in Figure 2. We then rebinned the power spectra to 64×64 pixels and saved them as arrays for fast access.Larger binned sizes were tested (e.g., 128 × 128) and showed no significant change in performance.
We ran several tests with the period axis of the wavelet power spectra, trying maximum periods of 128, 150, and 180 days, before settling on 180 days (i.e., half the observing window) for the final data products.In the tests with 128 and 150 days, the neural network had some success predicting periods longer than the maximum visible period in the periodogram, even for the noisy data.This suggests that neural networks can predict periods even when the period at maximum power is beyond the range of the plot, consistent with the results of Lu et al. (2020).This is encouraging for period predictions for stars outside the TESS continuous viewing zones, where observations are substantially less than a year in duration.In the end, we chose 180 days as the maximum value on the period axis to preserve the strongest rotation signals in as many objects as possible.
In addition to butterpy, our final data products will be made publicly available and include (1) the noiseless, synthesized light curves, (2) the normalized TESS galaxy light curves, and (3) the binned wavelet transform arrays for both the noiseless and noisy light curves.

CONVOLUTIONAL NEURAL NETWORK
We used a convolutional neural network to predict rotation periods from wavelet transforms.Table 2 outlines the CNN architecture.We used a sequence of 2D convolution layers with rectified linear activation ("rectifier" or "ReLU") followed by 1D max-pooling in the time dimension.The ReLU activation function has the form f (x) = max(0, x).Its nonlinearity allows the model to represent complex functions, and ReLU learns faster than other nonlinear activation functions.Max-pooling is used to down-sample input and impart a small amount of translational invariance.The shapes of the convolution and pooling kernels were chosen to impart equivariance in the frequency dimension (no pooling, since frequency is what we want to estimate) and translational invariance in the time dimension.This means that the CNN will treat periodic signals the same regardless of when they occur in the wavelet power spectrum (see Ch. 9 of Goodfellow et al. 2016).
The output of the convolution layers is then flattened to one dimension and fed into a series of three fully connected layers, also with ReLU activation.The final layer uses softplus activation, which has the form f (x) = ln(1 + e x ).A smooth approximation to the rectifier, softplus activation ensures positive output while preserving differentiability.The final layer outputs two numbers, which represent the rotation period and the period uncertainty.
The 2D wavelet power spectra were used as input to our neural network, while the corresponding model rotation periods served as the target output.Target periods were min-max scaled over the entire data set to the range [0, 1].Each power spectrum array was min-max scaled to the range [0, 1] separately-using the min and max over the entire data set suppressed lower-amplitude signals and substantially impaired performance.
Our full data set of one million examples was partitioned into three sets for model training, validation, and testing.The training set consisted of 80% and was used to fit the model weights.The validation set (10%) was used for early stopping (we stop training when the validation loss does not improve over a window of 10 training epochs to avoid overfitting) and for choosing the optimal hyperparameters.The test set (10%) was used for final model evaluation.
We used the Adam optimizer (Kingma & Ba 2014), which allows for a variable learning rate, with negative log-Laplacian likelihood as the loss function.This loss function allows us to predict both the rotation period and its error (a process known as heteroskedastic regression), indicating which period predictions are more reliable than others.It has the form where b is taken to represent the predicted uncertainty.
Maximizing the log-likelihood of the Laplace distribution is equivalent to minimizing the mean absolute error instead of the mean-squared error, or predicting the median period instead of the mean under uncertainty.This also means that in cases where the neural network cannot predict with high confidence, predictions will be biased toward the median of the period range.Formally, the Laplace distribution has variance 2b 2 and standard deviation √ 2b, but we use b to represent the uncertainty for simplicity since we use it only to determine the relative credence of predictions.Thus, our predicted uncertainties should not be considered statistically formal.
With 800,000 input-output pairs in the training set, our model takes roughly 3 hours until fully trained on a single NVIDIA RTX2080.Once trained, evaluation on the test input of 100,000 wavelet power spectrum plots takes less than a minute.

RESULTS
We trained and evaluated the neural network on yearlong simulations of both noiseless and noise-injected wavelet transform images.We additionally used Lomb-Scargle periodograms, autocorrelation functions, and wavelet transforms to obtain independent period estimates from the noisy data.Aigrain et al. (2015) performed blind injectionrecovery exercises on synthesized Kepler -like light curves to assess the reliability of conventional perioddetection methods.On average, the teams recovered periods with 10% accuracy in ∼70% of cases in which periods were obtained.We adopt this 10% accuracy threshold as our success metric, which we designate "acc10."In addition to acc10, we also quantify results with "acc20," mean absolute percentage error (MAPE), and median absolute percentage error (MedAPE), defined as follows.
If we define the absolute percentage error of example i to be i = |P pred,i − P true,i |/P true,i , then our recovery metrics are where H(x − i ) is the Heaviside or unit step function.Before commenting on our period recovery, it is important to note the differences in our light curve sample from that of Aigrain et al. (2015).The most important differences are in the range of activity level and the light curve pre-processing.Our sample spans a wider range of activity levels, ranging from 0.1 to 10 times Solar, as opposed to 0.3 to 3 times Solar in Aigrain et al. (2015).The logarithmic scale of the distribution ensures that the increase in range evenly adds examples to the high-and low-activity ends.Thus, despite having higher-amplitude examples in our sample, there should be enough lower-amplitude examples to compensate, preserving the comparability of our summary statistics to those of Aigrain et al. (2015).Light curve pre-processing differs because of the differences in the Kepler pipeline and our custom TESS FFI pipeline.In principle, the Kepler pipeline more aggressively removes systematics, so the Aigrain et al. ( 2015) simulated light curves are cleaner than ours.

Period recovery using conventional methods
We recovered periods from our sample of noiseinjected light curves using  3.In each panel, objects falling within 10% of the line y = x are successfully recovered according to our metric.All three methods struggle to recover periods longer than about 50 days.Longer than this, LSP and GWPS mistakenly recover signals approaching 30 days in period, which we suspect represents the TESS sector length of 27 days.LSP and GWPS are also susceptible to half-period aliases, which fall along the line y = 1 2 x.ACF is less susceptible to half-period aliases and is the most successful method overall.However, the ACF and GWPS often misidentify signals at 5, 13, and 27 days (all well-known frequencies associated with TESS telescope systematics) as the rotation period.Interestingly, the ACF has small pockets of higher success at integer multiples of 27 days, beginning at 54 days.These occur when a star rotates an integer number of times in an integer number of sectors.The sector-tosector stitching affects subsequent revolutions the same way, so the signal is preserved enough for the autocorrelation function to detect.
In general, recovery was better for targets with higher light curve amplitude for all three methods, as one would ."acc10" represents fraction of periods recovered to within 10% accuracy, while "acc20" is the recovery to within 20%.ACF has the highest overall success, but the recovery worsens significantly at periods longer than about 30 days.
expect.The recovery rates also improve when limited to shorter periods.We have assumed no rotation-activity relation, so the improved recovery at shorter periods occurs when more rotations are observed in the given baseline, resulting in higher power in the periodograms.Moreover, at shorter periods rotation signals are less easily lost in the telescope systematics.If we limit to periods between 0 and 50 days, acc10 and acc20 improve to 43% and 59% for LSP, 43% and 56% for GWPS, and 36% and 47% for ACF.Thus, for short periods, LSP achieved the highest rate of success.

CNN performance on noiseless data
Our neural network's predictions on noiseless test data are shown in the left panel of Figure 4.The predicted periods have a mean absolute percentage error of 14% and a median absolute percentage error of 7%.61% of periods were successfully recovered to within 10%, setting the bar for comparison to results for the noiseinjected data.

CNN performance on noisy data
We present the neural network predictions on the noise-injected test data in the right panel of Figure 4.The predicted periods have a mean absolute percentage error of 246% and a median absolute percentage error of 24%.Only 28% of periods are successfully recovered to within 10%.The horizontal band at predicted period of 90 days represents simulated stars for which the network could not predict the period at all, instead assigning it the median of the period range.
The addition of TESS -like noise severely inhibits the performance of the neural network.Like the conventional methods, the neural network predictions are more accurate at shorter periods.When limiting to periods of 50 days or less, the median absolute percentage error is 12%, and 44% of targets are recovered to within 10%.The introduction of noise to the light curves also affects the amplitudes at which the network is most reliably predictive.The left panel of Figure 5 shows network recovery rate as a function of amplitude R per (as defined by Basri et al. 2011) and equatorial rotation period.Here, "recovered" means the prediction is within 10% of the true period.As expected, the network performs better with higher-amplitude modulations, where the stellar signals are more easily picked out of the noise.
In addition to the rotation period, our choice of loss function allows us to predict the period uncertainty.This value is a metric for how well the network is predicting the period.The right panel of Figure 5 shows the predicted uncertainty versus period and amplitude.The predicted uncertainty, like the recovery rate, is better at higher amplitudes.Since the predicted uncertainty correlates with the recovery rate, the predicted uncertainty is a reliable metric for successful period recovery without already knowing the period.We can then use the predicted uncertainty to select a part of the sample recovered to a desired accuracy.
For our analysis, we selected the half of the test set with the lowest predicted fractional uncertainty.The median predicted uncertainty for the sample before the cut was σ pred /P pred = 0.35.The period recovery for the best-predicted half of the sample is shown in Figure 6.The cut removed the horizontal band at predicted period of 90 days, and all summary statistics were improved.46% of periods were correctly predicted to within 10%, and 69% were accurate to within 20%.The predicted periods had a mean absolute percentage error of 57% and a median absolute percentage error of 11%.A few targets with incorrectly predicted period between 100 and 150 days remained after the cut.These had low predicted fractional uncertainty due to their large predicted period compared, so they made the cut despite being poorly predicted.They accounted for about 4% of the sample after the cut.
As with all other methods, the CNN performed better on the noisy data when limited to targets with periods less than 50 days.For this subset of the sample, the median predicted fractional uncertainty was σ pred /P pred = 0.2.Making the same cut as before (using the median fractional uncertainty of 0.2), the recovery of short-period stars improved to acc10 of 58% and median absolute percentage error of 8%.Table 3 shows the complete summary of our recovery results.

DISCUSSION
We have demonstrated that convolutional neural networks are capable of extracting period information from noisy light curves or, more precisely, transformations of noisy light curves.Our model also predicts the uncertainty in the period estimate, enabling us to see where the network is most successful and determine which period predictions are most reliable.Here we discuss the strengths and weaknesses of our approach and compare them to those of conventional period detection methods.We then comment on the prospects of estimating rotation periods from real TESS light curves.

Strengths and weaknesses of deep learning approach
Our CNN outperformed conventional techniques in the recovery of rotation periods for the same underlying sample of synthetic light curves.Whereas the conventional methods failed to recover periods longer than ∼2 TESS sectors, our method successfully recovered simulated star periods across the full simulation range of  Predicted periods have mean absolute percentage error of 14%, median absolute percentage error of 7%, acc10 of 61%, and acc20 of 81%.Right: Period predictions from noise-injected data, where recovery is significantly worse.Predicted periods have mean absolute percentage error of 246%, median absolute percentage error of 24%, acc10 of 28%, and acc20 of 45%.The horizontal band at 90 days represents targets where the model struggled to predict the period.In these cases, the prediction was biased toward the distribution median, or 90 days.2014), to gauge where stars from Kepler would fall.Left: Period recovery rate as a function of period and amplitude for the noise-injected data.The neural network performs better at higher amplitudes, where the rotation signal overpowers instrumental noise.Right: The same data, now colored by the neural-network-predicted fractional uncertainty in rotation period.The prediction is more certain for higher amplitudes.Furthermore, the prediction is most certain in the region with the highest recovery rate, indicating the predicted uncertainty is a reliable metric for period recovery without already knowing the true period.Note-Recovery metrics for both the full 0.1-180-day period set, and for the subset with Prot ≤ 50 days.All methods perform better on shorter-period stars, but our neural network consistently outperforms the conventional techniques on simulated light curves with real TESS systematics. .Period recovery for the half of the test set with the lowest predicted fractional uncertainty.Predicted periods have mean absolute percentage error of 57%, median absolute percentage error of 11%, acc10 of 46%, and acc20 of 69%.The predicted error cut removes the cluster of predicted periods at 90 days, giving credence to our cut to remove spurious period predictions.The cloud of objects with short true periods and predicted periods between 100 and 150 days have low fractional error because their predicted periods are large, but they account for only 4% of the objects remaining after the cut.
0.1-180 days.Simulated stars in the range of periods yet unprobed by TESS -13.7 days up to 90 days and beyond-were recovered with the highest success rate.The recovery rate trails off at each end of the range (P rot < 10 days and P rot > 170 days) because of the choice of loss function: predicting the median under uncertainty biases predictions toward the median of the ensemble distribution and away from the ends of the range.
The challenge for classic period-recovery methods in TESS light curves is mostly due to sector-to-sector stitching and the presence of scattered moonlight (repeating every 27 and 13.7 days, respectively).Other effects such as temperature changes and momentum dumps appear at periods of 1.5, 2, 2.5, 3, 5, and 13.7 days (Vanderspek et al. 2018).All these effects combine to leave periodic imprints in the data that dominate stellar rotation signals and are difficult to remove.All three of our conventional method tests significantly misidentified 27-days as the rotation period.Different methods latch onto different signals as well.For example, ACF has significant misidentifications at 2.5 and 13.7 days and a weak twice-period alias, while WPS mistakes 1and 5-day signals as the rotation period.LSP mistakes these high-frequency signals less often, but often falls prey to half-period aliases, as does WPS.
It is noteworthy that, unlike with LSP and WPS, our neural network has no significant misidentification of half-period aliases or the high-frequency systematic aliases.This is especially encouraging since we use WPS as the basis for our training data.The fact that these aliases certainly exist in the training set but are not chosen as the period supports our claim that neural networks can learn and bypass systematic and false-period signals.At the very least, if the rotation period is ambiguous, the network will predict a large uncertainty, allowing us to throw away the prediction.Our results suggest that convolutional neural networks can learn systematic effects and regress rotation periods despite them, a significant step towards enabling large stellar rotation studies with TESS.

Comparisons to other period recovery attempts
Our results suggest that this method would be on par with or better than other recent attempts to estimate > 13-day rotation periods from TESS light curves.Canto Martins et al. (2020) used a combination of Fast Fourier Transform, Lomb-Scargle, and wavelet techniques to estimate periods for 1,000 TESS objects of interest.They obtained unambiguous rotation periods for 131 stars, but all were shorter than the 13.7-days TESS orbital period.Lu et al. (2020) trained a random forest (RF) regressor to predict rotation periods from 27-day sections of Kepler light curves coupled with Gaia stellar parameters.They then evaluated the trained model on single sectors of TESS data for the same stars.Despite the stark differences in light curve systematics, they were able to recover rotation periods up to ∼50 days with 55% accuracy, which is on par with the 57% mean uncertainty achieved by our model.There are caveats to this comparison, however.First, the RF regressor relied primarily on effective temperature and secondarily on the light curve variability amplitude; light curve periodicity was not used for the period regression.Second, Lu et al. (2020) used two-minute cadence, Pre-search Data Conditioned Simple Aperture Photometry (PDCSAP) TESS light curves, while our light curves were thirtyminute cadence, and our processing pipeline was more similar to simple aperture photometry (SAP).PDCSAP light curves are subjected to much heavier detrending than those produced by SAP methods.Finally, Lu et al. (2020) used real TESS data, while we used simulated light curves.Each set comprises different distributions of rotation period, amplitude, and other important parameters.Because of these caveats, we encourage the reader to take caution when comparing the results of these studies.

Prospects for measuring periodicity in TESS
We have so far demonstrated the ability to recover photometric rotation periods from simulated TESS -like stellar light curves using deep learning.But the biggest question remains: can we reliably measure long periods from real TESS data?This is a difficult question to answer definitively for several reasons.First, validation of any method requires a set of real stars for which rotation periods are already known.The ideal data set for comparison is Kepler, where tens of thousands of periods have been recorded (McQuillan et al. 2014;Santos et al. 2019).Unfortunately, the overlap between TESS and Kepler is small: most Kepler stars were observed for only a single sector at a time in TESS.With only a 27-day baseline, it is impossible to validate a method of obtaining long periods.Stars in the TESS CVZs were monitored continuously for almost a year, but only a handful of these stars have previously known rotation periods.
Despite the limitations, we attempted to recover rotation periods for a handful of stars observed by Kepler, the Kilodegree Extremely Little Telescope survey (KELT, Pepper et al. 2007), the MEarth Project (Berta et al. 2012), and the All-Sky Automated Survey for Supernovae (ASAS- SN Shappee et al. 2014;Kochanek et al. 2017).

Kepler
We targeted the few stars in the Kepler field that had two consecutive sectors in TESS, offering a baseline of roughly 50 days.We simulated an entirely new training set with periods spanning 0.1 to 50 days, using a sample of galaxies in the Kepler field as the noise model.With a 50-day baseline, only periods of up to 25 days might be recoverable, as timescales longer than this may be dominated by edge effects in the wavelet transform.Even so, our network was unable to recover Kepler periods reliably.

KELT
We similarly targeted 106 stars in the TESS SCVZ also observed by the KELT survey.Oelkers & Stassun (2018) obtained rotation periods for these stars using Lomb-Scargle periodograms of their KELT light curves.We specifically selected stars with a measured period greater than 13.7 days to test recovery of long periods.To maximize the chances of recovering rotation periods, we used the TESS Science Processing Operations Center (SPOC) FFI simple aperture photometry (SAP) light curves (Caldwell et al. 2020).At the time of writing, only sectors 1-6 were available, but we trained our CNN using year-long (13-sector) light curves.However, the construction of our wavelet power spectrum used the same vertical (frequency) axis regardless of light curve length, so any length of light curve could be used without needing to retrain the neural network.
Upon visual inspection of the TESS -SPOC light curves, we noticed that many did not show obvious rotational modulation.Selecting only those light curves with unambiguous rotational modulation, we were left with 26 light curves with KELT rotation periods spanning 13.7 to 47 days.We generated wavelet power spectra and evaluated our neural network on them.Figure 7 shows our predictions for these 26 KELT stars.We subjected these points to the same cut in predicted fractional uncertainty as in Figure 6.Stars that made the cut are displayed in black, while stars whose predicted uncertainties were too large are shown in red.We successfully recovered stars with rotation periods longer than 13.7 and 27 days, even when TESS systematics were the dominant source of power in the wavelet power spectra.Furthermore, using the predicted fractional uncertainty as a quality cut successfully removed stars whose predictions were unreliable or wrong while ensuring the most reliable predictions remained in the sample.

MEarth and TOI-700
Only one long-period target was observed by MEarth in the TESS SCVZ: TIC 149423103.Newton et al. (2018) measured a rotation period of 111 days for this target from MEarth data.Using our neural network on the FFI data from TESS, we obtained P rot = 116 ± 48 days.While the predicted period was within 5% of the "true" period, the relatively large uncertainty (41%) means this target would fail our quality cut, and an ensemble period recovery attempt would miss it.
TOI-700 is another well-characterized star in the SCVZ.Using ASAS-SN data, Gilbert et al. (2020) estimated a precise rotation period of 54.0±0.8 days.Hedges et al. (2020) used a systematics-insensitive periodogram of its TESS light curve to obtain a period of 52.8 days.With our model we predicted a period of 59±53 days.Our period prediction was accurate to within 10%, but the large uncertainty would cause this target to be missed as well.

General period recovery and improvements
While we successfully recovered the rotation periods of these few hand-picked stars, robust recovery on larger, statistical scales will require more work and vetting.Our method allows us to see beyond the 13.7-day barrier, but only the stars with the largest amplitudes were reliably recovered-we are still limited by the TESS noise floor.TESS is less precise than Kepler at all magnitudes ( Vanderspek et al. 2018), so spot modulations require higher amplitudes to rise above the noise.In our simulated test set, we recovered rotation signals with some success down to amplitudes of a few parts per thousand, but our model was most successful at and above amplitudes of 1%.Both panels of Figure 5 show the 10thand 90th-percentile envelopes of the rotating stars from McQuillan et al. (2014).The bulk of Kepler 's rotating population falls between amplitudes of 1 and 10 parts per thousand and lie in a region where our recovery was less successful.These kinds of stars will be more difficult to recover with TESS, whatever the method.
Still, we believe improvements to our method will maximize what is recoverable from TESS.There are several ways to extend the predictability of our neural network to lower amplitudes and enhance the predictability at high amplitudes.The first and perhaps most useful improvement will come through light curve processing.Our processing pipeline followed the regression corrector documentation of Lightkurve Collaboration (2020) using a magnitude-dependent aperture threshold.In practice, a more carefully developed pipeline should be preferred.At the time of writing, the FFIs of sectors 1-6 have been reduced by both the TESS Asteroseismic We applied the same fractional uncertainty cut applied to the simulation recovery results; the 21 stars that made the cut are in black, while the 5 stars with unreliable period predictions are in red.We successfully recovered periods longer than 13.7 and even 27 days from real TESS light curves.Even when TESS aliases are the dominant sources of power in the wavelet transform, our neural network was able to recover the correct rotation period.
Science Operations Center (TASOC) pipeline and the TESS Science Processing Operations Center (SPOC) pipeline.Once sectors 7-13 are processed, the Southern Continuous Viewing Zone (SCVZ) will be complete, providing year-long light curves for hundreds of thousands of targets.These light curves will feature more careful systematics removal and should contain cleaner examples to use as "pure noise" light curves in our sample.We leave the use of these light curves to a future paper.
Another improvement may come with the inclusion of observation metadata.At least with TESS, certain systematic effects are specific to particular cameras, CCDs, or even locations on a CCD.Including camera number, CCD number, and x-and y-pixel coordinates in the training data set will allow neural networks to learn where to expect certain features and more easily ignore them in favor of astrophysical signals.
Improvements can be made to the neural network as well.In its current form, our model assumes that all input signals have rotation signatures, but not all real light curves display rotational modulation.In the future we may include a classification step like Lu et al. (2020) to determine which signals contain rotational modulation.Adding this classification step will allow the regressor to focus on signals with recoverable rotation, making for more efficient training.
It is important to note that our implementations of the conventional period recovery techniques still perform better than in reality (e.g., Canto Martins et al. 2020;Avallone et al. 2021, in prep., who were unsuccessful in recovering anything past 13.7 days).This indicates that, despite all our attempts to create as realistic a training set as possible, our simulations are not perfectly representative of real stars.It could be that our stitching routine fails to suppress long-period signals as the real light curves do.Another possibility is that our spot model, while tuned to the Sun, may not be representative of real spots on other stars.Whatever the reason, we have demonstrated the ability to recover periods even when the systematics that are present in our simulations make conventional techniques fail.
Even though our spot evolution simulations include latitudinal differential rotation, we were unable to recover differential rotation in this study.In some wavelet power spectra of our simulated light curves, the differential rotation is apparent as a slope in the frequency of maximum power versus time.When binning the power spectra to 64 × 64 pixels, the slope is more difficult to resolve.While increasing the resolution of the wavelet power spectrum images should enable the recovery of differential rotation, it will come at the expense of longer training time.We will investigate the recovery of differential rotation, activity levels, and spot properties in future work.
If we can see beyond the complicated systematics, TESS will deliver the largest set of rotation periods yet.McQuillan et al. (2014) obtained rotation periods for 34,000 stars in the Kepler field.The TESS continuous viewing zones combine to cover 900 square degrees around the ecliptic poles, representing about eight times the sky coverage of Kepler during its primary mission.We can therefore expect hundreds of thousands of new stars with rotation period estimates from the TESS CVZs, and perhaps more from lower ecliptic latitudes.Because of TESS 's lower precision compared to Kepler, the true number will likely be somewhat smaller, but the prospect of hundreds of thousands of new periods is worth continued refinements of this technique.We leave the application of this tool to the full CVZ samples to a future paper.

SUMMARY AND CONCLUSION
We used a convolutional neural network to recover rotation periods and uncertainties from simulated light curves with real TESS systematics.Despite the systematics, we successfully recovered periods even for targets whose periods were longer than the 13.7-day barrier encountered by conventional period recovery methods.In the half of the simulated test data with the smallest predicted fractional uncertainty, we recovered 10%-accurate periods for 46% of the sample, and 20%-accurate periods for 69% of the sample.We also found no significant misidentification of half-period aliases, unlike the Lomb-Scargle and wavelet methods.While periods were retrieved more successfully from higher-amplitude signals, the ability to predict uncertainties allows us to probe lower-amplitude rotation signals as well.
In future work, we plan to use this method to produce a catalog of rotation periods from TESS full-frame image light curves.We will also add output options to our neural network to predict latitudinal differential rotation and understand more of the properties used to produce the training set.With deep learning, we hope to maximize the output of TESS in spite of the frustrations that arise from its systematics.The ability to recover rotation periods, especially long periods, from TESS data will finally enable large studies of rotation across diverse populations of stars in the Galaxy if only the systematics can be learned.

Figure 1 .
Figure 1.An example of a butterpy simulation of spot evolution and light curve generation.This figure is available online as an interactive figure.The online figure has an interactive slider and play/pause buttons that allow the user to move the figure through time and see changes in the light curve as spots rotate into and out of view.

Figure 2 .
Figure 2. Top: Morlet wavelet transform (a) of a noiseless light curve, shown with the light curve (b) and global wavelet power spectrum (c).Center: plots for the same light curve convolved with TESS noise and re-stitched.The dotted curve marks the cone of influence, below which the power spectrum is susceptible to edge effects.Bottom: example of a binned wavelet power spectrum we used to train our neural network.Neural networks can learn to ignore the noise and pick out stellar signals.

Figure 3 .
Figure 3. Period recovery using Lomb-Scargle periodogram (LSP, top), global wavelet power spectrum (GWPS, middle), and autocorrelation function (ACF, bottom)."acc10" represents fraction of periods recovered to within 10% accuracy, while "acc20" is the recovery to within 20%.ACF has the highest overall success, but the recovery worsens significantly at periods longer than about 30 days.

Figure 4 .
Figure 4. Left: Period predictions by our convolutional neural network trained on wavelet transforms of noiseless light curves.Predicted periods have mean absolute percentage error of 14%, median absolute percentage error of 7%, acc10 of 61%, and acc20 of 81%.Right: Period predictions from noise-injected data, where recovery is significantly worse.Predicted periods have mean absolute percentage error of 246%, median absolute percentage error of 24%, acc10 of 28%, and acc20 of 45%.The horizontal band at 90 days represents targets where the model struggled to predict the period.In these cases, the prediction was biased toward the distribution median, or 90 days.

Figure 5 .
Figure 5. Neural network performance across the full simulation space of periods and amplitudes.In both panels, the white lines represent the 10th and 90th percentiles of the distributions from McQuillan et al. (2014), to gauge where stars from Kepler would fall.Left: Period recovery rate as a function of period and amplitude for the noise-injected data.The neural network performs better at higher amplitudes, where the rotation signal overpowers instrumental noise.Right: The same data, now colored by the neural-network-predicted fractional uncertainty in rotation period.The prediction is more certain for higher amplitudes.Furthermore, the prediction is most certain in the region with the highest recovery rate, indicating the predicted uncertainty is a reliable metric for period recovery without already knowing the true period.
Figure6.Period recovery for the half of the test set with the lowest predicted fractional uncertainty.Predicted periods have mean absolute percentage error of 57%, median absolute percentage error of 11%, acc10 of 46%, and acc20 of 69%.The predicted error cut removes the cluster of predicted periods at 90 days, giving credence to our cut to remove spurious period predictions.The cloud of objects with short true periods and predicted periods between 100 and 150 days have low fractional error because their predicted periods are large, but they account for only 4% of the objects remaining after the cut.

Figure 7 .
Figure7.Period recovery of stars in both TESS and KELT for which rotational modulation is apparent in the light curve.We applied the same fractional uncertainty cut applied to the simulation recovery results; the 21 stars that made the cut are in black, while the 5 stars with unreliable period predictions are in red.We successfully recovered periods longer than 13.7 and even 27 days from real TESS light curves.Even when TESS aliases are the dominant sources of power in the wavelet transform, our neural network was able to recover the correct rotation period.
The technical support and advanced computing resources from the University of Hawai'i Information Technology Services -Cyberinfrastructure are gratefully acknowledged.
. acknowledges support from NASA through an Astrophysics Data Analysis Program grant to Lowell Observatory (grant 80NSSC20K1001).

Table 1 .
Distribution of Simulation Input Parameters

Table 2 .
Convolutional Neural Network Architecture Layer Type Number of Filters Filter Size Stride Activation Output Size

Table 3 .
Metrics of Period Recovery on Simulated Light Curves