The Benchmark M Dwarf Eclipsing Binary CM Draconis With TESS: Spots, Flares and Ultra-Precise Parameters

A gold standard for the study of M dwarfs is the eclipsing binary CM Draconis. It is rare because it is bright ($J_{\rm mag}=8.5$) and contains twin fully convective stars on an almost perfectly edge-on orbit. Both masses and radii were previously measured to better than $1\%$ precision, amongst the best known. We use 15 sectors of data from the Transiting Exoplanet Survey Satellite (TESS) to show that CM Draconis is the gift that keeps on giving. Our paper has three main components. First, we present updated parameters, with radii and masses constrained to previously unheard of precisions of $\approx 0.06\%$ and $\approx 0.12\%$, respectively. Second, we discover strong and variable spot modulation, suggestive of spot clustering and an activity cycle on the order of $\approx 4$ years. Third, we discover 163 flares. We find a relationship between the spot modulation and flare rate, with flares more likely to occur when the stars appear brighter. This may be due to a positive correlation between flares and the occurrence of bright spots (plages). The flare rate is surprisingly not reduced during eclipse, but one flare may show evidence of being occulted. We suggest the flares may be preferentially polar, which has positive implications for the habitability of planets orbiting M dwarfs.


INTRODUCTION
M dwarfs are the most common spectral type and very popular targets for exoplanet surveys, owing to a shorter period "habitable zone" that is easier to observe, and a high abundance of terrestrial planets (Dressing & Charbonneau 2015).Significant time and resources have been dedicated to this field, including new spectrographs sensitive to redder wavelengths, transit surveys including Mearth, TRAPPIST, SPECULOOS and even TESS, as well as significant JWST time for transmission spectroscopy.However, there remain two major challenges with M dwarfs.
First, the precision of our exoplanet parameters is a function of the precision in the host star parameters.For M dwarfs, empirical stellar mass-radius relations are poorly constrained (Torres et al. 2010) and often do not match theoretical models (e.g.radius "inflation", Parsons et al. 2018).This stems from a paucity of well-characterised M dwarfs in eclipsing binaries, which are the traditional calibrators (Parsons et al. 2018;Gill et al. 2019;Swayne et al. 2021;Sebastian et al. 2022;Maxted et al. 2023).
A second challenge with M dwarfs is heightened stellar activity, evidenced by frequent flares and a high spot coverage.This impedes planet detection with transits (e.g.Kipping et al. 2017;Feliz et al. 2019;Gilbert et al. 2021Gilbert et al. , 2022a) ) and radial velocities (RVs) (Aigrain et al. 2016).Activity also hinders characterisation of planets and their host star, and may be a cause of the discrepancy between models and observations of low-mass stars (López- Morales & Ribas 2005;Lubin et al. 2017).Even if we were to discover an "Earth-like" planet, its habitability may be questionable if the planet is frequently bombarded by flares and associated coronal mass ejections (CMEs) (Ranjan et al. 2017;Tilley et al. 2019;France et al. 2020;Bogner et al. 2022).The largest well-recorded solar storm (flare + CME) for the Sun was the Carrington Event of 1859.It released ≈ 10 32 erg, causing aurorae to be seen globally.If it occurred today it would devastate our power grids.Carrington-esque events can be a daily occurrence on M dwarfs.More knowledge of M-dwarf activity is needed to assess the viability of potentially habitable worlds.Fortunately, the availability of long-term photometry from observatories such as Kepler, TESS and the All Sky Automated Survey for SuperNovae (ASAS-SN) has revolutionised studies of flares and spots (Davenport 2016;Schmidt et al. 2019;Günther et al. 2020;Rodríguez Martínez et al. 2020;Feinstein et al. 2022;Mendoza et al. 2022).
In this paper we use 15 sectors of TESS photometry to study CM Draconis, a double-lined, spectroscopic eclipsing binary with twin M dwarfs, introduced in detail in Sect. 2. There are two fundamental aspects to our paper, related to the aforementioned two challenges with M dwarfs.First, we use the TESS photometry and archival radial velocities (Sect.3) to provide an order of magnitude improved fit to the eclipses (Sect.4), which we compare to stellar models.Second, we identify and characterise CM Dra's activity through flares (Sect.5) and star spots (Sect.6).We quantify the flare rate and place constraints on both the longitudinal and latitudinal distribution of the stellar activity.Through this, we quantify the connection between spots and flares and estimate the activity cycle, akin to the 11-year Solar cycle.All of our analysis is contained in Sect.7, before concluding in Sect.8.

THE CM DRACONIS SYSTEM
CM Draconis (CM Dra henceforth) is a benchmark system for the study of M dwarfs.It is a 1.26 day eclipsing binary of near twin fully convective dM4.5 stars.The two stars are remarkably well aligned to our line of sight with ( = 89.8• ,  = 0.12).With near identical radii ( A = 0.251 and  B = 0.238) this means that both primary and eclipses result in a near perfect occultation (50% drop in total flux).With near equal masses ( A = 0.225 and  B = 0.211) CM Dra is a double-lined spectroscopic binary, meaning that model-independent masses can be derived for both components.Indeed, the seminal Torres et al. ( 2010) stellar mass-radius relation, constructed using the best-characterised stars, contains only two fully convective stars: CM Dra A and CM Dra B. This is problematic not only due to a small sample size, but because both stars are "inflated" by ≈ 5%, meaning that their radius is higher than expected from theoretical models (Feiden & Chaboyer 2014b).The stars therefore may not be representative of typical low-mass stars.We will re-assess the inflation of CM Dra in Sect.7.6 (spoiler alert: they are still inflated).
CM Dra is actually a triple star system, with a bound white dwarf WD 1633+572 located 25.7" away, corresponding to a separation of at least 370 AU (Deeg et al. 2008).The white dwarf, by virtue of the known relationship between age and effective temperature, provides an age estimate of the system of 8.5 ± 3.5 Gyr (Feiden & Chaboyer 2014a).CM Dra also has sub-solar metallicity: [Fe/H]= −0.30±0.12from Terrien et al. (2012).We list observable, physical and orbital parameters in Table 1, including both our new values (Sect.4) and literature values from the previous reference (Morales et al. 2009) 1 .
CM Dra was the first target of circumbinary planet surveys (Schneider 1994;Jenkins et al. 1996;Deeg et al. 1998Deeg et al. , 2000)), with searches based on both eclipse timing variations and planetary This was interpreted by Muñoz & Lai (2015); Martin et al. (2015); Hamers et al. (2016) as evidence for a "violent", high-eccentricity evolution history of such tight binaries, which would destabilise planets.However, whilst circumbinary planets are objectively cool (Martin & Fabrycky 2021), they are not the focus of this paper.It has also been known for decades that CM Dra exhibits flares.Lacy 1977 conducted a dedicated photometric flare search on CM Dra with 18 hours of data spread over one year, covering all orbital phases with a high-speed cadence.No flares were observed in this part of the survey.However, there was also a shorter survey of differential near-infrared photometry with a different goal of improving the limb darkening and radius measurements.To quote Lacy (1977), "Due to the perversity of Nature, the only flare we have detected occurred during the infrared differential photometry when it was least expected and least desired!"Based on this single event, the flare rate was estimated to be between 0.48 and 1.2 flares per day.Subsequent studies by Metcalfe et al. (1996); Kim et al. (1997); Nelson & Caton (2007) calculated similar flare rates based on samples of typically less than ten.Stelzer et al. (2022) discovered 16 flares in one TESS sector.In our paper we will demonstrate the detection of 163 flares in 15 TESS sectors.

