Complex Modulation of Rapidly Rotating Young M Dwarfs: Adding Pieces to the Puzzle

New sets of young M dwarfs with complex, sharp-peaked, and strictly periodic photometric modulations have recently been discovered with Kepler/K2 (scallop shells) and TESS (complex rotators). All are part of star-forming associations, are distinct from other variable stars, and likely belong to a unified class. Suggested hypotheses include star spots, accreting dust disks, co-rotating clouds of material, magnetically constrained material, spots and misaligned disks, and pulsations. Here, we provide a comprehensive overview and add new observational constraints with TESS and SPECULOOS Southern Observatory (SSO) photometry. We scrutinize all hypotheses from three new angles: (1) we investigate each scenario's occurrence rates via young star catalogs; (2) we study the features' longevity using over one year of combined data; and (3) we probe the expected color dependency with multi-color photometry. In this process, we also revisit the stellar parameters accounting for activity effects, study stellar flares as activity indicators over year-long time scales, and develop toy models to simulate typical morphologies. We rule out most hypotheses, and only (i) co-rotating material clouds and (ii) spots and misaligned disks remain feasible - with caveats. For (i), co-rotating dust might not be stable enough, while co-rotating gas alone likely cannot cause percentage-scale features; and (ii) would require misaligned disks around most young M dwarfs. We thus suggest a unified hypothesis, a superposition of large-amplitude spot modulations and sharp transits of co-rotating gas clouds. While the complex rotators' mystery remains, these new observations add valuable pieces to the puzzle going forward.

are often fast rotators, with rotational periods ranging from hours to one or two days. Their rotation is one of the drivers of their magnetic dynamos and thus stellar activity (e.g., Moffatt 1978;Parker 1979;Browning 2008). This can be observed in terms of activity indicators, such as hydrogen and calcium H & K emission lines, frequent and strong flaring activity, and significant star spot coverage (e.g., Benz & Güdel 2010;West et al. 2015;Newton et al. 2017;Wright et al. 2018;Günther et al. 2020).
In photometric observations, young M dwarfs with spots often show smooth, semi-sinusoidal rotational modulation with amplitudes of a few percent. Their patterns are rather 'simple', manifesting only a few peaks in a Fourier transform, even in the presence of multiple spots and differential rotation ( Fig. 1 first row). Thus, even the most extreme rotational modulations discovered so far can be described by just a handful of spots (e.g., Rappaport et al. 2014;Strassmeier et al. 2017).
One of the first phenomena clearly standing out from this norm were dipper and burster stars. These show abrupt dips or bursts of light in a quasi-periodic or stochastic manner (Alencar et al. 2010;Morales-Calderón et al. 2011;Cody et al. 2014;Ansdell et al. 2016, Fig. 1 second row), and were grouped by their photometric morphology into seven distinct classes 1 .
Shortly after, Stauffer et al. (2017Stauffer et al. ( , 2018) discovered yet again three new morphology classes in Kepler/K2 data 2 . As these three share common features, we refer to them collectively as scallop shells throughout this paper ( Fig. 1 third row). These scallop shells differ from the dippers/bursters in two substantial ways: (1) the objects discussed by Stauffer et al. (2018) are strictly periodic; and (2) they rotate much more rapidly, typically on timescales of 2 days, compared to the timescales of multiple days to weeks for the dippers/bursters.
Most recently, Zhan et al. (2019) discovered ten very similar objects in TESS Sectors 1 & 2,, dubbed complex rotators ( Fig. 1  of all light curves per TESS orbit). All complex rotators and scallop shells show rapid rotation, strict periodicity, and dozens of harmonics in their frequency spectra indicating their sharp light curve features. We thus argue that they are likely the same class of objects, and any differences are only due to Kepler's and TESS' observing cadences (30 min vs. 2 min). This makes them strictly distinct from 'normal' spotted stars (even those with differential rotation), which only show one or two peaks in their frequency spectra, and from dippers, which are far less periodic and morphological stable.

Hypotheses for complex modulations and their limitations
All dipper and burster classes were linked to the presence of dusty disks and a viewing-angle dependency, suggested by observations of strong infrared excess (accreting dust disks; Bodman et al. 2017, Fig. 2 first panel).
The scallop shells were first suggested to arise from a patchy torus of material clouds at the Keplerian co-rotation radius periodically transiting the star (corotating clouds; Stauffer et al. 2017Stauffer et al. , 2018 Fig. 2 second panel) 3 . Such material might be warm coronal gas, dust, or a mixture of both.
When studying the complex rotators, however, Zhan et al. (2019) suggested a new idea: spotted, rapidly rotating stars that host an inner dust disk at a few stellar radii, and show a spin-orbit misalignment between their rotation axis and the dust ring (spots and misaligned disk ; Fig. 2 third panel). Spots might then pass behind the dust disk and get (partially) occulted, leading to sudden increases in photometric brightness.
In particular, Zhan et al. (2019) presented the following counterarguments to the previous hypotheses: • Spots only: even the superposition of multiple cold and hot stellar spots leads to smooth variations and can only explain 1-2 peaks in the frequency spectrum (also Kővári & Bartus 1997;Stauffer et al. 2017). • Accreting dust disk : (i) complex rotators' stable periodicity is in stark contrast to the semi-periodic and stochastic nature of morphologies caused by accreting disks; (ii) the absence of significant infrared excess in their spectral energy distributions (SEDs) contradicts the presence of accreting disks; (iii) their roation periods are much shorter than those of dippers (also Stauffer et al. 2017). • Co-rotating clouds: (i) if the material is gas, it is challenging to explain the large amplitudes of the modu-  Figure 1. A collage of different morphology classes of young M dwarfs with complex photometric variability. First row: a typical complex rotator (TIC 206544316) with a period of 0.32 days. Its phase-folded light curve highlights the complex modulation, and its frequency spectrum shows dozens of harmonics down to time scales of minutes. Second row: a typical scallop shell (EPIC 204897050) with a period of 0.26 days. It is apparent that these targets share the same morphology class as the complex rotators, even though the data are constrained by the Nyquist limit for K2 (30 min. cadence). Third row: a typical spotted M dwarf with a period of 0.58 days (GJ 1243; see e.g. Davenport et al. 2015Davenport et al. , 2020Günther & Daylan 2021), showing a simple modulation and only two peaks in the frequency spectrum (two peaks rather than one due to differential rotation). Fourth row: a typical dipper star observed with TESS (RY Lup; Bredall et al. 2020), usually showing no strict periodicity, no presence of multiple clear harmonics, and a wider range of periods (multiple days to weeks). Not shown are the more standard variability classes, such as single spot modulation and binary stars. 4 lation; (ii) the material cannot be dust, as it cannot be stably confined at the required distances of several stellar radii because the magnetic field at those large distances from the surface would be too weak; (iii) any material (be it gas or dust) trapped in the magnetic field closer to the stellar surface could not reproduce the observations of sharp features with amplitudes of several percent.
1.3. This paper This paper focuses on the ten targets discovered by Zhan et al. (2019) to further scrutinize the plausibility of several hypotheses by three means: 1. investigating their occurrence rates, 2. studying the morphologies' stability and longevity over one (non-continuous) year, and 3. probing the feature's chromaticity.
We give an overview of all observations in Section 2, and revise the stellar parameters in Section 3. Then, we study stellar flares and other brightenings in Section 4. Next, we scrutinize all hypotheses with respect to occurrence rates, stability and longevity, and color dependency in Section 5. Finally, we discuss our findings and present our conclusion in Sections 6 and 7.

TESS photometry
The ten complex rotators were discovered in TESS short-cadence data from Sector 1 (2018-07-25 to 2018-08-22) and Sector 2 (2018-08-22 to 2018-09-10) and observed as part of the cool dwarf catalog (Stassun et al. 2018;Muirhead et al. 2018). Here, we also add new data taken over the full first year of operations (Tab. 1). Light curves were prepared with the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016), a descendant of the Kepler mission pipeline (Jenkins 2002;Jenkins et al. 2010;Jenkins 2017;Stumpe et al. 2014;Smith et al. 2012). We use the pre-search data conditioned simple aperture (PDC-SAP) light curves, which are detrended for instrumental systematics.

