Galactic Kinematics and Observed Flare Rates of a Volume-Complete Sample of Mid-to-Late M-dwarfs: Constraints on the History of the Stellar Radiation Environment of Planets Orbiting Low-mass Stars

We present a study of the relationship between galactic kinematics, flare rates, chromospheric activity, and rotation periods for a volume-complete, nearly all-sky sample of 219 single stars within 15 parsecs and with masses between 0.1$-$0.3 M$_{odot}l$ observed during the primary mission of TESS. We find all stars are consistent with a common value of $\alpha$=1.984 $\pm$ 0.019 for the exponent of the flare frequency distribution. Using our measured stellar radial velocities and Gaia astrometry, we determine galactic UVW space motions. We find 64% of stars are members of the Galactic thin disk, 5% belong to the thick disk, and for the remaining 31%, we cannot confidently assign membership to either component. If we assume star formation has been constant in the thin disk for the past 8 Gyr, then based on the fraction that we observe to be active, we estimate the average age at which these stars transition from the saturated to the unsaturated flaring regime to be 2.4 $\pm$ 0.3 Gyr. This is consistent with the ages that we assign from galactic kinematics: We find that stars with Prot $<$ 10 days have an age of 2.0 $\pm$ 1.2 Gyr, stars with 10 $<$ Prot $\leq$ 90 days have an age of 5.6 $\pm$ 2.7 Gyr, and stars with Prot $>$ 90 days have an age of 12.9 $\pm$ 3.5 Gyr. We find that the average age of stars with Prot $<$ 10 days increases with decreasing stellar mass from 0.6 $\pm$ 0.3 Gyr (0.2 - 0.3 M$_{odot}$) to 2.3 $\pm$ 1.3 Gyr (0.1--0.2 M$_{odot}l$).