DATA ACQUISITION AND PREPARATION
In Fig. 1 (left) we show a cutaway of the light curve which highlights four phenomena: 1: primary eclipses; 2: secondary eclipses; 3: out of eclipse periodic variability (which we attribute to spots) and 4: flares.In Fig. 1 (right) we show the light curve phase-folded on the orbital period of 1.28 days.In the latter plot, the out of eclipse variability has been flattened but the flares and eclipses remain.
For radial velocities, we do not take new measurements, but we instead use Morales et al. (2009)'s RV measurements, which come from an improved analysis of the Metcalfe et al. (1996) spectroscopy.There are 233 measurements for each star, with a median precision of 1.2 km/s and 1.4 km/s for the primary and secondary star, respectively.With updated spectrographs and improved binary spectroscopy (e.g.Standing et al. 2022) we could improve upon these RV's, but we leave that to a future study.

COMBINED PHOTOMETRY & RV FIT
There are several steps to preparing the light curve.We provide online the light curve at every stage of the processing.
We first detrended and stitched the TESS data using the Wotan package (Hippke et al. 2019).We use the default Tukey's bi-weight filter, which Hippke et al. (2019) has shown to be the best at preserving the shape of short-timescale events (like eclipses and flares) whilst removing longer trends (e.g.spot modulation).We use  a window length of 0.228 days, which is four times the eclipse width.Two exceptions to this were for massive flare events at times 1953 and 2726 (BJD -2,457,000), for which this Wotan filter affects the flare shape.We account for this by taking the Wotan trend from the previous 1.268 days (the orbital period) and using that to detrend the light curve over this flare.The end result of this process is a flattened light curve with only eclipses and flares.
The next step was to fit the flattened photometry using the exoplanet software (Foreman-Mackey et al. 2021).This incorporates an eclipsing binary model from starry (Luger et al. 2019) with limb darkening parameters from Kipping et al. (2013).The purpose of this initial fit was not to derive stellar and orbital parameters, but rather to create a light curve model that we could subtract to remove the eclipses.Doing so produces a flat light curve with Within each TESS orbit is a parabolic trend, which we fit and remove (pink line).The remaining trends we attribute to spot modulation and possibly ellipsoidal variation.The inset shows a Lomb Scargle periodogram for the second TESS orbit in sector 51.The strongest periodicity is at  bin , followed by the longer period parabolic trend and a smaller peak at 1/2 bin .
nothing but flares.By model-subtraction of the eclipses rather than simply cutting them (as was the case in Stelzer et al. 2022), we can illuminate in-eclipse flares, as demonstrated in Fig. 2.This means we can compare the flare rate in and outside of eclipse (Sect.7.3.2) and search for occulted flares which would help us to constrain the location of the flaring region on the star (Sect.7.4).Subtraction of the eclipses also maintains ∼ 10% of the data that would otherwise be masked out.
We then used the stella neural network software (Feinstein et al. 2020a) to identify flares in this flattened, eclipse-subtracted light curve (detailed in Sect.5).We use this list of flare times and durations to cut them from the flattened light curve.We are finally left with a light curve that is flattened and contains nothing but eclipses and a few gaps where there used to be flares (both in and outside of eclipse).
On this light curve we create a joint photometry & RV fit with exoplanet.The free parameters in the model are the binary orbital period  bin " total baseline system flux  0 , the mass   and radius   of the primary star, primary-to-secondary mass ratio , radius ratio , and surface brightness ratio , epoch of first eclipse  0 , impact parameter , and eccentricity vectors ( √  sin , √  sin ).All prior distributions were chosen to be uninformative (i.e.uniform on  and ; normal on  0 ,  bin ,  0 , , , and ; log-normal on   and   ).Each star's quadratic limb darkening profile was parameterized using the uninformative prescription of Kipping (2013).
After first estimating the maximum a posteriori parameters, we use PyMC3 to derive a posterior distribution and 1 error bars, provided in Table 1.
Our radius measurement precision is 0.064% and 0.059% for the primary and secondary star, respectively.These measurements represent an order of magnitude improvement from the previous Morales et al. (2009) precisions of 0.75% and 0.63%, which were already some of the best known for M dwarfs.This is a product of the exquisite precision and long baseline of the TESS photometry.
Our precision on the stellar masses is 0.11% and 0.12% for A and B, respectively.This is an improvement with respect to the Morales et al. (2009) precisions of 0.39% and 0.47%.The improvement is not as large as it was with the radius, but that is because we improved the photometry but used the same RVs.The fact that the masses have improved despite using the same RVs is due to improved orbital parameters.For example, the precise eclipse phase from photometry more precisely constrain the eccentricity than the RVs alone.Since the conversion from RV semi-amplitude  to mass has an eccentricity dependence, an improved measurement of  improves the measurement of  A and  B .
Overall, these are the most precise M dwarf parameters ever measured and are likely amongst the most precisely-measured stars ever.There is some discrepancy between our measurements and the Morales et al. (2009) measurements.Based on the Morales et al. (2009) measurement error, our radii are both ≈ 1.3 smaller.Our masses are 6.6 and 4 smaller for A and B, respectively.Compared with Morales et al. (2009), we use new photometry but the same radial velocities.It may therefore seem surprising that we get different masses.However, as previously mentioned, changes to the photometry will still impact the masses.We have moved from ground-based photometry to an order of magnitude more space-based data.The source of the mass discrepancy might be a discrepancy in the derived impact parameter, since the radial velocity minimum mass becomes a true mass when the inclination (i.e.impact parameter) is measured.Accurate, unbiased and precise determination of the impact parameter is known to be challenging, in particular for grazing transits/eclipses, as is the case for CM Dra given the stars are practically twins (Gilbert et al. 2022a,b).This can be particularly the case when comparing derived values using different instruments, since one also has to factor in different prescriptions for limb darkening.Another possible explanation is that the stars have physically changed over the decades since the observations used in Morales et al. (2009).We know that CM Dra is active, with spots evolving on what we will derive to be a ≈ 4-year timescale (Sect.7.5).Subtle changes in the spot distribution inside and outside eclipse will change the derived physical parameters .A more complex physical model that includes both the eclipsing bodies and a modelled spot distribution, e.g. using PHOENIX (?), is being the scope of our paper.
Finally, as we will discuss in Sect.7.6, the degree to which the CM Dra stars have inflated radii with respect to models has not changed with our fit.Within each segment we create a fit to the data with a pair of sinusoids: one with a period near  orb and one with a period half that.This fitted curve is shown in dark blue.In both segments it is similar to the pink Wotan trend, indicating that our two-sinuisoid model is appropriate and well-fitted.By comparing the two sectors, it is clear that the amplitude and phase of both sinuisoids has changed significantly, indicating substantial evolution in the spots over this ∼ 900 day timespan.