SPECULOOS Southern Observatory photometry
The SPECULOOS Southern Observatory (SSO; Gillon 2018;Burdanov et al. 2018;Delrez et al. 2018) is located at ESO's Paranal Observatory in Chile and is part of the SPECULOOS network. The facility consists of four robotic 1-m telescopes (Callisto, Europa, Ganymede, and Io), each equipped with a near-infraredsensitive CCD camera with a resolution of 0.35 arcsec per pixel. We observed four targets, TIC 201789285, TIC 206544316, TIC 332517282, and TIC 425933644, each simultaneously in at least two wavelength bands (g', r', i', and z' band filters) for an entire observing night  (Bodman et al. 2017), where an accreting dust disk seen from different observing angles could explain dipper and burster stars. Second panel: Co-rotating clouds hypothesis (Stauffer et al. 2017(Stauffer et al. , 2018, where a patchy torus of gas clouds at the Keplerian co-rotation orbit periodically blocks out stellar light and might cause scallop shell modulation. Third panel: Spots and misaligned disk hypothesis (Zhan et al. 2019), where spot occultations might explain the complex rotator's modulation. (Table A1). We extracted light curves using the SSO pipeline (Murray et al. 2020), which uses the casutools software (Irwin et al. 2004) for automated differential photometry and detrends for telluric water vapor.

ANU spectroscopy
We also reuse the spectroscopic observations taken by Zhan et al. (2019). The low-resolution spectra covered four of the systems (TIC 177309964, TIC 206544316, TIC 234295610, and TIC 425933644) using the Wide Field Spectrograph (WiFeS; Dopita et al. 2007) on the Australian National University (ANU) 2.3 m telescope at Siding Spring Observatory, Australia, on January 18 and 19, 2019. The observations cover the 5200-7000Å band with a resolving power of R = 7000 and were reduced following Bayliss et al. (2013). All spectra reveal strong Hα emission features with equivalent widths of 4-7Å, typical of rapidly rotating young M stars (see Section 3). No signs of a binary nature of these four objects were found.

Revisiting the binary TIC 289840928 and TIC 289840926
Zhan et al. (2019) found two prominent rotation periods for TIC 289840928, which is in a spatially resolved binary system with TIC 289840926 and both stars are blended in a single TESS pixel. In our reanalysis, the TESS pixel-level data revealed the primary (TIC 289840928, M4V, 3100 K) only has a smooth spot modulation with a period of 15.625 h, while the secondary (TIC 289840926, M6V, 2800 K) is the actual complex rotator with 2.400 h. We thus update all corresponding values here.