INTRODUCTION
Main-sequence stars with convective envelopes display magnetically induced phenomena, such as star spots, coronal x-ray emission, and flares. Stellar flares are believed to be the result of the reconnection of magnetic field lines (Martens & Kuin 1989). Flares are observed at all wavelengths extending from x-ray to radio. With dedicated space-based photometric missions like Kepler/K2, and the Transiting Exoplanet Survey Satellite (TESS ), white-light flares, which are flares that that emit broadband emission extending from the nearultraviolet to optical wavelengths, have been studied extensively (e.g. Hawley et al. 2014;Davenport et al. 2014;Davenport 2016;Ilin et al. 2019;Davenport et al. 2019;Günther et al. 2020;Feinstein et al. 2020;Medina et al. 2020).
Previous studies have found that both the flare rate and the total luminosity emitted by all flares summed Corresponding author: Amber Medina amber.medina@austin.utexas.edu together for a given star over some length of time decreases with increasing rotation period (Davenport 2016;Ilin et al. 2019;Davenport et al. 2019;Howard et al. 2019;Ilin et al. 2021;Medina et al. 2020), but this relationship is not simple. Stars show a saturated level of magnetic activity that is independent of the rotation period up to a critical value; beyond this value, stars transition to the unsaturated regime where coronal xray emission (Wright et al. 2011(Wright et al. , 2018, ultraviolet emission (France et al. 2018), and chromospheric Hα emission (Douglas et al. 2014;Newton et al. 2017) decreases rapidly with increasing rotation period or Rossby number, Ro. Rossby number is a mass-independent proxy for rotation that is defined as the rotation period normalized by the convective turnover time of the star. In Medina et al. (2020), hereafter denoted M20, we studied a sample of 125 fully convective M dwarfs with masses between 0.1 and 0.3 M that reside within 15 parsecs and were observed in the southern ecliptic hemisphere during the first year of the TESS mission and found that this relationship extends to optical flare rates as well. Reiners et al. (2022) found that the underlying physical mechanism responsible for observed saturation of magnetic activity can be attributed to saturation of the magnetic dynamo itself. However, the age at which the stars transition from the saturated regime to the unsaturated regime, which could have a substantial impact on the planets orbiting these stars is not firmly established. Skumanich (1972) provided a model to explain observations that indicated stellar rotation periods increase as a function of age. This means that in theory the age of a star can be inferred if the rotation period is known; this method for age determination is known as Gyrochronology (Barnes 2003). By observing the rotational evolution of stars in clusters with different ages, previous studies (e.g. Barnes 2003;Mamajek & Hillenbrand 2008;Meibom et al. 2011;Agüeros et al. 2018) determined that stars with the same mass converge to a single rotation period value by an age that is dependent on stellar mass. For example, the rotation periods of solar-type stars converge to a single value by the age of 650 Myr, the age of the Hyades cluster, in contrast to the rotation periods of fully convective M-dwarfs which have not converged to a single value by this age. The spread in rotation periods makes age estimation using the technique of gyrochronolgy unreliable for the fully convective M dwarfs that are the subject of our study. However, galactic kinematics can provide a means to estimate the age of a population of M dwarfs based on their motion through and location in the Galaxy. The Galaxy is believed to be made up of three populations: the thin disk, thick disk, and halo, with each forming at different times throughout the evolution of the Galaxy. The stars in each population also have velocity dispersions that differ from each other due to their birth environments. Newton et al. (2016) used pre-Gaia parallaxes and proper motion information and found that most fully convective M dwarfs with rotation periods that exceed 70 days have ages in excess of 5 Gyr, while stars with rotation periods less than 10 days have ages that are less than 2 Gyr. Newton et al. (2017) showed that these kinematically younger stars also display increased chromospheric Hα emission while their slowly rotating counterparts show little or no Hα emission. Kiman et al. (2019) showed that fully convective M dwarfs with high levels of Hα activity have smaller velocity dispersions in the W component of their galactic velocity.
A better understanding of the relationship between age, rotation, and magnetic activity for fully convective M dwarfs is also needed to enable an understanding of their attendant planets. Transiting planets around these stars, such as the ones in our sample: GJ 1132b (Berta- Bonfils et al. 2018), LHS 1140bc (Dittmann et al. 2017;Ment et al. 2019), LHS 3844b (Vanderspek et al. 2019), TOI-540b (Ment et al. 2021), may be targeted in the near-term for atmospheric studies with the James Webb Space Telescope and next generation of extremely large telescopes. Although we often find planets orbiting inactive stars that have transitioned to the unsaturated regime, at earlier times these stars resided in the saturated regime and their planets suffered the consequences. It is important that we characterize the duration of the extended radiation environment associated with the saturated regime. It will allow us to reconstruct their stellar history, which is essential to understanding the current atmospheres, or lack thereof of these planets.
In this study we present a relationship between magnetic activity, kinematic age, and rotation for all single, or effectively single, fully convective M dwarfs with masses between 0.1 and 0.3 M that reside within 15 parsecs and were observed during the primary mission of TESS. We deem a star effectively single if it has a binary companion with a separation greater than 63 arcseconds, which corresponds to the width of three TESS pixels. We build upon the study presented in M20 in which we determined relationships between flare rates, rotation, and chromospheric emission among the subset of these stars in the southern ecliptic hemisphere (observed in Year 1 of the TESS Mission) to include the northern ecliptic sample of TESS stars observed during Year 2 of the TESS Mission.In § 2 we describe the stellar sample. In § 3, we describe our spectroscopic observations, our measurement of the equivalent widths of activity indicators, and our measurement of the radial velocities and estimates of the galactic space motions. We describe the photometric observations and measurement of new rotation periods from the MEarth project in § 4 and discuss the rotation period and flare analysis using TESS photometry in § 5. We present the results in § 6 followed by the discussion and conclusion in § 7.

STELLAR SAMPLE
In this paper we will be discussing two samples of stars. Both samples are subsets of the volume-complete sample of 512 M dwarfs with masses between 0.1-0.3 M that reside within 15 parsecs compiled by Winters et al. (2021).
The first sample we will discuss, which we refer to as the northern ecliptic sample are the stars observed during Year 2 of the TESS Mission, and indeed all of these stars fall within the northern ecliptic hemisphere. For these stars we present measurements of the flare rates, rotation periods, and spectroscopic activity indicators. To compile the northern ecliptic sample, we begin with the Winters et al. (2021) sample and retain only those stars observed in TESS sectors 14-26. We used the information provided in Winters et al. (2021) to remove any stars that are known eclipsing binaries, single-lined or doubled-lined spectroscopic binaries, and visual binaries with separations less than the width of three TESS pixels (63 arcseconds). Our northern ecliptic sample comprises 94 stars that are either single, or are separated by their companions by at least 3 TESS pixels. Winters et al. (2021) determined the distances to these stars using parallaxes from the second data release of Gaia (Gaia Collaboration et al. 2018) and their masses using the K s magnitude and the relations presented in Benedict et al. (2016). The uncertainties on the masses range from 4.7% to 14.0%.
The second sample of stars, which we will refer to as the all-sky sample is the union of stars from the northern ecliptic sample discussed above and the sample of 125 stars from the southern ecliptic hemisphere observed in TESS sectors 1-13 presented in M20. In total there are 219 effectively single fully convective M dwarfs that reside within 15 parsecs in the all-sky sample. In Table 1, we list the names, TIC identifiers, coordinates, masses, distances, proper motions, and TESS magnitudes for the all sky sample. Within that Table, we provide additional quantities, including spectroscopic activity indicators, radial velocities, rotation periods, Rossby numbers, flare rates, and spatial velocities in galactic coordinates. These quantities are either derived as described in the following sections, or reproduced from M20.

SPECTROSCOPIC OBSERVATIONS
We used the Tillinghast Reflector Echelle Spectrograph (TRES) located on the 1.5 meter Tillinghast Reflector telescope at the Fred Lawrence Whipple Observatory (FLWO) on Mount Hopkins, Arizona to obtain three or more high-resolution spectra for each target in the all-sky sample that had a declination greater than -15 degrees. We obtained spectra from September of 2016 to January of 2021. TRES covers the wavelength range 3900-9100Å and has a resolution R≈44,000. Exposure times ranged from 120s to 3 × 1200s to reach a signal to noise ratio of 10-40 at 7150Å. We reduced the spectra using the standard TRES pipeline (Buchhave et al. 2010).
For stars in the all-sky sample with declinations below -15 degrees we used the the HIgh ResolutiON (CH-IRON) spectrogragh (Tokovinin et al. 2013) located on the Cerro Tololo Inter-American Observatory (CTIO) SMARTS 1.5 meter Telescope (see M20, for details). CHIRON has wavelength coverage that extends from 4100-8700Å. We used the image slicer mode to achieve a resolution of R≈80,000. We gathered spectra from December of 2017 to March of 2021. We reached a signal to noise of 3-15 per pixel at the wavelength of 7150Å using exposure times that ranged from 180 seconds to 3 × 1800 seconds. We used the standard CHIRON reduction pipeline to extract each spectrum (Tokovinin et al. 2013;Paredes et al. 2021).

Spectroscopic Activity Indicators
In M20, we presented, for all stars in the southern ecliptic portion of the all-sky sample, the equivalent widths (EWs) of the spectroscopic magnetic activity indicators Helium I D3 at 5875.6Å, Hα at 6562.8Å, and one of the three Calcium infrared triplet lines, 8542.09 A. We did not report measurements of the other two members of the Calcium infrared triplet as they fell outside of the CHIRON spectral format. We reproduce these values in Table 1 for completeness.
We used FLWO/TRES to measure the EWs for these same three spectral features for all stars in the northern ecliptic portion of the all-sky sample. We follow the methods described in M20: we use the wavelength ranges for the feature and surrounding continuum regions as listed in Table 2. We take negative EW values to denote emission and state the maximum EW value (least amount of emission) of our multi-epoch observations as we believe this value provides a better representation of the quiescent level of emission. We report these values in Table 1.

Radial Velocities and Galactic Kinematics
We used cross-correlation methods similar to those described in Kurtz & Mink (1998) to measure radial velocities for every star in the all-sky sample. For our template, we used a high-SNR spectrum of Barnard's Star, observed with TRES on UT 2018 July 19 and with CHIRON on UT 2018 April 22. Barnard's Star is a slowly rotating (130.4 days, Benedict et al. 1998) M4.0 dwarf (Kirkpatrick et al. 1991) for which we adopted a Barycentric radial velocity of −110.3 ± 0.5 km s −1 , derived from presently unpublished CfA Digital Speedometer (Latham et al. 2002) measurements taken over 17 years. We saw negligible rotational broadening in our Barnard's Star template, in agreement with the v sin i of 0.07 km s −1 expected from its long photometric rotation period. This is also consistent with the v sin i upper limit of 2 km s −1 reported by Reiners et al. (2018).
We used methods similar to those described in Winters et al. (2018Winters et al. ( , 2020, which used only a single order for each spectrograph (order 41 for TRES and order 44   . We selected these orders from the red portion of the format (so that there was sufficient SNR), avoided the Hα feature, and were reasonably free of telluric lines, such that the peaks of the correlation functions generally exceeded 0.5 in wellexposed spectra.
We calculated the RVs in each order, with uncertainties calculated using the methods described in Bouchy et al. (2001). The RVs in each order were then combined using the inverse variance weighted mean, with the uncertainty being the standard error in the inverse weighted mean. We repeated this for each epoch. The final reported RV is the weighted mean from all available epochs. The final reported uncertainty is the weighted mean uncertainty from all available epochs scaled by 1.65 and 1.8 for TRES and CHIRON respectively to reach agreement with the observed velocity scatter and then added in quadrature with the 0.5 km s −1 RV uncertainty of the Barnard's star template spectrum; this is the RV uncertainty we present in Table 1. For six of our targets that were faint and/or rapidly rotating, the spectra had extremely low SNR (SNR ≤ 2), we use only the echelle order with the highest signal to noise and hence most reliable RVs to compute the cross-correlation; for TRES this was order 41 and for CHIRON this was order 44.
In the cases where there is measurable rotational broadening, we measured the v sin i per target spectrum as follows. We performed two nested grid searches for the maximum peak correlation, h, over a v sin i range of 0 − 100 km s −1 sampled at 1 km s −1 . We calculated a peak correlation against a broadened template spectrum. The template spectrum is broadened by convolving the template with a standard limb darkened rotational broadening kernel (Gray 2005) . We then used parabolic interpolation of the quantity h 2 to obtain a more precise estimate of the peak. To further refine the estimate, we repeated this procedure ± 1 km s −1 about the best value from the first grid sampled at 0.1 km s −1 . We again used parabolic interpolation to obtain the final value of v sin i from the 0.1 km s −1 grid results. We then averaged the v sin i per object and used this value to apply rotational broadening to the template spectrum when deriving the RVs.
In order to compute the spatial velocities of the allsky sample, we use proper motions and parallaxes as well as their associated uncertainties from the second data release of Gaia (Gaia Collaboration et al. 2018). For one star 2MASS J0943-3833, our RV analysis was inconclusive as we could not find a peak in the CCF and the spectrum was very low signal to noise. We remove this star from further kinematic analysis. For the remaining 218 stars, we used the radial velocities, proper motions, and parallaxes to compute the U,V, and W space motions for each star following the definition of Johnson & Soderblom (1987), where U describes the velocity towards the galactic center, V describes the velocity with respect to the rotation of the galaxy where positive is in the direction of galactic rotation, and W describes the vertical velocity with positive denoting in the direction of the North Galactic pole. We measure the velocities relative to the barycenter of the solar system. The uncertainties in U, V, and W are calculated using the method described in Johnson & Soderblom (1987). The dominant source of error in the calculation of the UVW space motion comes from our measured radial velocities. In Table 1, we list the proper motions, parallaxes, radial velocities, and UVW space motions measured relative to the local standard of rest as well as derived quantities defined in the following sections including the total space velocity, which is the magnitude of the combined U lsr , V lsr , and W lsr velocity components, and the probability ratio of each star residing in the thick disk versus the thin disk. To convert the component velocities to LSR, we add the LSR component velocities measured by Schönrich et al. (2010) of 11.10 km s −1 , 12.24 km s −1 , 7.25 km s −1 for U lsr , V lsr , W lsr .

NEW PHOTOMETRIC ROTATION PERIODS FROM MEARTH
Inhomogeneously distributed stellar spots rotating into and out of view can produce a photometric modulation that can be used to measure stellar rotation periods. The MEarth North and South Observatories each consists of eight 40-inch telescopes located atop Mt. Hopkins in Arizona and CTIO in Chile (Nutzman & Charbonneau 2008;Irwin et al. 2015) and have been used extensively since 2008 and 2014 respectively to determine photometric rotation periods for fully convective M dwarfs. For the 219 stars in our all-sky sample, 67 had previously published rotation periods measured from MEarth data (Newton et al. 2016Díez Alonso et al. 2019;Morales et al. 2019;Vanderspek et al. 2019). In M20, we presented 18 additional rotation periods; 17 of these were from MEarth-South and one was from MEarth-North. Here we followed the methods described in Newton et al. (2016Newton et al. ( , 2018 and used in M20 to measure an additional eight rotation periods, all from MEarth-North. We list the new and previously determined MEarth rotation periods and semi-amplitudes of variability in Table 1. In Table 1 we also list two additional rotation period measured from other groundbased data sets. Suárez Mascareño et al. (2016) used data from the All-Sky Automated Survey to determine the rotation period for GJ 54.1 (P = 69.2 days) and Díez Alonso et al. (2019) used data from the Super-WASP telescope to determine the rotation period of LEP 2211+4059 (P = 30.0 days). We independently confirm the rotation period of GJ 54.1 using MEarth data, but do not currently have sufficient data of LEP 2211+4059 to confirm its rotation period.

TESS PHOTOMETRY
We analyzed the TESS light curves of the 94 stars in our northern ecliptic sample (the data from year 2 of TESS ). We largely follow the methods described in M20. We summarize the methods below, and we refer the reader to that paper for the full description. We used the pre-search data conditioning simple aperture photometry (PDCSAP). We used all data with a quality flag equal to zero or 512. Quality flag 512 denotes a point that is an impulsive outlier. Feinstein et al. (2020) and M20 found that points during a flare can be assigned a quality flag of 512, and so we retain them.

Photometric Rotation Periods with TESS
We found seven new rotation periods with TESS ranging from 0.29 to 1.31 days. We used the methods described in Irwin et al. (2011);Newton et al. (2016Newton et al. ( , 2018; Medina et al. (2020) to measure these rotation periods. We also present the rotation period for one star, WT 887, from M20. WT 887 has the shortest measured rotation period (P=3.4h) of any star in the all-sky sample. The v sin i = 26 km s −1 we measure for this object from our CHIRON spectra confirms this rapid rotation. For the 24 stars with previously published rotation periods less than 14 days in the northern ecliptic sample that were initially presented in Newton et al. (2016), we measure their rotation periods using the TESS data to check for consistency. We find consistency between the rotation period we measure with the TESS data and the previously published period for 23 stars. For one star, LP 255-11, P rot = 8.04 days (Díez Alonso et al. 2019) the period determination using the TESS data is inconclusive. We do however confirm its rotation period using our own analysis of the MEarth North data. We confirm the previous findings presented in Newton et al. (2016) and M20 that the semi-amplitude of variability is not correlated with rotation period. In Figure 1, we show the semi-amplitude of variability as a function of rotation period for the 120 stars among the 219 in our allsky sample with measured rotation periods from TESS, MEarth or other ground-based efforts. We present the new rotation periods as well as the semi-amplitudes of variability in Table 1.

Stellar Flares with TESS
To search for flares we used the TESS two-minute cadence PDCSAP light curves. We used the same methodology that was described in M20, which we summarize here for completeness. We removed stellar and instrumental correlated variability using the python package exoplanet   Newton et al. (2016), M20, and this work. The purple circles show the rotation periods measured in this work and M20 using TESS data. The two blue plus symbols represent stars that did not use MEarth or TESS data for the determination of the rotation periods of Gl 54.1 and LEP 2211+4059. Red stars denote known transiting planet hosts. Magenta stars denote members of young moving groups discussed in section 6.6. We see no relationship between the rotation period and the semi-amplitude of variability.
for at least three consecutive three sigma outliers. Data that meet these criteria are deemed a stellar flare.

Flare Energies
We determined the energy of each detected flare following the methods described in M20, which we summarize here. We multiply the luminosity of the star in the TESS bandpass by the equivalent duration (ED) of the flare, which is measured in seconds. Physically the ED is defined as the time that it takes for the star in a non-flaring quiescent state to emit the same energy that was radiated by the flare. We determined the ED by integrating the area (i.e. summing the values) during the flare in the light curve after subtracting the Gaussian process model. As in M20, we assign a 5% uncertainty to our flare energies. This uncertainty arises because of the uncertainty in where the flare ends (see M20 for details on how the uncertainty is determined).
In total we find 1959 flare events in the northern ecliptic sample. These events are a mixture of classical single peaked flare morphologies and multi-component morphologies. In this study we were interested in the total energy of each flare event and thus do not differenti-ate between complex and classical flare morphologies because we simply integrate the total area under the light curve. We found that 72 of the 94 stars in our northern ecliptic sample flare at least once during sectors 14-26. We provide a catalog of the Barycentric Julian date at which the peak of the flare occurred, flare energies, amplitudes, durations, and equivalent durations for the stars in the northern ecliptic sample in Table 3.

Flare Completeness Correction
We computed a completeness function for each star to understand over what energy range of flares our algorithm is sensitive. We followed the methods presented in M20, which we summarize here for completeness. We injected flares at 30 logarithmically spaced energies ranging from 10 28 to 10 33 ergs. We used the flare template described in Davenport et al. (2014) to inject a total of 100 flares at each energy into each light curve by injecting 10 flares into a given light curve 10 times. We assigned each flare a random start time, but we took care to not inject a flare at a time where a real flare was located. The flare template is described by the time of peak flux, the amplitude, and the full width at half maximum of the flare, which is a proxy for the flare duration. As discussed in M20, the Davenport et al. (2014) template does not describe all flare morphologies. For the purposes of this study however, we are concerned with only flare energies. Thus as long as the flaring event shows at least three positive, three sigma outliers our algorithm will detect it, regardless of its morphology. We then determined the fraction of flares recovered for each energy, which we model using the error function.
We show the completeness function for each star in the northern ecliptic sample in Figure 2 as well as the completeness function where the results have been scaled to account for the inverse square law to a common distance of 10 parsecs. Figure 2 demonstrates that distance is the dominant variable in setting flare detection sensitivity, followed by stellar mass. At the 50% completeness level, which is the energy at which our algorithm is able to detect at least 50% of injected flares, our algorithm is able to detect flares with energies above Log 10 (E/ergs) = 30.35 on GJ 1111 (distance = 3.58 parsecs). This is in comparison to the 50% completeness energy of Log 10 (E/ergs) = 31.46 for GJ 1171, which is the most distant star in our sample (distance = 14.98 parsecs). We use log to denote natural log; log base 10 is denoted as log 10 . Lacy et al. (1976) showed that the frequency of flares as a function of energy (known as the flare frequency distribution (FFD)) can be described by,

Flare Frequency Distribution
where Ω is a normalization constant and α is power-law exponent of the FFD. In M20, we found that the fully convective stars in our southern elliptic sample shared the same value of α = 1.98 ± 0.02. This value is consistent with, but more precisely determined than previous flare studies of much smaller samples of fully convective M dwarfs Lurie et al. 2015;Silverberg et al. 2016), or those that used data with a much longer cadence of 30 minutes (Ilin et al. 2019).
We determined α for the all-sly sample following the methods described in M20: first, we used Markov Chain Monte Carlo methods and the python package emcee (Foreman-Mackey et al. 2013) to determine α for 39 stars that showed five or more flares during TESS observations of sectors 14-26. We then combined these with the 38 values from M20 finding the error weighted mean for the all-sky sample to be α = 1.984 ± 0.019. For the remainder of this paper, we will present results from the analysis of this combined sample, the all-sky sample. Figure 3 shows the Z-score, which is the difference between the individual α i and the error-weighted average value α, divided by the individual uncertainties σ i . This figure shows that all stars are consistent with a single value of α = 1.984 ± 0.019.
We then assumed this value for all stars, and proceeded to determine either a flare rate (for stars that had at least one flare), or an upper limit, as described in M20. We chose to evaluate the flare rate at an energy of 3.16 × 10 31 ergs or log 10 (E/ergs) = 31.5 in the TESS bandpass. As in M20, this is the energy at which all stars in our sample are 50% complete or above i.e. the energy at which our algorithm is able to detect at least 50% of the flares we inject in the light curve of the star for which we are the least sensitive to flares in the sample. We denote as R 31.5 the number of flares per day at or above an energy of 3.16 × 10 31 ergs. We list the log R 31.5 values and their uncertainties in Table 1. We emphasize that log R 31.5 is the natural log of the flare rate not log base ten. 6. RESULTS

Flare Rate and Stellar Rotation Period
Rossby number is dimensionless proxy for stellar rotation that is defined as Ro = P rot /τ , where P rot is the stellar rotation period and τ , is the mass dependent convective turnover time. Rossby number is dimensionless . The right panel shows the completeness function for each star in our northern ecliptic sample scaled to a common distance of 10 parsecs to account for the inverse square law. After correcting for distance, we see that stellar mass, a proxy for luminosity, is a secondary factor: We are less sensitive to flares of a given energy around stars of greater mass due to a smaller contrast between the flare and the intrinsic stellar luminosity. We use these curves to correct for this changing sensitivity on a star-by-star basis.  quantity that serves as a proxy for stellar rotation. It is defined as Ro = P rot /τ , where P rot is the stellar rotation period and τ , is the mass-dependent convective turnover time. The convective turnover time has not been measured directly for any star other than the Sun. We used empirical relations from Wright et al. (2018) to compute the convective turnover time for each star. In Figure 4, we show the flare rate as a function of rotation period and Rossby number for the 122 stars in the all-sky sample with a measured rotation period. We found that there are two different regimes of flaring behavior. Stars that have Rossby numbers less than 0.5 show a saturated flare rate equal to log R 31.5 = -1.32 ± 0.06, whereas stars with Rossby numbers greater than 0.5 displayed flare rates that could decrease six orders of magnitude or more (since, for some stars, we only have an upper limit on the flare rate). From Figure 4, we note that stars in the saturated regime appear to show an intrinsic scatter of 0.79 in the value of log R 31.5 about the mean of -1.32.
Instead of a sharp drop in the flare rate for Rossby numbers above 0.5, there exists a region of Rossby numbers above 0.5 for which stars with the same Rossby number can show either a saturated or unsaturated flare rate. This dichotomy in flare rates for a given Rossby number is reminiscent of a finding by Morin et al. (2010) concerning a dichotomy in magnetic field topologies of fully convective M dwarfs. They find that stars with non-axisymmetric fields with quadrupolar and octupo-lar components have a weak global magnetic field compared to stars of the same rotation period with a dipolar asymmetric fields and a strong global magnetic field. However, by construction, due to the limitations of Zeeman Doppler Imaging, the technique used to measure magnetic field topologies, the stars in the Morin et al. (2010) sample were rapidly rotating, with rotation periods less than 1.5 days. The stars we observe in this region all have rotation periods significantly longer that this. Brown (2014) provides a possible explanation by proposing two different regimes of angular momentum loss separated by an abrupt transition. At young ages, the stellar dynamo, which is believed to be responsible for the generation of a magnetic field, is weakly coupled to stellar wind making angular momentum losses minimal. At later times, the dynamo abruptly transitions to be strongly coupled to the stellar wind leading to a time of greater angular momentum loss and the star spins down. Perhaps in this region above Ro = 0.5, we are observing a mix of stars that are in the process of this abrupt transition.
We find that there are four exceptions to this relationship between rotation and flaring: GJ 1187, WT 887, SCR J1855-6914 and 2MASS J23303802-8455189. These stars have a rotation periods of 0.83, 0.143, 1.40, and 6.43 days and log R 31.5 = -6.77 ± 0.94, -4.64 ± 0.84, and upper limits of -10.06, -11.9, respectively. All of these stars also show Hα in emission. We postulate a few reasons for these anomalously low flare rates. One is that these stars may have undetected binary companions, which may lead to incorrect mass estimates erroneously placing them in our sample or perhaps that the flares and rotational modulation are coming from different members of the system. Using a bisector analysis of the spectral line profiles of these four stars (E. Pass private communication) there is no line variation as a function of time that would indicate of a companion that has a flux ratio of at least 0.4 and less than 0.8 and an RV offset of at least 3.0 and less than 8.0 km s −1 (E. Pass private communication). This analysis does not rule out equal mass binaries or ones with an RV offset less than 3.0 km s −1 . However it is unlikely that GJ 1187, WT 887, and SCR-J1855-6914 are unresolved equal mass binaries because they would be over-luminous and we find that their photometric distances and trigonometric distances are in agreement. For 2MASS J23303802-8455189 the photometric distance is greater than the trigonometric distance indicating this star is actually underluminous. Furthermore, we do not believe any of these stars are astrometric binaries as the Gaia astrometric noise, which can be used as indication of an astrometric binary if the value is greater than two, is less than one for all of these stars. We would like to note that WT 887 is the most rapidly rotating star in our sample with a rotation period of 0.11 days and thus may be demonstrating the phenomenon of super-saturation, wherein extremely rapidly rotating stars show reduced levels of magnetic activity (Doyle et al. 2022).
If a binary companion is not involved, it could be the case that there is some aspect of the evolution of the star that is suppressing magnetic activity. To investigate this, we explored the location of these stars on color magnitude diagrams (CMD). We obtained colors for our sample of stars from Gaia DR2. We find, when examining the absolute K band magnitude as a function of G-K, that WT 887 and SCR J1855-6914 are some of the lowest mass stars in the sample. This same CMD shows that GJ 1187 and 2MASS J23303802-845518 have relatively blue colors and thus, may be metal poor which could effect the rotational evolution of the star and thus its activity (Amard et al. 2020). In any case, we do not include GJ 1187, WT 887, SCR J1855-6914, or 2MASS J23303802-8455189 in the calculation of the saturated value of log R 31.5 = -1.32 ± 0.06.

Flare Rate and Spectroscopic Activity Indicators
In this section we explore the relationship between two magnetically induced phenomena; chromospheric emis-sion and stellar flares. In Figure 5, we show the flare rate as a function of the equivalent widths of Hα, He I D 3 , and one of the calcium infrared triplet lines at 8542.09 A. As in M20, we find both a saturated and an unsaturated regime, and in all three cases there exists a range of EWs in which the saturated and unsaturated populations overlap: for Hα this occurs between EWs of -1.0 and 0.0Å. We find Hα to be the most well-separated where for Hα EW values less than -1.0Å, there is a distinct region where the flare rate increases with increasing Hα emission. Because of this, we have focused on exploring the relationship between Hα and the flare rate.
Previous studies have found a saturated relationship between various activity indicators and rotation (Wright et al. 2018;Douglas et al. 2014;Newton et al. 2017). Instead in M20, we found that a model where the flare rate in the saturated regime increases modestly with increased Hα emission was favored over a model where the flare rate showed a constant value in the saturated regime. We examined this particular relationship fur-ther using the piece-wise function presented in M20 to fit the relationship between Hα and R 31.5 , log R 31.5 (x Hα ) =    γ(x Hα − X Hαc ) + C x Hα ≤ X Hαc β(x Hα − X Hαc ) + C x Hα > X Hαc (2) where X Hαc is the critical value, C is the rate at which the two regimes intersect, x Hα is the Hα EW, γ is the slope for EW values below the critical value, and β is the slope above the critical EW value. As we find that there is significant overlap in values of the flare rate for a given Hα EW value between -1.5 and 0.0Å, we fit this function twice at two critical values X Hαc1 and X Hαc2 to reflect this region. We determine the value of X Hαc1 using χ 2 minimization by fitting equation 2 to each Hα EW value and taking the Hα corresponding to the minimum χ 2 as X Hαc1 . We found that X Hαc1 = -0.183Å. We decided to place X Hαc2 at -1.5Å as this is where the lower limit appears to fall.
Using these critical values for X Hαc1 and X Hαc2 , we used emcee to determine the values γ, β, and C for each critical value labeled with subscripts 1 and 2 respectively. We place uniform priors on each of these parameters. For X Hαc1 , we found the value for C 1 is -2.58 ± 0.23, γ 1 is -0.26 ± 0.13, and β 1 is -19.85 ± 1.16. For X Hαc2 we fix β 2 = β 1 and determine the values of C 2 = -1.69 ± 0.09 and γ 2 = -0.11 ± 0.03. The fits are shown Figure 6.

Fraction of Active Fully Convective M dwarfs
We explore the dependence of magnetic activity on stellar mass. To do this we examine the fraction of stars that are magnetically active in three mass bins (0.1 < M ≤ 0.15), (0.15 < M ≤ 0.20), and (0.20 < M ≤ 0.30) where each mass bin has 73, 68, and 78 stars respectively. We define a star being magnetically active if it has a flare rate, Log R 31.5 ≥ -5. We find the fraction of active stars is 33 ± 5%, 53 ± 6%, and 31 ± 5% for the (0.1 < M ≤ 0.15), (0.15 < M ≤ 0.20), and (0.20 < M ≤ 0.30) mass bins respectively. In addition, we also determine the fraction of stars that are active as defined by having an Hα EW < -1.0Å. We find the fraction of Hα active stars is 36 ± 5%, 42 ± 6%, and 25 ± 5% for the (0.1 < M ≤ 0.15), (0.15 < M ≤ 0.20), and (0.20 < M ≤ 0.30) mass bins respectively. In Figure 7, we show the fraction of active stars as defined by the flare rate and Hα EW in each mass bin. The error bars were computed using binomial statistics. We find that the fraction of active stars defined by the flare rate peaks at intermediate mass ranging from 0.15 -0.20 M . When considering the Hα indicator, we find that the activity fraction increases from the highest mass bin to the intermediate, but we find that the there is not a statistically significant decrease from the intermediate mass bin to the lowest mass bin. We find the behavior between the fraction of active stars as defined by Hα as a function of stellar mass to be in agreement with previous studies (e.g Hawley et al. 1996;Gizis et al. 2000;West et al. 2004West et al. , 2015 where the fraction of field stars that are active as defined by having Hα EW < -1Å increases with decreasing stellar mass. However, the fact that the fraction of stars that are in the saturated flare regime does not increase monotonically with decreasing mass is worthy of note.

Kinematics in the Galactic Disk
The Galactic disk contains three populations of stars showing different kinematic properties; the thin, thick disk, and halo (see, e.g., Bensby et al. 2003). Compared to the thin disk, stars in the thick disk show larger velocity dispersion, an increased scale height, and come from an older population. The halo is generally made up of the oldest stars in the Galaxy with high velocity dispersions and a spherical spatial distribution. Fantin et al. (2019) finds that the thin disk has had approximately constant stellar formation for the past 8 Gyr, the thick disk had an epoch of stellar formation that peaked 10 Gyr ago, and the halo had its peak star formation epoch approximately 11 billion years ago. In this study we determine, the ratio of the probability a given star is a member thick disk over the probability it is a member of the halo, and the probability a star is in the thick disk over the probability it is in the thin disk. Newton et al. (2016) showed that fully convective M dwarfs with rotation periods less than 10 days were more likely to be dynamically cold thin disk members than their intermediate or slowly rotating counterparts. We followed the methods described in Bensby et al. (2003) to extend the Newton et al. (2016) study to include updated kinematic information and flare rates. Bensby et al. (2003) estimated disk membership based on the U, V, and W space velocities of each star, the characteristic velocity dispersions σ U lsr , σ V lsr , σ W lsr , asymmetric drift V asym in each component, and the number density of stars present in each component. We used the values presented in Table 1 of Fantin et al. (2019) for σ U lsr = (33, 15, 15) km s −1 , σ V lsr = (40, 32, 28) km s −1 , σ W lsr = (131,106, 85) km s −1 and asymmetric drift V asym = (-12, -85, -226) km s −1 for the (thin disk, thick disk, halo) respectively. We assumed the local number density of stars to be 83.0% for the thin disk, 17.0% for the thick disk, and that 0.1% of stars reside in the halo of the galaxy (Fantin et al. 2019). In

Fraction Active
Flare Active H Active Figure 7. The fraction of active stars as defined by Log R31.5 ≥ to -5 (purple points) and as defined by Hα EW < -1.0Å (maroon triangles) as a function of stellar mass. The purple and maroon points are offset in stellar mass for clarity.
≥ 10, stars that are indeterminate are defined as 0.1 < P(thick)/P(thin) < 10, and stars located in the thin disk as P(thick)/P(thin) ≤ 0.1. We also compute P(thick)/P(halo), where P(thick)/P(halo) > 10 indi-cates a high probability a star is a member of the thick disk. We found that with the exception of two stars LHS 288 and GJ 299 where P(thick)/P(halo) = 7.5 and 7.9 respectively, P(thick)/P(halo) > 100. LHS 288 and GJ 299 are the outliers at P(thick)/P(thin) > 10 6 in Figure  8. We found that 5.05% of stars reside in the thick disk, 30.73% in the indeterminate region, and the remaining 64.22% reside in the thin disk.
In Figure 9 we show the Toomre Diagram which is used to distinguish thin, thick, and halo populations. We see that thin disk stars, and stars for which the assignment is indeterminate with a value P(thick)/P(thin) < 1 are all contained within the V tot = 75 km s −1 constant velocity contour. Furthermore thin disk magnetically active stars with Log R 31.5 > -4 are all confined within the V tot = 50 km s −1 constant velocity contour indicating that they are kinematically cold. We find that indeterminate stars with 0.1 < P(thick)/P(thin) < 10 reside in between V tot = 25-125 km s −1 . The thick disk stars, P(thick)/P(thin) > 10 all reside outside of the V tot = 75 km s −1 indicating that this population of stars is kinetically warmer. Newton et al. (2016) studied a similar sample of stars in the pre-Gaia era and found that the velocity dispersion of the individual velocity components, (U lsr ,V lsr ,W lsr ), and the total space velocity increased as a function of increasing rotation period. We show the individual (U lsr ,V lsr ,W lsr ) space motions as a function of rotation period with a colorbar indicating the flare rate in Figure 10. We found that in general stars with rotation periods exceeding 90 days show a higher velocity dispersion in each component and have log R 31.5 that can be 6 orders of magnitude or more below the saturated value than stars with rotation periods < 90 days. At longer rotation periods, we also see evidence of asymmetric drift in the V component where the velocities become increasingly negative at rotation periods longer than 90 days (Eggen & Iben 1989;West et al. 2015;Newton et al. 2016). In Figure 11 we plot the total space motion as a function of stellar rotation. Here again, we observed that stars with rotation periods that exceed 90 days have low flare rates and a higher total space velocity dispersion.
Previous studies have used a relationship between age and velocity dispersion (Wielen 1977;Schmidt et al. 2007;Reiners & Basri 2009;Faherty et al. 2009;Aumer & Binney 2009;Aumer et al. 2016;Newton et al. 2016;Yu & Liu 2018;Kiman et al. 2019;Lu et al. 2021) to estimate ages of stars in the solar neighborhood. We followed the methods first outlined in Yu & Liu (2018) and described further in Lu et al. (2021), which uses the velocity dispersion in the Z direction σ vz , (towards the galactic north pole) as defined by the Galactocentric coordinate system to determine ages. We use Galpy to transform heliocentric UVW coordinates into the Galactocentric coordinate system (Bovy 2015). We calculate σ vz using the median absolute deviation (MAD) multiplied by 1.48 to account for fact that we assume the underlying distribution of Z component velocities is Gaussian. To compute the uncertainties in σ vz , we perturbed the Z component of the velocity using 1000 Gaussian deviates with the standard deviation set to the error bar of the Z velocity component value. We broke our sample into bins based on rotation period. A rapidly rotating bin with stars showing P rot < 1 day, a moderately rapidly rotating bin 1 < P rot < 10 days, an intermediate rotation period bin 10 < P rot < 90 days, and a slowly rotating bin with P rot > 90 days. We computed σ vz for each bin finding that σ vz = 8.63 ± 1.52 km s −1 , 9.42 ± 2.99 km s −1 , 17.56 ± 4.97 km s −1 , and 29.71 ± 4.93 km s −1 , respectively. Lu et al. (2021) provides the following equation to compute the age based on σ vz , where β = 1.58 ± 0.19 and α = -2.80 ± 0.53. As noted in Lu et al. (2021), Yu & Liu (2018) computed the age velocity relation (AVR) for a sample of metal-poor thick disk stars and a sample of metal-rich thin disk stars, but did not provide the value of the intercept for the metal-rich thin disk sample. As such Lu et al. (2021) re-estimated equation 3 and this is the one we use in this study. We found that the rapidly rotating bin has an age of 1.8 ± 0.5 Gyr and mean log R 31.5 = -1.69 ± 0.07, the moderately rotating bin has an age of 2.1 ± 1.1 Gyr and mean log R 31.5 = -2.11± 0.12 , the intermediate rotators have an age of 5.6 ± 2.7 Gyr and mean log R 31.5 = -6.54 ± 0.15, and the slow rotators have ages of 12.9 ± 3.5 Gyr and mean log R 31.5 = -9.65 ± 0.12. We computed the uncertainty in the age by repeating the age calculation for each σ vz value plus and minus its uncertainty. The AVR presented in Yu & Liu (2018) is not the only relation presented in the literature. Aumer & Binney (2009) andAumer et al. (2016 present AVRs that also use only the vertical velocity dispersion and a power law with a different functional form than Yu & Liu (2018) to describe the relationship between age and velocity dispersion. Wielen (1977) presents a relation using the total velocity dispersion σ tot = σ 2 u + σ 2 v + σ 2 w . We determined ages using the methods described in Wielen (1977) and Aumer et al. (2016) finding consistency between the ages we computed and those computed using other methods. We chose to use the AVR presented in . The probability P that a star resides in the thick disk, P (thick) divided by the probability that a star resides in the thin disk P (thin) as a function of the V lsr . The dashed lines denote the thin disk, P(thick)/P(thin) ≤ 0.1, the indeterminate region defined as 0.1 < P(thick)/P(thin) < 10, and the thick disk defined as P(thick)/P(thin) ≥ 10. The color of the points denotes the flare rate, Log R31.5, value given by the color bar. Red stars denote known transiting planet hosts in our sample. Magenta stars denote members of young moving groups discussed in section 6.6. Yu & Liu (2018) as it incorporates updated kinematic information from the Gaia mission. In Table 4, we present the bins, number of stars in each bin, the mean period, σ vz , and their assigned ages.
For completeness we also determined σ vz and estimate an age for stars with no measured rotation period. We found that σ vz = 26.01 ± 3.01 km s −1 leading to an age of 10.5 ± 1.9 Gyr. The age and σ vz of these stars is consistent with the age of stars with rotation periods exceeding 90 days. Thus we postulate that the majority of the stars without a measured rotation period are old and rotate slowly. This statement is further supported by the fact that all of these stars do not show Hα in emission.
Although we broke up our sample into four distinct rotation period bins, σ vz for the stars with rotation pe-  . The Toomre Diagram which is the defined as U 2 lsr + W 2 lsr as a function of V lsr for all the stars in our sample. The circles represent stars that are likely members of the thin disk, squares represent stars for which the assignment to the thin or thick disk is indeterminate, and star-shaped symbols represent stars that are likely members of the thick disk. The color of the points denotes the flare rate, Log R31.5, value given by the color bar. Where we have broken the colorbar up in discrete regions of the Log R31.5 = greater than -4, between -4 and -7, and less than -7. Lines of constant Vtot are indicated with the dashed grey lines and are in steps of 25 km s −1 . Red stars denote known transiting planet hosts. Magenta stars denote members of young moving groups discussed in section 6.6.  Log R31.5 Figure 11. The total space velocity, Vtot as a function of stellar rotation period. The color of the points denotes the flare rate,Log R31.5, given by the color bar. Red stars denote known transiting planet hosts. Magenta stars denote members of young moving groups discussed in section 6.6.
riods < 1 day and σ vz for stars with rotation periods between 1 and 10 days are statistically indistinguishable and thus their inferred ages are indistinguishable. As such, we combined the rapidly rotating and moderately rotating bins and determined an age of 2.0 ± 1.2 Gyr. With this, we posit that stars can reside in the very active phase associated with the saturated regime for at least 2 Gyr.

σ vz as a function of Stellar Mass
We calculated the vertical velocity dispersion as a function of rotation period for the three mass bins presented in Section 6.3 for stars with rotation periods less than 10 days. We implemented the 10 day rotation period cut-off because we were interested in examining the relationship between age and stellar mass in the saturated regime in particular. We found that σ vz = 8.83 ± 2.34 km s −1 , 11.09 ± 2.60 km s −1 , and 3.11 ± 1.30 km s −1 for the (0.1 < M ≤ 0.15), (0.15 < M ≤ 0.20), and (0.20 < M ≤ 0.30) mass bins respectively. We then determined the ages of each respective mass bin using the relations from Lu et al. (2021). We found ages of 1.9 ± 0.9 Gyr, 2.7 ± 1.0 Gyr, and 0.6 ± 0.3 Gyr for the low, intermediate, and high mass bins respectively. From this, we observe that the stars in our lowest and intermediate mass bin are older than stars in our highest mass bin implying that lower mass stars stay in the saturated regime of magnetic on average longer than high mass stars. In Figure 12, we present the V z component of the velocity as a function of rotation period for the low, intermediate, and high mass bins.

Membership in Young Associations
We used BANYAN Σ (Gagné et al. 2018) to determine whether any stars in our all-sky sample are members of young stellar associations. The BANYAN Σ tool takes as inputs the coordinates, radial velocity, proper motion in right ascension, µ α , and declination, µ δ and parallax of the star. We found the majority of our stars (216) to have the highest probability of being field stars. However, we found that three stars have probabilities greater than 99.99% of being associated with young moving groups. G 7-34 has a 99.997 % probability of being a member of the AB Doradus moving group that has an age of 149 Myr and has already been noted in Bell et al. (2015) as being a member of this moving group. AP Col and G 161-71 have probabilities of 99.996% and 99.994% respectively of being members of the Argus Association which has an age of 40 Myr (Torres et al. 2008;Zuckerman 2019). AP Col has been noted as being a pre-main sequence star and being a likely member of the Argus Association by Riedel et al. (2011) and G 161-71 was previously noted as being a member of Argus by Bartlett et al. (2017). G 161-71 has the greatest amount of Hα in emission with an EW = -11.45Å of all stars in sample. AP Col shows the second highest Hα EW = -8.61Å, and G 7-34 shows the 5th highest with an Hα EW = -7.82Å. All of these stars have flare rates and rotation periods that place them in the saturated regime of Figure 4.

DISCUSSION AND CONCLUSIONS
We used high-resolution spectroscopic data to measure chromospheric magnetic activity indicators, photometry to measure flare rates and rotation periods, and galactic velocities to constrain the ages of a volume complete sample of 219 fully convective M dwarfs with masses between 0.1 and 0.3M that reside within 15 parsecs. The flare rates, chromospheric activity indicators, and rotation periods for 125 of these stars were initially presented in M20. In this study, we add 94 additional stars that were observed in the northern ecliptic hemisphere during the second year of the primary mission of TESS. We measured seven new rotation periods with TESS photometry that ranged from 0.29 to 1.31 days and eight new rotation periods with MEarth North Photometry ranging from 70 to 160 days. We found that fully convective M dwarfs fall largely into two populations: 29% have Hα in emission, a saturated flare rate equal to log R 31.5 = -1.32 ± 0.06, have Rossby numbers less than 0.50 and form a kinematically cold population with an average age of 2.0 ± 1.2 Gyr. The remaining 71% show little to no Hα in emission, have log R 31.5 ≤ -3.00 ± 0.67, rotation periods, when measured, that exceed exceed 90 days, and have a large galactic velocity dispersion with a mean age of 12.9 ± 3.5 Gyr. For 99 of the 152 stars in the second group, the photometric rotation period has not been determined, we expect all of these stars to have rotation periods greater of 90 days based on their activity level and kinematics.
We found that the stars in our all-sky sample show a common power-law exponent for the flare frequency distribution of α = 1.984 ± 0.019. This is consistent with the value of α ≈ 2 that has been found in previous spaced based studies using Kepler, K2, or TESS for smaller samples of fully convective M dwarfs (Haw-ley et al. 2014;Silverberg et al. 2016;Davenport et al. 2020;Ilin et al. 2019Ilin et al. , 2021, see M20 for an overview of those studies) with our values being more precise than in previous works. However, Feinstein et al. (2020) in a study of stars with temperatures below 4000 Kelvin in young moving groups with ages up to 750 Myr found that α = 1.58 ± 0.03. This disagreement brings up the question whether α may change as a function of stellar age. For the three stars that are probable members of young moving groups with ages less than 150 Myr: G 7-34, G 161-71, and AP Col we found α = 1.90 ± 0.12, 1.86 ± 0.10, and 1.89 ± 0.13. Ilin et al. (2019) and Ilin et al. (2021) in studies of flare stars in different clusters with known ages observed by K2 found of that α appears to be approximately 2 with no dependence on age or spectral type. We found using the methodology and software STELLA presented in Feinstein et al.
(2020) on G 7-34, G 161-71, and AP Col that many spurious lower-energy flares were being detected. The detection of these spurious signals could account for the discrepant value of α that is obtained in that study. Furthermore, we found no indication of a change in α for stars with different kinematic ages. An α value of 2 or greater means means that low-energy flares occur frequently enough that the energy they deposit to the stellar corona may account for observed coronal temperatures ≈ 10 6 Kelvin. Furthermore, this common value of α points to a common mechanism for the production of flares on fully convective M dwarfs and perhaps all stars with convective envelopes.
We compute the mean ages for groups of stars in different rotation period bins. We find that stars with rotation periods less than 10 day have ages of 2.0 ± 1.2 Gyr, stars with 10 < P < 90 days have ages of 5.6 ± 2.7 Gyr, and stars with rotation periods greater than 90 days have an age of 12.9 ± 3.5 Gyr. We find that stars with a saturated flare rate and Rossby number less than 0.5 days likely belong to a younger, kinematically colder population. Through our kinematic analysis, we found that two stars in our sample G 161-71 and AP Col are members of the young stellar association Argus with an age of 40 Myrs and one star, G 7-34 is a member of the AB Doradus Moving Group with an age of 149 Myrs. We found that these stars are distinct in that they have some of the highest Hα emission with EWs of -8.61Å, -11.45Å, and -7.82Å for AP Col, G 161-71, and G 7-34 respectively. Medina et al. (2022) recently showed that G 7-34 differed from the other 9 active stars in that sample in that its Hα emission varied in phase with the stellar rotation period of the star. This in phase variation indicated that the dominant source of Hα variability on G 7-34 originates from constant emission from fixed magnetic structures, such as starspots and plage. This is in contrast to the other 9 active, but older stars in their sample which showed no correlation between Hα variability and rotational phase; the dominant source of Hα variability in these stars is likely due to stellar flares. It is not yet clear why the dominant source of Hα variability on this young star is different than older active stars, but probing the Hα variability timescale of G 161-71 and AP Col may provide insight into how the interplay between emission from magnetic structures and flares evolves with time.

At What Age Do Fully Convective Stars Spin
Down?
Our results, and those of others concerning coronal xray emission (Stauffer & Hartmann 1987;Pizzolato et al. 2003;Wright et al. 2018), Hα emission (Douglas et al. 2014;Newton et al. 2017), and magnetic field strength (Reiners et al. 2022) show that fully convective M dwarfs fall largely into two categories: rapidly rotating stars that show a saturated magnetic activity that is independent of rotation up to a critical value, and stars that fall into the unsaturated activity-rotation regime. Theoretical studies posit that stars occupying the saturated regime of magnetic activity have complex multipolar magnetic field topologies, few open field lines (Garraffo et al. 2018), and a dynamo that is weakly coupled to the stellar wind (Brown 2014). These theoretical studies suggest that stars occupying the unsaturated regime have dipolar fields, many open field lines, and have a dynamo that is strongly coupled to the stellar wind (Garraffo et al. 2018;Brown 2014). A recent observational study by Reiners et al. (2022) found that saturation of the average surface magnetic field is due to a saturation of the magnetic dynamo rather than a saturation of the stellar surface by magnetic features. They suggest that the surface magnetic field generated is limited by the available kinetic energy and thus can increase only until it reaches the kinetic field. It is currently unclear what role, if any, the topology of the field plays. Furthermore, the duration of the saturated phase and at what age the transition from one regime to the other takes place is currently unknown.
Here we estimate, using a very simplified model, on an age that the transition may occur based on our results and the picture put forth by Brown (2014). As discussed previously, Fantin et al. (2019) showed that the star formation in the thin disk has been approximately constant for past 8 Gyr with a stark drop off in the star formation rate at older ages. We show that in Figure 8 that the majority of our stars are members of the thin disk. As such, we assume that stellar formation of the stars in the all-sky sample has been constant for the past 8 Gyr. If we take these assumptions to be true, coupled with the findings that 29% of stars reside in the saturated regime and 71% in the unsaturated regime, we posit that the age at which stars spontaneously transition from the saturated to unsaturated regimes occurs at approximately 2.4 ± 0.3 Gyr. This estimation is consistent with the kinematic ages we measured in this study where we find that stars with rotation periods less than 10 days have ages of 2.0 ± 1.2 Gyr, stars with 10 < P rot < 90 days have ages around 5.6 ± 2.7 Gyr, and stars with rotation periods greater than 90 days have an age of 12.9 ± 3.5 Gyr. The estimate of 2.4 ± 0.3 Gyr fits in where expected, between the rotation period bin containing stars with rotation periods between 10 and 90 days. We do however, note that this is a lower limit as some of our stars do reside in the thick disk where the stellar formation history has not been constant. We also find that when we examine the mass dependence of the transition for stars in the saturated regime with Prot < 10 day, we find that the average age of stars with Prot < 10 days increases from 0.6 ± 0.3 Gyr for stars with masses between 0.2-0.3M to 2.3 ± 1.3 Gyr for stars with masses between 0.1-0.2 M .
With the exception of planets around active M stars like K2-25 b, which is not in our sample, (Gaidos et al. 2020) and TOI-540 b (Ment et al. 2021) which is in our sample, the planets we discover largely orbit inactive stars that reside in the unsaturated regime. However if it is true that these stars sustain high levels of extreme ultra-violet and x-ray radiation associated with the saturated regime for typically 2.4 Gyr, this would have a profound effect on the atmospheres of planets orbiting this stars. These atmospheres are those that will be observable in the near-term with the next generation of space based telescopes and ground based extremely large telescopes. In this study, we characterized the activity in the saturated regime and the age the transition occurs. These constraints are pertinent to the planets orbiting fully convective stars; in order to fully understand the results of upcoming atmospheric studies, the stellar radiation environment and its history must be reconstructed.