FLARE DETECTION
We detect flares on the flattened, eclipse-subtracted light curve using the stella software (Feinstein et al. 2020a).This is a flare detection algorithm based on a Convoluted Neural Network (CNN), and has been applied to TESS detection in Feinstein et al. (2020bFeinstein et al. ( , 2022)).We use the default training set.stella assigns a "flare likelihood probability" to each flare, for which the default threshold is 50%.These probabilities are provided in our online table of flares.We detect 175 flares above this default threshold.However, within this paper we instead use a more conservative threshold of 75%, from which we detect 163 flares.There are several reasons for a stricter threshold.First, by eye many of the lower probability flares look unconvincing above typical noise.Second, when calculating the flare rate in Sect.7 we are calculating the flare frequency distribution, with accounts for incomplete flare detection at low energies, i.e. the flares likely to have low probability from stella.Finally, when calculating the distribution of flares (Sect.7.3) we prefer a conservative approach where the flare statistics may be weaker but we are more confident that they are actually flares.
We follow Shibayama et al. ( 2013); Günther et al. (2020) to calculate the energy of each flare.The flare is modelled as black body radiation with  = 9000 K. Given we do not know which star a given flare is occurring on, and they are close to identical, we take the stellar radius and temperature to be the average between stars A and B. In general, the luminosity of a blackbody with in a given observing bandpass (i.e.TESS) is calculated by where  is the surface area of the blackbody,   is the TESS response function, combining both the detector filter transmission and the quantum efficiency, and () is the Planck function as a function of the temperature of the blackbody temperature .This equation, evaluated for the flare and star, is where the above equation demonstrates that the flare luminosity and hence area will be time-dependent, compared with the constant, quiescent stellar luminosity.
In the normalised light curve the change in flux is related to the change in luminosity by The factor of 2 differs from Shibayama et al. ( 2013); Günther et al. (2020) to account for the fact that CM Dra contains two stars of nearly equal brightness.For example, a 10% increase in the light curve flux corresponds to a 20% increase in the brightness of the individual flaring star.We therefore solve for the area of the flare region: The bolometric flare luminosity is given by where  SB is the Stefan-Boltzmann constant.By integrating the luminosity over the flare time we obtain the total flare energy: Our 163 detected flares span energies between 9.1 × 10 30 and 2.4 × 10 33 erg, as discussed in Sect.7. The largest flare area is ≈ 7 × 10 12  2 , which corresponds to a diameter that is ≈ 10% that of the star.This is a much larger spot-to-star ratio than we see on the Sun, but smaller than the highly active M-dwarfs studied by Ilin et al. (2021).