Ages and activity corrections
Zhan et al. (2019) estimated the ages of all ten complex rotators via their probabilistic membership of young stellar associations, using the banyan sigma software , with input from the TESS Input Catalog version 8 (TICv8; Stassun et al. 2018) and Gaia Data Release 2 (Gaia DR2 Collaboration et al. 2018). They found that all ten targets have a high probability of belonging to young associations (see Table 1). 4 We here conduct an independent estimate of their stellar age. Young M dwarfs are pre-main-sequence stars and as such have larger radii than their main-sequence counterparts of the same mass. Additionally, they have high levels of activity, leading to activity induced radius inflation and temperature suppression (Stassun et al. 2012). This could be due to strong chromospheric activity, presumably arising from rapid rotation. We can thus use the ANU spectra of four targets (Section 2.3) to explore if these stars have larger radii than expected for the main-sequence.
For example, TIC 234295610 shows an Hα equivalent width of 6.8Å in the ANU spectrum. According to the empirical relations from Stassun et al. (2012), this equivalent width (EW) predicts a radius inflation of 15.6% and a temperature suppression of 7.1%. Next, we performed a spectral energy distribution (SED) fit to better constrain the apparent radius to R apparent = 0.415 ± 0.029 R and the apparent temperature to T apparent eff = 3075 ± 100 K (see Fig. A2). Assuming that this radius and temperature represent the activityinflated radius and activity-suppressed temperature values, then in the absence of activity we would have R w/o activity ≈ 0.36 R and T w/o activity eff ≈ 3300 K. Comparing the corrected values with models for lowmass pre-main-sequence stars (Baraffe et al. 2015), we find that they are fully consistent with a star of mass of 0.20 M and age of 40 Myr. 5 This leads to a stellar mass that is only half of the mass of a main-sequence star with the same radius and effective temperature. We consider this a strong affirmation of the young age suspected from its young association membership.
We perform the same revision of stellar parameters for all systems. First, we perform SED fits to refine all their T apparent eff and R apparent (see Fig. A2). Next, for those with ANU spectra (TIC 177309964, TIC 206544316, TIC 234295610, and TIC 425933644), we measure their Hα EWs and use them to compute the temperature suppression and radius inflation factors, providing the corresponding values 'without activity'. For the other six targets, we provide provisional corrections assuming they have Hα EWs comparable to the average of the four measured ANU spectra. Fig. 3 and Table 1 summarize the revised stellar parameters and their isochrone matches. The 'corrected' T eff and radii place nine of the ten complex rotators in the range of ∼20 or ∼50 Myr, consistent with the respective ages of their young associations. Only TIC 332517282 falls closer to the ∼150 Myr isochrone, consistent with its membership in AB Doradus (150 Myr), making it our single 'oldest' star.  Figure 3. Revised stellar parameters of the ten complex rotators after correcting for activity-induced radius inflation and temperature suppression. Solid curves are isochrones from the Baraffe et al. (2015) models for ages of 20, 50, and 150 Myr, with gray circles marking masses at 0.1 M increments. Error bars show apparent effective temperatures (T eff ) and radii measured from SED analyses. Arrows show where the stars would fall if they were inactive, i.e., after they are corrected for activity. Filled blue circles mark the four targets with Hα measurements from ANU (Section 2.3), whereas empty blue circles mark the remaining six targets for which we provide estimates only. The 'corrected' values place most of the stars near isochrones consistent with the respective ages estimated from their young association memberships, suggesting inferred masses of 0.1-0.3 M (Table 1).

