The Ca II H&K Rotation-Activity Relation in 53 mid-to-late type M-Dwarfs

In the canonical theory of stellar magnetic dynamo, the tachocline in partially convective stars serves to arrange small-scale fields, generated by stochastic movement of plasma into a coherent large-scale field. Mid-to-late M-dwarfs, which are fully convective, show more magnetic activity than classical magnetic dymano theory predicts. However, mid-to-late M-dwarfs show tight correlations between rotation and magnetic activity, consistent with elements of classical dynamo theory. We use data from Magellan Inamori Kyocera Echelle (MIKE) Spectrograph to detail the relation between Ca II H\&K flux and rotation period for these low-mass stars. We measure $R'_{HK}$ values for 53 spectroscopically identified M-dwarfs selected from the MEarth survey; these stars span spectral classes from M5.0 to M3.5 and have rotation periods ranging from hours to months. We present the rotation--activity relationship as traced through these data. We find power law and saturated regimes consistent to within one sigma of previously published results and observe a mass dependence in $R'_{HK}$.


INTRODUCTION
M-dwarfs are the most numerous stars in our galaxy. Some planet search campaigns have focused on M-dwarfs due to the relative ease of detecting small planets in their habitable zones (e.g. Nutzman & Charbonneau 2008); however, spun-up M-dwarfs are more magnetically active when compared to larger and hotter stars (Hawley & Pettersen 1991;Delfosse et al. 1998;Schmidt et al. 2014). The increase in activity may accelerate the stripping of an orbiting planet's atmosphere (e.g. Owen & Mohanty 2016), and may dramatically impact habitability (Shields et al. 2016). Therefore, it is essential to understand the magnetic activity of M-dwarfs in order to constrain the potential habitability and history of the planets that orbit them. Additionally, rotation and activity may impact the detectability of hosted planets (e.g. Robertson et al. 2014;Newton et al. 2016;Vanderburg et al. 2016).
Robust theories explaining the origin of solar-like magnetic fields exist and have proven extensible to other regions of the main sequence (Charbonneau 2014). The classical αΩ dynamo relies on differential rotation between layers of a star to stretch a seed poloidal field into a toroidal field (Parker 1955;Cameron et al. 2017). Magnetic buoyancy causes the toroidal field to rise through the star. During this rise, turbulent helical stretching converts the toroidal field back into a poloidal field (Parker 1955). Seed fields may originate from the stochastic movement of charged particles within a star's atmospheres.
In non-fully convective stars the initial conversion of the toroidal field to a poloidal field is believed to take place at the interface layer between the radiative and convective regions of a star -the tachocline (Noyes et al. 1984;Tomczyk et al. 1996;Dikpati & Charbonneau 1999). The tachocline has two key properties that allow it to play an important role in solar type magnetic dynamos: 1), there are high shear stresses, which have been confirmed by astroseismology , and 2), the density stratification between the radiative and convective zones serves to "hold" the newly generated toroidal fields at the tachocline for an extended time. Over this time, the fields build in strength significantly more than they would otherwise (Parker 1975). This theory does not trivially extend to mid-late M-dwarfs, as they are believed to be fully convective and consequently do not contain a tachocline (Chabrier & Baraffe 1997). Moreover, fully convective M-dwarfs are not generally expected to exhibit internal differential rotation (e.g. Barnes et al. 2004Barnes et al. , 2005, though, some models do produce it (Gastine et al. 2014).
Currently, there is no single accepted process that serves to build and maintain fully convective M-dwarf magnetic fields in the same way that the α and Ω processes are presently accepted in solar magnetic dynamo theory. Three-dimensional magneto anelastic hydrodynamical simulations have demonstrated that local fields generated by convective currents can self organize into large scale dipolar fields. These models indicate that for a fully convective star to sustain a magnetic field it must have a high degree of density stratification -density contrasts greater than 20 at the tachocline -and a sufficiently large magnetic Reynolds number 1 .
An empirical relation between the rotation rate and the level of magnetic activity has been demonstrated in late-type stars (Skumanich 1972;Pallavicini et al. 1981). This is believed to be a result of faster rotating stars exhibiting excess non-thermal emission from the upper chromosphere or corona when compared to their slower rotating counterparts. This excess emission is due to magnetic heating of the upper atmosphere, driven by the underlying stellar dynamo. The faster a star rotates, up to some saturation threshold, the more such emission is expected. However, the dynamo process is not dependent solely on rotation; rather, it depends on whether the contribution from the rotational period (P rot ) or convective motion -parameterized by the convective overturn time scale (τ c ) -dominates the motion of a charge packet within a star. Therefore, the Rossby Number (Ro = P rot /τ c ) is often used in place of the rotational period as it accounts for both.
The rotation-activity relation was first discovered using the ratio of X-ray luminosity to bolometric luminosity (L X /L bol ) (Pallavicini et al. 1981) and was later demonstrated to be a more general phenomenon, observable through other activity tracers, such as Ca II H&K emission (Vilhu 1984). This relation has a number of important structural elements. Noyes et al. (1984) showed that magnetic activity as a function of Rossby Number is well modeled as a piecewise power law relation including a saturated and non-saturated regime. In the saturated regime, magnetic activity is invariant to changes in Rossby Number; in the non-saturated regime, activity decreases as Rossby Number increases. The tran-sition between the saturated and non-saturated regions occurs at Ro ∼ 0.1 (e.g. Wright et al. 2011). Recent evidence may suggest that, instead of an unsaturated region where activity is fully invariant to rotational period, activity is more weakly, but still positively, correlated with rotation rate (Mamajek & Hillenbrand 2008;Reiners et al. 2014;Lehtinen et al. 2020;Magaudda et al. 2020).
Previous studies of the Ca II H&K rotation-activity relation (e.g. Vaughan et al. 1981;Suárez Mascareño et al. 2015;Astudillo-Defru et al. 2017;Houdebine et al. 2017) have focused on on spectral ranges which both extend much earlier than M-dwarfs and which do not fully probe late M-dwarfs. Other studies have relied on v sin(i) measurements (e.g. Browning et al. 2010;Houdebine et al. 2017), which are not sensitive to the long rotation periods reached by slowly rotating, inactive mid-tolate type M dwarfs (70-150 days: Newton et al. 2016). Therefore, these studies can present only coarse constraints on the rotation activity relation in the fully convective regime. The sample we present in this paper is focused on mid-to-late type M dwarfs, with photometrically measured rotational periods, while maintaining of order the same number of targets as previous studies. Consequently, we provide much finer constraints on the rotation-activity relation in this regime.
We present a high resolution spectroscopic study of 53 mid-late M-dwarfs. We measure Ca II H&K strengths, quantified through the R HK metric, which is a bolometric flux normalized version of the Mount Wilson S-index. These activity tracers are then used in concert with photometrically determined rotational periods, compiled by Newton et al. (2017), to generate a rotation-activity relation for our sample. This paper is organized as follows: Section 2 provides an overview of the observations and data reduction, Section 3 details the analysis of our data, and Section 4 presents our results and how they fit within the literature.

OBSERVATIONS & DATA REDUCTION
We initially selected a sample of 55 mid-late M-dwarfs from targets of the MEarth survey (Berta et al. 2012) to observe. Targets were selected based on high proper motions and availability of a previously measured photometric rotation period, or an expectation of a measurement based on data available from MEarth-South at the time. These rotational periods were derived photometrically (e.g. Newton et al. 2016;Mann et al. 2016;Medina et al. 2020). For star 2MASS J06022261-2019447, which was categorized as an "uncertain detection" from MEarth photometry by Newton et al. (2018), including new data from MEarth DR10 we find a period of 95 days. This value was determined following similar methodology to Irwin et al. (2011) and Newton et al. (2016Newton et al. ( , 2018, and is close to the reported candidate period of 116 days. References for all periods are provided in the machine readable version of Table 1. High resolution spectra were collected from March to October 2017 using the Magellan Inamori Kyocera Echelle (MIKE) spectrograh on the 6.5 meter Magellan 2 telescope at the Las Campanas Observatory in Chile. MIKE is a high resolution double echelle spectrograph with blue and red arms. Respectively, these cover wavelengths from 3350 -5000Å and 4900-9500Å (Bernstein et al. 2003). We collected data using a 0.75x5.00" slit resulting in a resolving power of 32700. Each science target was observed an average of four times with mean integration times per observation ranging from 53.3 to 1500 seconds. Ca II H&K lines were observed over a wide range of signa-to-noise ratios, from ∼ 5 up to ∼ 240 with mean and median values of 68 and 61 respectively.
We use the CarPy pipeline (Kelson et al. 2000;Kelson 2003) to reduce our blue arm spectra. CarPy's data products are wavelength calibrated, blaze corrected, and background subtracted spectra comprising 36 orders. We shift all resultant target spectra into the rest frame by cross correlating against a velocity template spectrum. For the velocity template we use an observation of Proxima Centari in our sample. This spectrum's velocity is both barycentrically corrected, using astropy's SkyCoord module (Astropy Collaboration et al. 2018), and corrected for Proxima Centari's measured radial velocity, -22.4 km s −1 (Torres et al. 2006). Each echelle order of every other target observation is cross correlated against the corresponding order in the template spectra using specutils template correlate function (Earl et al. 2021). Velocity offsets for each order are inferred from a Gaussian fit to the correlation vs. velocity lag function. For each target, we apply a three sigma clip to list of echelle order velocities, visually verifying this clip removed low S/N orders. We take the mean of the sigma-clipped velocities Finally, each wavelength bin is shifted according to its measured velocity.
Ultimately, two targets (2MASS J16570570-0420559 and 2MASS J04102815-5336078) had S/N ratios around the Ca II H&K lines which were too low to be of use, reducing the number of R' HK measurement we can make from 55 to 53.
Since the early 1960s, the Calcium Fraunhaufer lines have been used as chromospheric activity tracers (Wilson 1963). Ca II H&K lines are observed as a combination of a broad absorption feature originating in the upper photosphere along with a narrow emission feature from non-thermal heating of the upper chromosphere (Catalano & Marilli 1983). Specifically, the ratio between emission in the Ca II H&K lines and flux contributed from the photosphere is used to define an activity metric known as the S-index (Wilson 1968). The S-index increases with increasing magnetic activity. The S-index is defined as where f H and f K are the integrated flux over triangular passbands with a full width at half maximum of 1.09Å centered at 3968.47Å and 3933.66Å, respectively. The values of f V and f R are integrated, top hat, broadband regions. They approximate the continuum ( Figure 1) and are centered at 3901Å and 4001Å respectively, with widths of 20Å each. Finally, α is a scaling factor with α = 2.4. Following the procedure outlined in Lovis et al. (2011) we use the mean flux per wavelength interval,f i , as opposed to the integrated flux over each passband when computing the S-index. This means that for each passband, i, with a blue most wavelength λ b,i and a red most wavelength λ r,i ,f i is the summation of the product of flux (f ) and weight (w i ) over the passband.
where w i represents the triangular passband for f H & f K and the tophat for Additionally, the spectrograph used at Mount Wilson during the development of the S-index exposed the H & K lines for eight times longer than the continuum of the spectra. Therefore, for a modern instrument that exposes the entire sensor simultaneously, there will be 8 times less flux in the Ca II H&K passbands than the continuum passbands than for historical observations. This additional flux is accounted for by defining a new constant α H , defined as: Therefore, S-indices are calculated here not based on the historical definition given in Equation 1; rather, the slightly modified version: The S-index may be used to make meaningful comparisons between stars of similar spectral class; however, it does not account for variations in photospheric flux and is therefore inadequate for making comparisons between stars of different spectral classes. The R HK index (Middelkoop 1982) is a transformation of the S-index intended to remove the contribution of the photosphere. R HK introduces a bolometric correction factor, C cf , developed by Middelkoop (1982) and later improved upon by Rutten (1984). Calibrations of C cf have focused on FGK-type stars using broad band color indices, predominately B-V. However, these FGK-type solutions do not extend to later type stars easily as many midlate M-dwarfs lack B-V photometry. Consequently, C cf based on B-V colors were never calibrated for M-dwarfs as many M-dwarfs lack B and V photometry. Suárez Mascareño et al. (2016) provided the first C cf calibrations for M-dwarfs using the more appropriate color index of V − K. The calibration was later extended by Astudillo-Defru et al. (2017), which we adopt here.
Generally R HK is defined as where K is a factor to scale surface fluxes of arbitrary units into physical units; the current best value for K is taken from Hall et al. (2007), K = 1.07 × 10 6 erg cm −2 s −1 . S phot is the photospheric contribution to the S-index; in the spectra this manifests as the broad absorption feature wherein the narrow Ca II H&K emission resides. σ is the Stephan-Boltzmann constant. If we define then we may write R HK as We use the color calibrated coefficients for log 10 (C cf ) and log 10 (R phot ) presented in Table 1 of Astudillo-Defru et al. (2017).
We estimate the uncertainty of R HK as the standard deviation of a distribution of R HK measurements from 5000 Monte Carlo tests. For each science target we offset the flux value at each wavelength bin by an amount sampled from a normal distribution. The standard deviation of this normal distribution is equal to the estimated error at each wavelength bin. These errors are calculated at reduction time by the pipeline. The R HK uncertainty varies drastically with signal-to-noise; targets with signal-to-noise ratios ∼ 5 have typical uncertainties of a few percent whereas targets with signal-tonoise ratios ∼ 100 have typical uncertainties of a few tenth of a percent.

Rotation and Rossby Number
The goal of this work is to constrain the rotation activity relation; therefore, in addition to the measured R HK value, we also need the rotation of the star. As mentioned, one of the selection criteria for targets was that their rotation periods were already measured; however, ultimately 6 of the 53 targets with acceptable S/N did not have well constrained rotational periods. We therefore only use the remaining 47 targets to fit the rotation-activity relation.
In order to make the most meaningful comparison possible we transform rotation period into Rossby Number . This transformation was done using the convective overturn timescale, τ c , such that the Rossby Number, Ro = P rot /τ c . To first order τ c can be approximated as 70 days for fully-convective M-dwarfs (Pizzolato et al. 2000). However, Wright et al. (2018) Equation (5) presents an empirically calibrated expression for τ c . This calibration is derived by fitting the convective overturn timescale as a function of color index, in order to minimize the horizontal offset between stars of different mass in the rotation-activity relationship. The calibration from Wright et al. (2018) that we use to find convective overturn timescales and subsequently Rossby numbers is: We adopt symmetric errors for the parameters of Equation 8 equal to the larger of the two anti-symmetric errors presented in Wright et al. (2018) Equation 5.

ROTATION-ACTIVITY RELATION
We show our rotation-activity relation in Figures 2 & 3. Note that errors are shown in both figures; however, they render smaller than the data point size. Ca II H&K is also known to be time variable (e.g. Baroch et al. 2020;Perdelwitz et al. 2021), which is not captured in our single-epoch data. There is one target cut off by the domain of this graph, 2MASS J10252645+0512391. This target has a measured vsini of 59.5 ± 2.1 km s −1 (Kesseli et al. 2018) and is therefore quite rotationally broadened, which is known to affect R HK measurements (Schröder et al. 2009, figure 8). The data used to generate this figure is given in Table 1. Table 1 includes uncertainties, the R' HK measurements for stars which did not have photometrically derived rotational periods in MEarth, and data for 2MASS J10252645+0512391 We find a rotation activity relationship qualitatively similar to that presented in Astudillo-Defru et al. (2017). Our rotation activity relationship exhibits both the expected saturated and unsaturated regimes -the flat region at Ro < Ro s and the sloped region at Ro ≥ Ro s respectively. We fit the rotation activity relation given in Equation 9 to our data using Markov Chain Monte Carlo (MCMC), implemented in pymc (Salvatier et al. 2016 Ro s is the Rossby number cutoff between the saturated and unsaturated regime. R s is the maximum, saturated, value of R HK and k is the index of the power law when Ro ≥ Ro s . Due to the issues measuring R HK for high vsini targets discussed above, we exclude 2MASS J10252645+0512391 from this fit. All logarithms are base ten unless another base is explicitly given. We find best fit parameters with one σ errors: Note that R HK for one of six of these targets (2MASS J15165576-0037116) is consistent to within 1σ of the saturated value; therefore, the reported Ro for this target should only be taken as an upper bound. The remaining five targets have measured R HK values consistent with the unsaturated regime. Estimated periods are consistent with previous constraints. Of the six stars, two were listed as non-detections in Newton et al. (2018), and the remaining four as uncertain (possible) detections. Of the four classed as uncertain, 2MASS 12384914-3822527 and 2MASS 19204795-4533283 have candidate periods > 100 days and non-detections of H-alpha emission (Hawley et al. 1996). These two stars and the two non-detections have Ca II H&K activity levels suggesting very long periods. 2MASS 13464102-5830117 has a candidate period      2017) and Mamajek & Hillenbrand (2008). Note both that Mamajek & Hillenbrand (2008) focuses their work on earlier spectral classes and fits the rotation activity relation in linear space.
of 45 days, and 2MASS 15165576-0037116 of 0.8 days, both consistent with their higher levels of Ca II H&K emission.
As a test of the proposed weak correlation between activity and rotation in the "saturated" regime seen in some works (Mamajek & Hillenbrand 2008;Reiners et al. 2014;Lehtinen et al. 2020;Medina et al. 2020)though not in others (Wright et al. 2011;Núñez et al. 2015;Newton et al. 2017) -we fit a second model whose power law index is allowed to vary at Ro < Ro s . We find a saturated regime power law index of −0.052 ± 0.117, consistent with 0 to within 1σ. Moreover, all other parameter for this model are consistent to within one σ of the nominal parameters for the model where the index is constrained to 0 below Ro = Ro s . We can constrain the slope in the saturated regime to be between -0.363 and 0.259 at the 3σ confidence level. Ultimately, we adopt the most standard activity interpretation, a fullysaturated regime at Ro < Ro s .
We investigate whether our lack of detection of a slope for Ro < Ro s is due to the limited number of observations in that region when compared to other works (e.g. Medina et al. 2020, 93 targets Ro < Ro s ) through injection and recovery tests. We inject, fake, rotationactivity measurements into the saturated regime with an a priori slope of -0.13 -the same as in Medina et al.. These fake data are given a standard deviation equal to the standard deviation of our residuals (12%). We preform the same MCMC model fitting to this new data set as was done with the original dataset multiple times, each with progressively more injected data, until we can detect the injected slope to the three sigma confidence level. Ultimately, we need more than 65 data points -43 more than we observed in the saturated regimeto consistently recover this slope. Therefore, given the spread of our data we cannot detect slopes on the order of what has previously been reported in the literature.
We observe a gap in rotational period over a comparable range to the one presented in Newton et al. (2016) Figure 2. Namely, that M-dwarfs are preferentially observed as either fast or slow rotators, with a seeming lack of stars existing at mid rotational periods. This period gap manifests in the Rossby Number and can be seen in Figure 3 as a lack of our targets near to the knee-point in the fit. This period gap likely corresponds to that seen by Browning et al. (2010), who found a paucity of M dwarfs at intermediate activity levels in Ca II H&K and note the similarity to the Vaughn-Preston gap established in higher mass stars (Vaughan & Preston 1980). Magaudda et al. (2020) also identify a double-gap in xray activity for stars in the unsaturated regime; it is not clear that the gap we see is related. As a consequence of this period gap, there exists a degeneracy in our data between moving the knee-point and allowing the activity level to vary in the saturated regime. In the following, we adopt the model of a fully saturated regime.
We wish to compare our best fit parameters to those derived in Astudillo-Defru et al. (2017); however, the authors of that paper do not fit the knee-point of the rotation-activity relation. They select the canonical value for the rotational period separating the saturated regime from the unsaturated regime (P rot,s = 10 days) and use a fixed convective overturn timescale (τ c = 70 days). To make our comparison more meaningful we use the P rot and V − K colors presented in Astudillo-Defru et al. (2017) to re-derive Ro values using τ c (Wright et al.  Table 3 and fitting the same piecewise power law as before, we find best fit parameters of Ro s = 0.17 ± 0.04, log(R s ) = −4.140 ± 0.067, and k = −1.43 ± 0.21. Compared to the best fit parameters for our data, Ro s and the unsaturated regime's index, k, are consistent to within one sigma, while the saturated value, R s , differs. The mass ranges of our respective samples explain the differences in saturation values between our work and that of Astudillo-Defru et al. (2017). Our work focuses on mid-to-late M-dwarfs and includes no stars above a mass of 0.5 M ( Figure 5). The strength of Ca II H&K emission is known to decrease as stellar mass decreases (Schrijver & Rutten 1987;Rauscher & Marcy 2006;Houdebine et al. 2017). As Rauscher & Marcy (2006) note, this is the opposite as the trend seen in Halpha; the latter primarily reflects the increasing length of time that lower M dwarfs remain active and rapidly rotating (West et al. 2015;Newton et al. 2016).
A mass dependence can be seen in Figure 10 in Astudillo-Defru et al. (2017), consistent with expectations from the literature. If we clip the data from Astudillo-Defru et al. (2017) Table 3 to the same mass range as our data-set (M * < 0.5M ) and fit the same function as above, we find that all best fit parameters are consistent to within one sigma between the two datasets.
We also compare our best fit Ro s to both those derived in Newton et al. (2017) using H α as an activity measure and those derived in (Wright et al. 2018;Magaudda et al. 2020) using L X /L bol as an activity measure. Works using L X /L bol identify a similar, yet not consistent to within one sigma result for Ro s ; while, the value of k we find here is consistent between all four works. Therefore, we find similar results not only to other work using the same activity tracer, but also a power-law slope that is consistent with work using different tracers.

CONCLUSIONS
In this work we have approximately doubled the number of M-dwarfs with both empirically measured R HK with M * < 0.5M . This has enabled us to more precisely constrain the rotation-activity relation. This relationship is consistent with other measurements using R HK , and L X /L bol ; our data does not require a slope in the saturated regime. Finally, we identify a mass dependence in the activity level of the saturated regime, consistent with trends seen in more massive stars in previous works. This work has made use of the NASA Astrophysical Data System. We thank Dan Kelson for assistance with the CarPy data reduction suite. Additionally, we would like to thank Zach Berta-Thompson, Brian Chaboyer, and John Thorstensen for their support. We acknowledge support for an NSF AAPF awarded to Elisabeth R. Newton (award # 1602597). This material is based upon work supported by the National Science Foundation under grant AST-1616624, and by a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. N.M. was supported by the National Science Foundation Graduate Research Fellowship Program, under NSF grant No. DGE1745303. N.M. thanks the LSSTC Data Science Fellowship Program; his time as a Fellow has greatly benefited this work.