SPOT DETECTION
CM Draconis is a so-called "BY Draconis" variable, which is a main sequence variable (typically a K or M dwarf) exhibiting variations on the order of roughly a magnitude.The variation is caused by the presence of star spots.Since spots evolve and may exist at different latitudes (and hence different rotation rates), this class of variable is only "semi-regular".This is in contrast with say cepheid and ellipsoidal variables, which are "fully regular".
Since CM Dra is such a tight binary, we expect the two stars to be 1:1 spin-orbit synchronised due to tides.We also expect the effect of spots to be visible on CM Dra given M dwarfs are highly active and it is a bright target.Overall, we expect an out of eclipse variability with a period of 1.268 days.
In Fig 3 we plot sectors 25 and 51 of TESS.They are separated by approximately two years.Each sector contains two of TESS's 13.7-day orbits around the Earth.There is a gap in the middle of each sector for data download 2 .For clarity the eclipses have been removed.Several things are visibly apparent.First, there is a modulation of the light curve with a periodicity a bit longer than 1 day, as expected.However, it is more complicated than a simple sinusoid.Within each TESS orbit there is a longer term parabolic modulation of the light curve (fitted pink line).This changes from TESS orbit to orbit.Given this occurs on a TESS orbital timescale, we attribute it to systematics in the data and not something physical with CM Dra.
Even if we account for these TESS systematics, it is visible that the shape of the photometric modulation is not strictly a single sinusoid.This is emphasised in the Lomb-Scargle periodogram (Fig. 3 inset) for the second TESS orbit in sector 51, where we see the three most significant powers are i)  bin , ii) ≈ 14 days due to the parabolic trend of TESS systematics, and iii) 1/2 bin .Finally, the amplitude of the modulation is visibly larger in sector 51 than it is in 25, suggesting that the spot coverage has become more pronounced and/or uneven over time.
To quantify the out-of-eclipse light curve modulation, and its time evolution, we invoke the following procedure: (i) Cut (not subtract) all eclipses from the light curve.(ii) Split the data up into separate TESS orbits and fit and subtract a parabola (like Fig. 3).In some TESS orbits there is a sharp change to the flux near the start or end of the orbit, beyond the simple parabolic trend, in which case we manually cut these edge data.
(iii) Further split the data up into 4 × 1.268 day segments, where 1.268 days is CM Dra's orbital period.The choice of spitting the data up into 4 ×  bin was an ad hoc compromise between having a sufficiently long timeseries to make a good sinusoid fit and not extending the timeseries too long such that we would blur what we will see is fairly rapid spot evolution.
(iv) Within each segment, we fit two sinusoids using scipy.optimize.curve_fit,using 1.268 and 0.634 days as the initial period guesses (like Fig. 4).
In Fig. 5 we show the period (top), phase (middle) and amplitude 2 This gap is larger in sector 51 because of a change to the TESS mission to take more short cadence data, which takes longer to download to Earth.(bottom) of the two sinusoids fitted to each 4 ×  bin segment.The periods are consistently close to either  bin (1.268 days, blue circles) or 1/2 bin (0.634 days, pink squares).The phases and amplitudes of both the  bin and 1/2 bin signals are seen to change, both between cycles and within cycles.We analyse what we believe to be the source of these two signals in Sect.7.1.

What Causes The Light Curve Modulation?
We attribute the 1.268 day signal (Fig. 5) to star spots, as expected from two tidally locked M dwarfs with a 1.268-day orbital period.We do not know which star is producing the spot signal, but it is likely that both stars have a high spot coverage and the 1.268 signal is an amalgamation of the two given that the two stars are so similar in size.The spot modulation changes both in phase and amplitude, which we analyse in Sect 7.5 in the context of a possible activity cycle.The trickier question is what causes the 1/2 bin -day frequency?
For a tight binary we would typically ascribe a frequency at 1/2 bin to ellipsoidal variation.This is the tidal deformation of the star such that at different orbital phases we see different surface areas, and hence different brightnesses.The ellipsoidal signal should have three attributes: 1) a period at exactly 1/2 bin , 2) a phase such that the signal is maximised at orbital phases 0.0 and 0.5 (i.e. it is minimised during eclipses at phases 0.25 and 0.75) and 3) a constant amplitude: where  ellip = 0.15 is calculated using the stellar gravity darkening coefficient  and the limb darkening coefficient  within the TESS bandpass.Equation 9is taken from Mazeh & Faigler (2010).We added a factor of 2 to account for the fact that we have two identical stars inducing ellipsoidal variation in each other, as opposed to the star-exoplanet binaries considered in Mazeh & Faigler (2010).We calculate  ellip ≈ 0.0007.With this in mind, we consider the following possibilities for this 1/2 bin -day signal: (i) Ellipsoidal variation: Within Cycle 2 we see the 0.634-day signal looks roughly as expected, with respect to its period, phase and constant amplitude.In Cycles 4 and 5, however, things look peculiar.The period remains clearly 0.634 days, yet the average amplitude has roughly doubled.Furthermore, the phase has changed to an average of roughly 0.6.This may seem close to the expected 0.5, but as demonstrated in Fig. 4 the minimum is noticeably offset from the eclipses.We conclude that the power of the 1/2 bin frequency is not (at least predominantly) coming from ellipsoidal variation.
(ii) Reflection or Doppler beaming: These other binary-specific effects occur with a period of  bin , not 1/2 bin , so they do not work.
(iii) Dilution from a neighbouring star: One might suggest that a change in amplitude of the variation could be caused by varying dilution from a neighbouring star which we have not accounted for.However, we can rule this out based on CM Dra being fairly bright (Jmag = 8.5) and any such time-dependent dilution would also affect the eclipse depths, which is not seen.Dilution could also not explain a change in phase.
(iv) Each star has a different rotation rate: If both stars are spotted but have different rotation rates then we would expect multiple frequencies.However, we are unaware of any tidal mechanisms in a tight twin binary that would produce a 1 : 1 spin-orbit ratio in one star and a 2 : 1 spin-orbit ratio in the other.
(v) Differential rotation: Spots can manifest as multiple frequencies if they occur at different latitudes in a star with differential rotation.This is seen on the Sun.However, to have one latitude of CM Dra rotate at  bin and another at precisely 1/2 bin seems implausible.
(vi) Spot clusters separated by 180 • in longitude on a single star: If spots are not uniformly distributed in longitude but rather in two distinct clumps on opposite sides of the star then this would manifest as a 1/2 bin signal.
(vii) Spot clusters at sub-stellar points: If both stars have an overabundance of spots near the substellar point then we would predict two things.First, this would produce a 1/2 bin periodicity since the observer would see each cluster once during an orbital period, separated by half an orbital period.Second, we would expect the 1/2 bin signal to have a maximum flux near orbital phases 0.0 and 0.5 (i.e. a minimum flux near eclipse at phases 0.25 and 0.75).This is close to what we see in Fig. 5.
The two most plausible explanations are (vi) and (vii), which we

A B
Spot clusters separated by 180° on a single star (vi)