FLARES AND SUDDEN BRIGHTENINGS
The complex rotators and scallop shells show frequent and large-amplitude flaring, along with other sudden brightenings whose shapes are distinct from usual M dwarf flare profiles (Stauffer et al. 2017(Stauffer et al. , 2018Zhan et al. 2019). In particular, brightenings of the entire modulation often appeared right after strong flares, sometimes even followed by changes in the overall morphology, underlining strong magnetic activity.
It is still disputed whether flares on stars other than our Sun correlate with the rotational phase and are linked to localized clusters of spots on the stellar surface. Many previous studies found that superflares were distributed randomly uniform over the rotational phase for main sequence dwarfs (Roettenbacher & Vida 2018;Doyle et al. 2018Doyle et al. , 2019Doyle et al. , 2020, and young stars (Vida et al. 2016;Feinstein et al. 2020b). Doyle et al. (2018) reason that depending on the viewing geometry, polar spots could be seen at all phases, and their interaction with emerging active regions can thus cause continuously visible flaring. However, other studies found flares to be more prominent at certain rotational phases, and hence potentially bound to the locations of star spots for the Sun (Zhang et al. 2008), Sun-like stars (Notsu et al. 2013), and the smallest flares on main sequence dwarfs (Roettenbacher & Vida 2018). Hence, a unifying idea is that superflares occur over the entire surface while small flares are tied to spots.
Here, we utilize the extended coverage by TESS (up to one year for TIC 177309964) to study whether complex rotators' flares and other brightenings correlate with the phase of the modulation. We searched the light curves in two ways. In the first approach, we ran the stella software, a convolutional neural network developed for probabilistic flare detection in TESS 2 min. cadence data (Feinstein et al. 2020b,a). As the algorithm was trained on a large sample of stars with smooth spot modulation, many of the initial flare candidates were misidentified (often spikes of the complex rotators). By visually vetting, we then selected reliable criteria of a probability ≥0.9 and amplitude ≥5%, and identify a confirmed sample of at least 67 flares on TIC 177309964. In the second approach, we independently inspected the entire light curve by eye, and found a total of ∼70 confirmed flares, agreeing with the machine learning results.
We find that flares on the complex rotators are distributed randomly uniform in phase, showing no clear dependency of their location or amplitude (see example of TIC 177309964 in Fig. 4). There is also no clear correlation between flare amplitudes and the one year time span. It is peculiar that the three largest flares (amplitudes of 1.5 to 4 in normalized flux) all occur within two few weeks from one another, but the sample size is too low to rule out mere coincidence. Most flares are described by the same profiles as their main sequence counterparts, suggesting similar origins and processes driving them. We also observed somewhat more complicated 'outbursts' of flares, which again resemble those of main sequence M dwarfs; these can be explained as superpositions of multiple flare events (e.g. Günther et al. 2020).
Finally, we also find that some sudden brightenings do not resemble typical M dwarf flare profiles (Fig. 5); instead, they seem like amplified versions of the complex periodic morphology. This was also pointed in Zhan et al. (2019) and Stauffer et al. (2017), often (but not always) following superflares. On TIC 177309964, these alterations also occur without any preceding flare observation. It is possible that they were triggered by a superflare that was (i) not visible in the visual, but would have been classified as such in the UV or X-ray spectrum (e.g. Wolter et al. 2007;Loyd et al. 2020), or (ii) was located outside of the visible hemisphere.

TESTING THE HYPOTHESES
The limitations of previous hypotheses leave us with two remaining possibilities for the complex rotators and scallop shells: (i) the idea of a patchy torus of clouds of gas at the Keplerian co-rotation radius (co-rotating clouds hypothesis) and (ii) the idea of spots being periodically occulted behind a spin-orbit-misaligned dust disk (spots and misaligned disk hypothesis). In the following analyses, we hence focus on these two cases.

Occurrence rates
We can estimate the expected yield of complex rotators in TESS Sectors 1 & 2 from the co-rotating clouds and spots and misaligned disk hypotheses (see Fig. 2), respectively, as and N SMD comp.rot. = N yMrr · P spots · P disk · P misal · P geom . (2) Here, N yMrr is the number of young M dwarfs with rapid rotations ( 2 days  probability that a given star has a cloud of dust/gas orbiting it at the Keplerian co-rotation radius. P spots is the probability that a given star has at least one large spot. P disk is the probability that these stars have a dust disk orbiting them at a few stellar radii. P misal is the probability that a given star shows a spin-orbit misalignment between the stellar rotation axis and the dust disk. Finally, P geom is the geometric probability that an existing structure (cloud or disk) falls in the line of sight between the star and the observer. To first order, we can estimate N yMrr as the number of stars in known open clusters and associations that Complex Rotators II 9 have effective temperatures below 3900 K and radii below 0.6 R . For this, we cross-match the TESS short cadence target lists 6 of Sectors 1 & 2 with three young star catalogs: 1. a catalog by Feinstein et al. (2020b), which was assembled through a combination of searching the TESS Guest Investigator proposals and the data from Faherty et al. (2018). All targets were vetted with the banyan sigma software ) and the ones with ≥ 50% probability and available TESS 2 min data were included. 2. a catalog by Bouma et al. (2019), collecting targets from numerous literature lists, including members of open clusters, moving groups and young associations. 3. a catalog we created by matching all TESS shortcadence targets with the banyan sigma software . The algorithm uses a Bayesian model to predict whether a given target is part of one of 27 young associations within 150 parsec. The target is identified by its coordinates, proper motion, parallax, distance, and radial velocity measurements, which we retrieve from TICv8 and Gaia DR2. The results of our cross-matches are shown in Figure 6. We find a total of 290 young M dwarfs from open clusters and young associations. Notably, their real number might be even higher, as new clusters and associations are still being discovered. For example, Gagné et al. (2020)  To estimate N yMrr from this, we investigate the rotation periods of those 290 M dwarfs using all available TESS data from Sectors 1 through 13. We flag a target as a 'young M dwarf with rapid rotation' if we find a rotation period below 2 days using Lomb-Scargle periodograms (Lomb 1976;Scargle 1982). We measure rotation periods for 269 out of the 290 targets, with the remaining 21 targets not showing measurable photometric variability (false alarm probability > 0.01). This might just be due to low signal-to-noise, or the stars might simply be seen pole-on; both leaves open that they could actually be rapidly rotating. All measured rotation periods and effective temperatures are shown in Fig. 7, and we conclude that N yMrr 175. Out of these, 175 targets show rotation periods shorter than 2 days (blue shaded area; NyMrr: number of young M dwarfs with rapid rotations), comparable to the ten known complex rotators from this sample (red squares).

How many complex rotators should we expect?
We here put the finding of N yMrr 175 into the perspectives of Eqs. 1 and 2.
Co-rotating clouds hypothesis -For this hypothesis, we also have to account for (i) P clouds , the probability of clouds of dust/gas being present at the Keplerian corotation radius, R cr , as well as (ii) P CC geom the geometric alignment probability of an edge-on alignment to the observer. We can estimate 7 P CC geom ≈ (R /R cr ), with R cr derived as: with the rotation period P rot , gravitational constant G and stellar mass M . For our targets, this yields R cr ≈ 1 − 6 R (Table 1), which leads to P CC geom ≈ (R /R cr ) ≈ 0.2 − 1. Putting all pieces together, we can estimate from Eq. 1: Given that we found 10 complex rotators in Sectors 1 & 2, this would only require 30% of all rapidly rotating young M dwarf to have material clouds trapped at their Keplerian co-rotation radius.
Spots and misaligned disk hypothesis -For this hypothesis, we already searched all known young stars in TESS Sectors 1 & 2 for photometric rotation periods and signs of spots, and found that, on average, P spots ∼ 1 (see above). We estimate the geometric probability P SMD geom ≈ (R /d), where R is the stellar radius and d ≈ 5 − 15 R is the outer edge (or gap) of a possible disk (Zhan et al. 2019). This leads to P SMD geom ≈ 0.07 − 0.2 and: Considering our 10 complex rotators in Sectors 1 & 2, this would imply that P disk · P misal is on the order of 0.3 − 1. Hence, for this hypothesis to hold true, a large fraction of 30% of rapidly rotating young M dwarfs would need to have an inner dust disk with a slight spinorbit misalignment to their rotation axis.

Time dependency
The photometric modulations of the complex rotators appear stable over timescales of at least one year ( Fig. 8  and Fig. 9). The complex rotator TIC 177309964 fell into TESS' continuous viewing zone and was observed for the consecutive Sectors 1-13, spanning an entire year of data from July 2018 until July 2019 (Fig. 8).
We can also see this stability on the examples of TIC 201789285, TIC 206544316, TIC 332517282, and TIC 425933644 (Fig. 9) when combining TESS and SSO photometry. The original TESS data were taken in Aug-Oct 2018, while the SSO observations were taken in Nov/Dec 2019, over one year later. Despite the large time span, the modulation profiles still follow the same pattern and periodicity. All major features remain the same, while only some minor evolution of the morphology is evident in the SSO light curves compared to TESS. In the case of TIC 201789285, it appears that a minor feature has increased in amplitude (near BJD TDB 2,458,795.70), while another has decreased in amplitude (near BJD TDB 2,458,795.75).
The stability and longevity of these morphologies is extraordinary. Considering the spots and misaligned disk hypothesis, this is well compatible with the life times of dust disks and the persistence of stellar spots on young M dwarfs. While stellar spots on most stars only last for a few weeks (e.g., on the Sun they last only for 3-6 rotations Gaizauskas et al. 1983), we have examples like GJ 1243, which had remarkably constant spot modulation over 10 years observed with Kepler and TESS (Davenport et al. 2020). Notably, GJ 1243 is a member of a young association with an age of around 30-50 Myr and, as such, quite comparable to our complex rotators. Furthermore, a 200-day photometric monitoring campaign of the open cluster Blanco 1 (∼115 Myr) with the Next Generation Transit Survey (NGTS) suggests that most young M dwarfs display generally stable spot modulation patterns over this baseline, while F, G and early-K dwarfs show moderate-to-significant evolution in their light curve morphologies (Gillen et al. 2020).
In contrast, considering the co-rotating clouds hypothesis, such a stability and longevity would seem surprising. The idea does suggest a rather fine-tuning problem with clouds being confined strictly at the co-rotation radius. Any separating clouds would slowly drift away and slightly alter their orbital period. Even for small drifts, a year-long time-span might lead to noticeable changes in the morphologies. This would lead to a certain amount of material away from co-rotation at any given time, which would blur out the strictly periodic signals over year-long time spans.
TIC 177309964 Period = 10.88 h Figure 8. Example of TIC 177309964, illustrating the stability of major photometric features over year-long timescales (TESS orbits from Sectors 1-13, from July 2018 until July 2019). Each orbit is phase-folded to show two periods of the modulation and vertically offset for clarity. Flares have been clipped before plotting. The major features of the modulation remain stable over the full year, while minor features appear and disappear over a few weeks (dashed box).

Color dependency
We obtained a total of nine telescope nights worth of SSO observations (Section 2.2) to capture four of the complex rotators in simultaneous multi-color bandpasses. We compare all these light curves with phasefolded TESS observations (taken one year earlier) in Fig. 9. Evidently, the sharp-peaked features are more prominent in bluer bandpasses, and less expressed in the reddest bandpasses. This matches the expectations from both the co-rotating clouds and spots and misaligned disk hypotheses: (i) for the co-rotating clouds hypothesis, the material's extinction would have to be stronger in the blue, leading to deeper features. (ii) for the spots and misaligned disk hypothesis, the contrast between the stellar surface (∼3000 K) and a cool spot is stronger in bluer wavelengths. The disk material could be a gray absorber or could have a color dependency, which would add a secondary effect.
Combining TESS data with SSO r, i and g band observations, our total data span more than one year. We find that the same features are still present in the data at the predicted phases (accounting for uncertainties in the period estimation). There is no doubt that the modulation is still the same, and thus stable and long-lived over year-long time spans. 6. DISCUSSION

Could spots and pulsations be another hypothesis?
Pulsations on low mass stars have long been theoretically explored and predicted (e.g., Gabriel 1964;Noels et al. 1974;Palla & Baraffe 2005;Rodríguez-López et al. 2012; Rodríguez-López 2019) but, despite observational efforts, have so far proven elusive. In theory, fully non-adiabatic models of M dwarfs suggest that they could excite (i) radial modes, (ii) low-order, lowdegree non-radial modes, and (iii) solar-like oscillations (Rodríguez-López et al. 2012. This requires the models to be completely convective or have large convective envelopes. This would match the spectral types of all complex rotators and scallop shells, placing the stars beyond the fully convective limit. Periods of the pulsations are predicted to range from 20 min to 3 h, again agreeing with the typical time scales we see for the sharp-peaked features of our targets. However, the amplitudes of M dwarf pulsations are expected to be in the range of 1 ppm to 1 ppt, while our targets typically show amplitudes of several percent. The only bypass to this caveat could be a superposition of the effects from spots (creating the smooth, large-amplitude modulations) and synchronous pulsations (creating the sharp-peaked, lower-amplitude features). Yet, this would require spots and pulsations to be synchronized, which seems implausible. Rodríguez-López et al. (2014) identified the theoretical instability strips of M dwarfs, which resulted in two islands of 'instability' in the parameter space. Stars falling into one of these islands would, in theory, be capable of showing pulsations. With the revised and activity-corrected stellar parameters (Section 3), some of the stars fall near the lower island, yet remain outside of it. Again, this makes pulsations seem implausible.

The pros and cons of various hypotheses: a summary
We here briefly summarize the pros and cons of the various hypotheses introduced and scrutinized throughout this paper (see Sections 1 and 6.1 for short overviews). Fig. 10 additionally contrasts all observations with all hypotheses, highlighting which aspects can and cannot be explained through a given hypothesis.
Spots only -Pro: Spot modulations can be strictly periodic and stable over many years, even in the presence of differential rotation (e.g., Davenport et al. 2020). The temperature difference between the surface and the spots causes a color dependence, and spots would not cause any infrared excess. Most young stars are spotted and are often accompanied by strong signs of mag- Figure 9.
Comparison of multi-color light curves of four of the ten complex rotators: TIC 201789285, TIC 206544316, TIC 332517282, and TIC 425933644. For each target, SSO observations were taken simultaneously in at least two of the g', r', i', and z' bandpasses (shown as green, orange, red and dark-red curves, respectively). For TESS observations, flares have been removed, and light curves are averaged over all available Sectors, binned in 5 min intervals, and slightly shifted in phase to correct for imprecision in the period measurement (shown as dark-gray curves). The lower panel compares the normalized transmission functions of all respective bandpasses. There is a clear color dependency of the light curve features visible in the simultaneous SSO observations, with features being much more prominent in bluer bandpasses. Additionally, the general shape and largest features in the SSO light curves are still comparable to the TESS light curve, even though the SSO data were taken about 1 year later. This suggests that the overall mechanism causing these patterns also causes a color-dependency (e.g., spots, non-gray dust, or pulsations), and that it is stable over long times.
netic activity, such as flaring. Con: Spots alone cannot explain the sharp-peaked features (Zhan et al. 2019). However, spots could still play a major role in combination with other factors (e.g, pulsation or circumstellar material; Sections 6.1 and 6.3).
Accreting dust disk -Pro: Accreting dust disks can lead to the morphologies for dipper/burster stars and occur frequently enough. They might show color dependency depending on the absorption and scattering properties of the material. Con: Accretion is a rather stochastic process, and thus neither strictly periodic nor stable. The dippers and bursters also show strong infrared excess due to the large disks.
Co-rotating clouds of material -Pro: Clouds of material at the Keplerian co-rotation radius could qualitatively explain sharp features and strict periodicity (Stauffer et al. 2017(Stauffer et al. , 2018. Depending on the material, a color dependency is possible, and small enough clouds would cause no infrared excess. As young stars are often surrounded by material, they could also occur at high enough rates. Con: If the material is gas, the absorption would likely not be able to explain percentage-scale amplitudes (Zhan et al. 2019). If the material is dust, these clouds are likely not stable at the required distances (d/R 3; Zhan et al. 2019). Another challenge might be the stability and longevity of the morphology over year-long time spans, i.e., over hundreds of orbital periods. With some parts of the clouds slowly drifting away from co-rotation, the signals would blur out and evolve, which does not seem to be the case.
Material trapped near the surface -Pro: Material trapped in the magnetic field and bound to the stellar rotation would remain strictly periodic, and could survive over many years near the stellar surface (d/R ∼ 1; Zhan et al. 2019). Depending on the material's properties, a color dependence is possible, and in small amounts it might not cause any infrared excess. Con: Any material that close in cannot explain the sharppeaked, percentage-scale amplitudes of the modulation but would instead produce a rather smooth variation similar to spots; this can only be explained by material at larger distances (d/R 3; Zhan et al. 2019).
Spots and spin-orbit-misaligned dust disks -Pro: The patterns can be strictly periodic and stable over many years. Spots induce a color dependency, and disk material might add to this effect. There are enough young and rapidly rotating M dwarfs in the sample to explain the high occurrence rates (with caveats, see below). Lastly, the spots are a sign of magnetic activity, and agree with the frequent flaring found on the complex rota-tors and scallop shells. Con: the scenario would require most young M dwarfs to have close-in dust disks with spin-orbit misalignments. There is no obvious formation mechanism that would explain this behavior. Also, the one M dwarf for which we have disk and rotation measurements, Au Mic, does appear co-planar. However, the misalignment does not need not to be very large. A 10 • obliquity between the spin and magnetic axes of T Tauri stars is reasonable, based on Zeeman studies and recent work by McGinnis et al. (2020). If the disks are confined by the magnetic field, this slight misalignment could already be enough to mitigate this caveat and cause the observed morphologies. A potential driver for such misalignment might be perturbations from nearby passing stars, either dynamically or through radiation pressure (e.g. Rosotti et al. 2014).
Spots and pulsations -Pro: In theory (qualitatively), a superposition of spots and pulsations could lead to a smooth large-amplitude trend (due to spots) superposed by a sharp-peaked small-amplitude pattern (due to pulsations). Both features can be stable over long times, show color dependencies, a lack of infrared excess, and frequent flaring due to their activity. Con: This would require a strict synchronization between rotation and pulsation, which seems implausible. Further, the stars do not lie within any of the theoretical instability islands.
6.3. Towards a unified hypothesis: could it be spots and co-rotating clouds of material?
So far, none of the hypotheses stand out as a definite answer, and each come with limitations. The two most promising hypotheses have their own caveats. On the one hand, the co-rotating clouds hypothesis, with clouds of gas at the Keplerian co-rotation radius could remain stable and allow the required viewing angle, but gas likely cannot explain the large, percentage-scale absorption features. On the other hand, the spots and misaligned disk hypothesis could tick all boxes, but has the implicit requirement that a large fraction of rapidly rotating young M dwarfs must have misaligned disks. That said, if the disks are confined by the magnetic field, a small misalignment of ∼10 • seems not to be uncommon (McGinnis et al. 2020) and could have been induced by nearby passing stars (e.g. Rosotti et al. 2014).
Rapidly rotating young M dwarfs are known to be magnetically active and generally show high spot coverage rates. Spots alone were ruled out easily, but what was not explored so far is: how much can spots explain?
The superposition of smooth, large-amplitude spot modulations and sharp, sudden features from transits of co-rotating clouds of gas could represent a unified hy- pothesis (Spots and co-rotating clouds). This explanation could expand the co-rotating clouds hypothesis by requiring only a minimal amount of circumstellar material to cause the overall morphology, making gas clouds a more plausible candidate. It would also mitigate the occurrence rate caveats which challenge the spots and misaligned disk hypothesis.

Toy models
We developed simplified forward models for the three most promising ideas, the (i) co-rotating clouds, (ii) spots and misaligned disk, and (iii) spots and co-rotating clouds hypotheses. We took the TESS light curve of TIC 201789285 as an example, and tried to imitate its morphology as closely as possible while keeping the models simple, using a hybrid of statistical inference and manual parameter selection. We model the star as a fine grid in spherical coordinates. The rotation axis of the star is left as a free parameter, and a quadratic limb darkening effect is applied to each cell based on its orientation relative to the observer.
For the spots and misaligned disk hypothesis, we model spots with four parameters: two angles describing the location on the star, its size, and its temperature. The disk is parametrized by its inner and outer radius, inclination, and opacity. The observed flux is computed for each grid cell by integrating the Planck function across the TESS bandpass, accounting for geometric effects, spots and the disk. We then try to mimic the TESS light curve of TIC 201789285 by choosing a simple model with three cold spots and a fully opaque disk, optimizing their parameters using nested sampling via dynesty (Speagle 2020).
For the co-rotating clouds and spots and co-rotating clouds hypotheses, we model the torus of clouds as a series of orbiting spheres whose orbital period matches the rotational period of the star. For this, we used the processing package (https://processing.org/) to draw the model and count the flux in each pixel. We then manually evaluated different scenarios of spots and sparse to dense tori of clouds.
We find that all three toy models can replicate the typical morphology of complex rotators (Fig. A3). Additionally, spots as drivers for an underlying largeamplitude modulation can explain a large portion of the signal, easing the constraints on circumstellar material. 6.5. What can TESS short-cadence do for us?
The TESS complex rotators were all found in shortcadence (2 min) observations. In contrast, the K2 scallop shells were discovered using 30 min cadence. However, the longer cadence means that the data are limited by the Nyquist frequency, and numerous harmonics will be missed in the frequency spectrum (see Fig. 1). We tested whether we would have discovered the complex rotators in the same way from TESS long-cadence (30 min) observations, by extracting light curves directly from the full frame images (FFI). For this, we used a TESS FFI photometry pipeline based on the one developed by Pál (2012). It is designed to extract photometry of faint stars in crowded fields (TESS-mag>15) by combining difference and aperture photometry. We selected nearby stars from the Gaia DR2 catalog that could contaminate the target's light curve and applied a principle component analysis (PCA) detrending.
We find that the complex rotators morphologies are most apparent in the 2 min cadence data, and that many might have been missed due to their sharp-peaked features being blurred out in 30 min cadence data (Fig. 11). The 2 min data will thus be the best source to search for more complex rotators throughout Sectors 1-23. With the start of the TESS extended mission, 20 s cadence data will be enabled. Additionally, both the original TESS sectors as well as the K2 fields will be revisited. Re-observing complex rotators and scallop shells with such short cadence could greatly increase our sensitivity to the sharpest features.

Implications for exoplanet systems
We know thanks to Kepler and other missions that, in average, each early-to mid-M dwarf has at least one small exoplanet (Dressing & Charbonneau 2015). All complex rotators and scallop shells are young mid-M dwarfs (most around 10 -45 Myr), raising the question of how the effects causing their morphology might impact young, recently formed exoplanets. After all, for the co-rotating clouds, spots and misaligned disk, and spots and co-rotating clouds hypotheses, occurrence rates suggest that most young M dwarfs would go through the same processes, we just cannot see their morphologies (Section 5.1).
At these ages, most processes driving planet formation have likely been concluded. Terrestrial planet formation is thought to be completed after at most 10-30 Myr, regardless of the driving processes (e.g., Chambers 2010). Particularly, gas giant planets around mid-to late M dwarfs are rare; even if they formed around any complex rotators, they require a substantial amount of gas in the proto-planetary disk to be present at the late stage of their formation, yielding formation time spans of less than 10 Myr (the maximal lifetime of gas discs; e.g., D'Angelo et al. 2010).
The material in question for causing the complex morphologies, however, is likely much closer to the star (near 3-10 stellar radii) than any forming or migrating exoplanets. Exoplanet systems hosted by mid-M dwarfs are widely studied, and if the effects at play for complex rotators are indeed ubiquitous at early ages, they seem to not have caused a strong effect on their planets.
If the spots and misaligned disk hypothesis were correct, more than half of all young M dwarfs would have a slight spin-orbit misalignment between their rotation axis and remnant dust disk. An interesting question is whether this would imply the proto-planetary disk to also be misaligned. However, there are currently no strong signs that a large fraction of exoplanets forming around M dwarfs have spin-orbit misalignments. Limited data suggest there to be a wide variety. For example, Au Mic b is proven to be aligned (e.g, Hirano et al. 2020a), the TRAPPIST-1 system might be slightly misaligned (Hirano et al. 2020b) and GJ 436 b is inclined (Knutson et al. 2011;Bourrier et al. 2017).

CONCLUSION
Recently, at least one new distinct morphology class of young stars has been discovered in white-light photometry from Kepler/K2 and TESS (Stauffer et al. 2016(Stauffer et al. , 2017Zhan et al. 2019), adding to the seven classes established by Cody et al. (2014) (which include dippers and bursters). Here, we added three puzzle pieces to unveil the physical nature of these complex rotators and probe whether several hypotheses could hold true given these new observational constraints.
The tested hypotheses include spots only, accreting dust disks, co-rotating clouds of material, magnetically constrained material, spots and a spin-orbit-misaligned disk, and spots and pulsations. We particularly focused on the co-rotating clouds and spots and misaligned disk hypotheses, as others are ruled out more easily.
First, we investigated if their occurrence rates make sense in the light of the total number of rapidly rotating young M dwarfs in the given field of view. We find that TESS Sectors 1 & 2 harbor at least 175 young M dwarfs with rotation periods below 2 days, rendering the finding of 10 complex rotators plausible. However, for both hypotheses this comes with a caveat. For the co-rotating clouds hypothesis to work, this would mean that almost every such star must have clouds of dust/gas trapped at the Keplerian co-rotation radius. For the spots and misaligned disk hypothesis to work, this would imply that a large fraction of these stars must have an inner disk and show a spin-orbit-misalignment. If the latter held true, it could have consequences for exoplanet systems around mid-to-late M dwarfs, which might have formed under the same conditions. Second, we studied the longevity of these features over one year, and find that they remain remarkably stable over these time spans. While the major features of the complex rotators remain unchanged, we find evidence for additional small features building up and decaying over a few weeks. In the Co-rotating Clouds hypothesis, this would imply subtle changes in the dust/gas cloud structures. In the Spots and Misaligned Disks hypothesis, this can very likely be caused by smaller spots appearing and disappearing, major spots changing size, or spots wandering along the surface.
Third, we probe the color-dependency of the complex rotators photometric features. We indeed find the expected behavior predicted by both the co-rotating clouds and spots and misaligned disk hypotheses. The features are more pronounced in bluer wavelengths, which could be explained by either chromatic absorption by the circumstellar material or the smaller spot-to-star brightness contrast in the red / infrared.
All new clues to the case -occurrence rates, longevity and color-dependency -could in principle match any of the hypotheses shown in Fig. 2. It is well possible that the truth will lie somewhere between these hypotheses. Rapidly rotating young M dwarfs are known to be magnetically active, so the final answer will likely have contributions from both spots and circumstellar material, leading to their complex photometric morphologies.  Figure A1. All available TESS light curves of young M dwarfs with complex photometric variability, spanning multiple Sectors of TESS observations. Each light curve shows the data of one TESS orbit (14 days of observations) phase-folded onto the rotational period of each star, with, e.g., "S1 i" denoting the first orbit of Sector 1 and "S13 ii" denoting the second orbit of Sector 13. Flares have been removed prior to phase-folding and plotting. All complex rotators' features show remarkable stability and longevity over many weeks to a year (at least). See also Fig. 8 for a closer view of TIC 17730996.  Figure A2. All spectral energy distribution (SED) fits for the ten complex rotators. The results are used to determine the stars' apparent effective temperatures and radii. The near-UV and blue optical excesses for some targets are consistent with chromospheric activity, another indicator of strong magnetic activity for these stars. Notably, no infrared excesses are seen. Flux Spots and co-rotating clouds Figure A3. Toy models of a complex rotator derived from three different hypotheses, aiming to replicate the morphology of TIC 201789285.  Table A1. Summary of all SSO observations for four of the ten complex rotators.