A B
Source of the 1/2 Pbin Spot Signal AU Mic-like effect Spot clusters at sub-stellar points (vii) Figure 6.The two most plausible explanations for a spot modulation signal at 1/2 bin .These are explanations (vi) and (vii) from our list of seven in Sect.7.1.The primary (A, blue) and secondary (B, pink) stars are assumed to have spots across many longitudes, which will produce the most prominent signal at  bin .Additionally, there is an over-abundance of spots clustered in longitude, denoted in this diagram by a black strip.Top: one star (we do not know which) has concentrations of spots separated in longitude by 180 • .This effect was seen prominently in the single star AU Mic, shown in the inset.Bottom: both stars have spot clusters near the sub-stellar point.The sub-stellar point is constant due to tidal locking.
illustrate in Fig. 6.Explanation (vi), with spot clusters separated by 180 • on a single star, is plausible because this effect has been seen on single stars such as AU Mic, which is another BY Draconis variable (Fig. 6 inset).There is clearly a primary periodicity at 4.8 days, with a secondary periodicity at half that.Martioli et al. (2021); Szabó et al. (2021) also attributed this to 180 • -separated spots.Skarka et al. (2022) noted this effect in 31 magnetic rotator (ROTM) variables out of a few hundred A-F spectral type stars.It was also noted that this effect may be confused with ellipsoidal variables of non-eclipsing binaries.With TESS we also see this 1/2 periodicity effect in BY Draconis itself, which we will study in a future paper.Explanation (vii), which is only for binary stars, is also plausible.It was proposed by Simon et al. (1980); van den Oord (1988); Gunn et al. (1997) that there may be magnetic reconnection field lines that connect two stars in a tight, tidally locked binary, and this may increase spot coverage near the sub-stellar points, and hence a reduced flux near eclipse.We note that this effect would be easily confused with ellipsoidal variation.The main difference is that whilst ellipsoidal variation must have a minimum exactly during eclipse, it is reasonable that the spot cluster may not be precisely on the sub-stellar point, and hence the minimum can be just near eclipse.
The fact that CM Dra is a binary does indeed confuse things.In the future if we see the phase of the 1/2 bin signal change significantly then this would favour explanation (vi), since this explanation has no restriction on the spot phase whereas (vii) dictates the phase with respect to the eclipses.We ultimately leave it to a future study to construct a physical evolving spot model that is a good fit to the CM Dra data.

Flare Rate
Our time-series spans 1198 days.By cutting out gaps between and within sectors, there are 328 days worth of observations.With a total of 163 flares, we deduce a simple flare rate of 0.5 flares per day.We recall that Lacy (1977) calculated a flare rate between 0.48 and 1.2 flares per day based on a single flare, which is compatible with our result from 163 flares.
One important subtlety is that this is for the CM Dra binary as a whole.If we are to compare this with flare rates of single M dwarfs we must account for the fact that it is a binary and hence two potential sources of flares.We will argue in Sect.7.3.2that the TESS data demonstrate flares occur on both stars.This is not surprising given that the stars are practically twins.We therefore approximate the flare rate on each individual M dwarf as 0.25 flares per day, i.e. half the total rate.
A more sophisticated way of quantifying the flare rate is to calculate the flare frequency distribution (originally proposed by Lacy et al. 1976): where  is the flare frequency and  is the flare energy.For our fit  = 40.40 and  = −1.28.This is plotted in Fig. 7.This accounts for flares being produced with a distribution of energies.It also accounts for detection limitations, i.e. our flare sample will not be complete at low flare energies.To assist interpretation of Fig. 7 the horizontal lines demarcate the energy of flares that occur on a daily, weekly and monthly basis.The pink line is a power law fit.It is fitted to the higher energy flares (in this case > 10 32 erg) because they are the easiest to detect, and hence we believe the detections are complete in this parameter space.At lower energies the curvature of the flare frequency distribution indicates that the flare sample is not complete because of limits in detection sensitivity when it comes to low energy flares.The flare frequency distribution allows for comparisons with other stars.MacGregor et al. (2021) studied our nearest stellar neighbour, Proxima Centauri, and derived a flare frequency distribution of log 10 () = 27.2 − 0.87 log 10 ().Compared with CM Dra, Proxima Cen has a smaller  but a larger , i.e. a less steep negative slope.This means that CM Dra has more frequent low-energy flares, but Proxima Cen has more frequent high-energy flares.

Flares as a Function of Phase
We are interested in the flare distribution with respect to two phases: the orbital phase and the spot modulation phase.Even though the orbital and main spot periods are the same, the spot phase changes over the TESS timeseries.We therefore analyse the two phases separately.That is to say, flares with an energy of 10 31.7 erg and above occur on a daily basis, flares over 10 32.3 erg occur on a weekly basis, flares over 10 32.7 occur on a monthly basis, and so on.

Spot Phase
Since flares and spots are expected to be associated, one might expect that more flares will occur when the spot modulation flux is negative.This is because a negative flux implies a higher-than-average concentration of spots facing the observer, and hence associated flares would also face the observer.However, we see in Fig. 8 (top) that in fact the opposite is true: 104 flares (64%) coincide with a positive spot modulation flux.The spot modulation flux is taken as the wotan trend fitted in Sect. 3 with the parabolic TESS-orbit trends removed according to Fig. 3.There is a slight asymmetry in the spot modulation flux overall -52% of all data points are positive -but this is not enough to account for the asymmetry in Fig. 8 (top).
We quantify the significance of our apparent correlation between flare occurrence and positive spot flux using a simple bootstrap test.We randomly draw 163 flares across our time series and test if 104 or more flares correspond to a positive flux.We repeat this process 10 4 times.From this derive a p-value of 0.0083.This implies good evidence of a positive correlation between the flare rate and positive spot modulation.
As a second test, we conduct a one-sample Student t-test, comparing the distribution of the spot flux of the 163 flares to a normal distribution with a mean of zero.We return  stat = 1.942 and a p-value of 0.054.We interpret this as moderate evidence that the flares are not normally distributed around zero spot modulation flux.We also note a caveat that an assumption of the Student t-test is that events are independent, whereas flares may exhibit "sympathetic flaring", where the production of one flare may increase the likelihood of subsequent flares shortly after (Moon et al. 2002;Feinstein et al. 2022).We do not analyse sympathetic flaring in CM Dra.
Overall, our data provide preliminary evidence that on CM Dra flare rates are inversely proportional to spot occurrence, which is contrary to the expectation.However, we recommend repeating this analysis with future TESS data before drawing strong conclusions.Fortunately, TESS will observe CM Dra for 7 sectors in Cycle 6, starting March 2024.As previously stated, one might expect more flares when there is a negative spot modulation (i.e.negative flux means more spots are visible and flares are expected to be more likely near spots).If there truly are more flares during positive spot modulation, then this would have interesting implications for the models of stellar activity and surface inhomogeneities.Some recent literature studies.Dal & Evren (2011); Hawley et al. (2014); Doyle et al. (2018) found found no correlation between flare rate and rotational phase.Roettenbacher & Vida (2018) found no correlation for the stronger flares and only weak correlation for the weaker flares.One possible explanation is that the flares are indeed correlated with dark starspots, but they are also correlated with bright spots known as faculae/plages.Starspots are dark spots on the stellar photosphere, resulting in a flux deficit.Faculae are conversely bright spots on the photosphere.Plages are corresponding bright spots observed higher up in the stellar chromosphere.Our Sun, for example, actually gets brighter when there are more sunspots.Whilst this might seem unintuitive, this is because the presence of spots is accompanied by the presence of faculae/plages.For the Sun, the excess brightness of plages/faculae outweighs the flux defecit of the spots, resulting in a net increased brightness of 50% during peak stellar activity in the Solar cycle (Chapman et al. 1997).It would be of interest to test this trend in other M-dwarf stars, but we leave that for future studies.

Orbital Phase
In Fig. 8 (middle) we plot a histogram of the flare count as a function of orbital phase.The choice of bins is such that the bins centered at phase 0.25 and 0.75 (highlighted in aqua) correspond to an eclipse duration (1.37 hours).The mean flare count across all non-eclipse bins is 7.10.By eye, there is no obvious phase preferred by flares.In Fig. 8 (bottom) we calculate the cumulative distribution function (CDF) across all non-eclipse histogram bins.We compare this with a Poisson distribution for five different  parameters.We see that the non-eclipse flares can be modelled roughly by a Poisson distribution with  = 7.5.We therefore conclude that the flare distribution as a function of orbital phase is flat aside for Poisson noise, i.e. there is no favoured orbital phase.A flat orbital phase distribution of flares is in contrast to the spot distribution, which seemingly has some clustering (Sect.6).
We treat the flare counts within the primary and secondary eclipse separately because they should not be representative of the rest of the orbit.Within each eclipse one of the stars is almost completely obscured.Depending on the distribution of flares (longitude, latitude and across both stars) we should expect some reduction during eclipse.We do not see a reduced flare rate during eclipse.
Seven flares are seen during primary eclipses and nine are seen during secondary eclipses.These counts are higher than the mean of 7.10 across non-eclipse phases.If we take the 22 bins and fold them across a phase of 0.5, there are 21 flares in eclipse, whereas there out of eclipse bins have 14.2 ± 3.5 flares, where  = 3.5 is the standard deviation of the bin counts.We discuss several possible flare distributions and their compatibility with our surprisingly high in-eclipse flare rate: (ii) Flares concentrated at sub-stellar points.Simon et al. (1980);van den Oord (1988); Gunn et al. (1997) proposed that there may be magnetic reconnection field lines that connect two stars in a tight, tidally locked binary, and this may increase flare rates near the sub-stellar points.This might be connected to the possible spot clustering near the sub-stellar points, illustrated in Fig. 6.One would therefore expect more flares at phases 0.25 and 0.75, hence explaining our high in-eclipse flare count even if some of them are obscured.However, a concentration of flares near sub-stellar points should not produce the flat distribution of flares seen outside eclipse, so this scenario seems to be unlikely.

Flares During Primary Eclipse
Flares During Secondary Eclipse (iii) Only one star flares.In this case one eclipse would have no flare reduction and the other would have a 50% reduction (marginalised over the entire eclipse bin).Since there is a similar number of flares during each eclipse (seven and nine) and both are above the mean out of eclipse flare count, this scenario also seems unlikely.
(iv) Confusion with spot crossings.If a spot is occulted during an eclipse then there is an upwards bump in the light curve since a relatively dark region is being blocked.The morphology would be different to a typical flare, which has a sharp rise and exponential decay, but we speculate that stella could possibly confuse the events.This would imply that some of our 16 in-eclipse "flares" are in fact spot crossings, and hence the true flare rate may actually be reduced.We do not favour this scenario because most of the events in Fig. 9 truly look like flares.
(v) Polar flares.Since CM Dra has such a small impact parameter ( = 0.11), if the flares are predominately at high latitudes, and the two stars are spin-orbit aligned3 , then there will be less blocking of the flares during eclipse.Therefore, in this scenario the in-eclipse flare count should be the same as the out of eclipse flare count.This explanation is the most compatible with our data, but we do not consider this solution conclusive.A polar distribution of flares would affect the timing of flares observed during eclipse, i.e. we would see more flares at the start and end of the eclipse compared with in the middle.However, conducting a significant analysis on our small sample of 21 flares observed during eclipse is not possible, particularly when you consider the fact that both stars could flare and hence a flare in the middle of an eclipse could be on the nearside of the foreground, occulting star.
Evidence for the latitudinal distribution of flares on M dwarfs is currently minimal.Ilin et al. (2021) discovered six flares on four M dwarfs with latitudes above 50 • .This discovery came using rapidly-rotating stars ( rot ≈ hours) with energetic flares on a longer timescale.The rotation brings the flare in and out of view, modulating the light curve and allowing the latitude to be isolated.Our discovered flare distribution is consistent with polar flares, and hence the results of Ilin et al. (2021).However, our evidence is only marginal and is a less direct method of constraining flare latitudes.
Huang & Ormel (2022) studied 12 M dwarf eclipsing binaries and found that only some showed a relationship between flare rate and orbital dependence.There were also a couple of targets for which these events.Spot crossings could also possibly be confused with in-eclipse flares.the flare rate seemed to drop during eclipse.These results were not interpreted within the context of flare latitudes.

Occulted Flares
In Sect.7.3.2we counted 21 flares during primary and secondary eclipses.Depending on the timing of flare, it is possible that the flare will be occulted by the foreground star, imprinting a distinct signature on the light curve.We will investigate this effect in a dedicated future paper (Armitage et al., under review), but for now we conduct a brief investigation into CM Dra.From Sect. 5, the median surface area of the flares (at max flux) is 6.9 × 10 13 m 2 .Naively assuming circular flare regions, this corresponds to a diameter of 4.7 × 10 6 m, which is 1.3% of CM Dra A's diameter.During eclipse, the relative motion of the two stars is ≈ 150 km/s, meaning that the entire flare region will be passed over a timespan of ≈ 30 seconds.With 120-second cadence data the flare ingress/egress of the occultation would therefore be effectively instantaneous.
In Armitage et al. (under review) we detail the phenomenon of flare occultations, including the possible light curve morphologies.In the simplest cases though, if a flare is on the decay phase and is then occulted you will see a sharp and noticeable drop in flux, assuming the occultation occurs when the flare is still sufficiently bright.The main alternative is that the flare starts when it is being occulted by the foreground star, and subsequently the flare is revealed as the star passes by.In such a scenario the rapid egress of the occultation may be indistinguishable from the rapid rise of a standard flare, and hence we may not be able to deduce that the flare was occulted.
In Fig. 9 we plot all 21 flares discovered during eclipse.Most of the flares occuring within the eclipse do not show signs of being occulted.For this to be the case they must either occur on the nearside of the foreground star, or on nearside of the background star but on an uncovered region of its surface.Flare #121, on a primary eclipse, may be occulted because there appears to be a sharp drop in flux, as opposed to the typical gradual drop in flux (most clearly evidenced in flares #1 and #79).We highlight this flare in green and then plot the 20 second data separately in Fig. 10.For flare #121 we also plot the 20 second data in (highlighted in green), where in this plot the .Bottom: semi-amplitude of the out of eclipse spot modulation, calculated as half of the maximum peak to peak difference of the wotan trend within a given TESS orbit.This trend encompasses both the 1.268 and 0.634-day frequencies, but parabolic intra-orbit trends (Fig. 3) have been removed.We see a significant increase in the amplitude of spot modulation over time, which is accompanied by a modest increase in the flare rate.The mean flare energy, however, does not vary.Based on the spot modulation we estimate a 4.1-year activity cycle.
eclipse model is shown but not subtracted.There appears to be a typical sharp rise in the flare and then a decay over ≈ 200 seconds.This is a little longer than the ≈ 30 second calculation from before, but we emphasise that that was a rough estimate.If we compare flare #121 to say flare #1 we see that both stars have a similar peak amplitude, yet flare #1 takes a much longer ≈ 3500 seconds to decay.
One challenge with CM Dra though is that since effectively all latitudes of each star are blocked at some point of the eclipse, the observation of an occulted flare does not necessarily identify its latitude.However, since the potential occultation of flare #122 occurs early in the eclipse, the flare must have been near the star's equator (assuming spin-orbit alignment).This is in contrast with our deduction of polar flares in Sect.7.3.2,but ultimately neither piece of evidence is conclusive.
An occulted flare from a transiting exoplanet more clearly identifies the flare latitude because planets only cover a narrow transit chord.Although, one does still have the question of the spinorbit alignment of the system.Flares may also have complex shapes irrespective of any occultations (e.g. if there are multiple small flares simultaneously).We leave a thorough analysis of occulted flares to Armitage et al. (in prep).

Activity Cycle
On the Sun the occurrence of flares and spots varies on an 11-year cycle.We test this in CM Dra in Fig. 11 by plotting the flare rate, mean flare energy and spot modulation semi-amplitude over time.The data is split up to individual 13.7-day TESS orbits (two per sector).
The spot modulation has the strongest variation.This is essentially the same signal as shown in Fig. 5 (bottom) except the 1.268 and 0.634-day signals are combined.There is a slight increase within cycle 2.The jump from Cycle 2 to Cycle 4 corresponds to a large increase in flux, followed by a decrease throughout Cycle 4 and into Cycle 5. We can fit a 4.1-year activity cycle to these data, but we stress that a longer time-series is needed to properly constrain this.There is also ambiguity arising from the fact that we cannot assign spots to individual stars, and it is not guaranteed that CM Dra A and CM Dra B would have synchronised and identical activity cycles.The flare rate is seen to increase slightly over time, but the mean flare energy does not seen to change.

Is CM Draconis Still Inflated?
Yes, with our more precise parameters, the CM Dra radii are still inflated relative to stellar models4 .This is demonstrated in Fig. 12, where we plot the mass and radius of CM Dra A and B from both our work and Morales et al. (2009), along with some other precise M dwarf measurements from the literature (Parsons et al. 2018;Duck et al. 2023).The comparative model is from the Mesa Isochrone and Stellar Track (MIST, Dotter 2016).Several well-characterised M dwarfs have larger (inflated) radii than expected at given masses.
This MIST model is with a metallicity [Fe/H] = 0.One may argue that deviations from the model are a result of a different metallicity.However, in the case of CM Dra this argument exacerbates the inflation; Terrien et al. (2012) derived [Fe/H] = −0.30± 0.12 for CM Dra, and at a sub-solar metallicity we would expect an even smaller radius.Feiden & Chaboyer (2014b) suggest that the magnetic activity of CM Dra may be the cause of its inflation.Whilst we do not explicitly test this effect in this paper, this is further motivation to study CM Dra's stellar activity.

Probing Activity Through Spectroscopy
We have solely used photometry to observe flares, but they can also be observed spectroscopically by detecting emission lines.Metcalfe et al. (1996) observed emission lines for Mg I and Fe II for CM Dra.Based on the radial velocity, the flares could be attributed to the primary star A, which is advantageous compared with photometry for which it is typically ambiguous which star caused any given flare, particularly for near twins like CM Dra.One may also use spectroscopy to calculate spot filling factors (Somers & Pinsonneault 2015;Cao & Pinsonneault 2022), where the spectrum is fitted with an additional spot component.Combined with a measurement of the projected rotational velocity  sin , this technique may reveal the spot latitude.Whilst we do not assess CM Dra's activity spectroscopically, we encourage future studies, in particular if they are simultaneous with future photometry from TESS or the PLAnetary Transits and Oscillations of stars (PLATO).either an anti-correlation between flares and dark sunspots, or a positive correlation between flares and bright plages/faculae.
We also discover changes to the spot activity over a 4.1-year, which may be indicative of an activity cycle similar to the Sun.The flare rate is also seen to change slightly over this activity cycle, but no variation is seen in the average flare energy.Additional data from TESS, including 7 sectors in the upcoming Cycle 6, will improve our flare and spot statistics and provide a better probe of any activity cycles.

Figure 1 .
Figure 1.Left: segment of PDCSAP photometry, illustrating the four major components of the CM Dra light curve: primary eclipses, secondary eclipses, flares and out-of-eclipse modulation, the last of which is predominantly due to spots on the rotating stars.Right: CM Dra light curve that has been detrended (i.e.flattened) and phase-folded on the orbital period of 1.268 days.The primary and secondary depths are similar, almost corresponding to a 50% reduction in flux because the orbital alignment is almost perfectly edge-on, and the two stars are almost identical.Flares are seen to occur at all orbital phases (analysed in depth in Sect.7.3.2).

Figure 2 .
Figure 2. Eclipse model fitted with exoplanet (pink line) to the flattened TESS data (blue data).Since there are well over 100 eclipses, the occasional flare that coincides with an eclipse does not affect the model.Therefore, when the eclipse model is subtracted the flare stands out (inset).

Figure 3 .
Figure 3. TESS PDCSAP photometry of CM Dra where the only processing is the eclipses have been cut out.On the left is sector 25 and on the right is sector 51, spaced roughly two years apart.Each sector contains two 13.7-day orbits of the TESS spacecraft.The central gap is when the data is downloaded to Earth.Within each TESS orbit is a parabolic trend, which we fit and remove (pink line).The remaining trends we attribute to spot modulation and possibly ellipsoidal variation.The inset shows a Lomb Scargle periodogram for the second TESS orbit in sector 51.The strongest periodicity is at  bin , followed by the longer period parabolic trend and a smaller peak at 1/2 bin .

Figure 4 .
Figure 4. Two example lightcurve segments: one from TESS sector 19 (left) and one from TESS sector 52 (right), over 900 days later.Each segment has a length equal to four orbital periods of CM Dra (4 × 1.268 = 5.072 days).Both types of eclipses have been cut out, with solid and dashed vertical lines used to denote when the primary and secondary eclipses occur, respectively.The pink curves indicate the lightcurve trend calculated using a Tukey's biweight filter in Wotan.Within each segment we create a fit to the data with a pair of sinusoids: one with a period near  orb and one with a period half that.This fitted curve is shown in dark blue.In both segments it is similar to the pink Wotan trend, indicating that our two-sinuisoid model is appropriate and well-fitted.By comparing the two sectors, it is clear that the amplitude and phase of both sinuisoids has changed significantly, indicating substantial evolution in the spots over this ∼ 900 day timespan.

Figure 5 .
Figure 5. Activity variation over time.Each data point comes from the two sine curve fit to a 5.072-day segment (= 4 ×  bin ) of the light curve, as demonstrated in Fig. 4. Errorbars are not visible for most points.Top: period of the fitted sinusoids.The strongest frequency is at  bin , shown as blue circles.There is a weaker but significant frequency at 1/2 bin , shown as pink squares.Middle: orbital phase of the binary where the sinusoid is at a maximum.At each time step the 1/2 bin frequency has two max phases within one CM Dra orbital period.Bottom: semi-amplitude of the fitted sinusoids.For comparison, we show the expected ellipsoidal variation amplitude (Eq.9), which matches in TESS cycle 2 but not in cycles 4 and 5.

Figure 7 .
Figure7.Flare frequency distribution, with both horizontal and vertical axes in log scales.The vertical axis is a cumulative rate of flares up to a certain energy.That is to say, flares with an energy of 10 31.7 erg and above occur on a daily basis, flares over 10 32.3 erg occur on a weekly basis, flares over 10 32.7 occur on a monthly basis, and so on.

FlaresFigure 8 .
Figure 8. Top: flares as a function of the spot modulation flux at the time of flare.The spot modulation flux is taken from the Wotan fit after removing parabolic TESS orbit trends.One might expect more flares during negative flux, because that means more spots are facing the observer.Instead, the opposite is seen.Middle: flares as a function of orbital phase.One might expect less flares during the primary and secondary eclipses, but again, the opposite is seen.Bottom: cumulative distribution function (CDF) of flares within each orbital phase histogram bin from the middle plot.Flares within the eclipses are excluded.For comparison, we show different poisson distributions, parameterised by .

( i )
Uniform flares across both stars.If flares are spread across all latitudes and longitudes then we should expect a 50% reduction in the flare rate at eclipse totality, since the background star is completely blocked and the only flares can come from the nearside of the foreground star.If we marginalise over the entire V-shaped eclipse (Fig. 1, right), accounting for ingress and egress, then there should be a 25% reduction in the flare rate within the entire eclipse bin.Therefore, we should expect 10.65 flares across both eclipses.Our measurement of 21 flares is incompatible with this by   3.

Figure 9 .
Figure 9. Flares where the peak of the flare (as detected by stella) occurs within a primary (top) or secondary (bottom) eclipse.The vertical lines demarcate the eclipse duration.The numbers in the top left corners indicate the flare number, out of the 125 flares discovered.

Figure 10 .
Figure 10.A candidate occulted flare as seen with 20 second cadence TESS data.This is flare #121 out of 163 detected.It occurs at the start of a primary eclipse and has a very short duration (≈ 200 seconds), potentially indicative of it being occulted by the foreground secondary star.This flare is highlighted in green in Fig. 9.

ActivityFigure 11 .
Figure 11.Activity cycle of CM Dra.All of the data points are calculated within a given 13.7-day TESS orbit.Top: average flare rate.Middle: average flare energy (erg).Bottom: semi-amplitude of the out of eclipse spot modulation, calculated as half of the maximum peak to peak difference of the wotan trend within a given TESS orbit.This trend encompasses both the 1.268 and 0.634-day frequencies, but parabolic intra-orbit trends (Fig.3) have been removed.We see a significant increase in the amplitude of spot modulation over time, which is accompanied by a modest increase in the flare rate.The mean flare energy, however, does not vary.Based on the spot modulation we estimate a 4.1-year activity cycle.