[a]Hiroshi Ohki

# Hybrid stochastic method for the tensor renormalization group

Erika Arai  Masaaki Tomii

###### Abstract

We propose a hybrid stochastic method for the tensor renormalization group
(TRG) approach. TRG is known as a powerful tool to study the many-body systems
and quantum field theory on the lattice. It is based on a low-rank
approximation of the tensor using the truncated singular value decomposition
(SVD), whose computational cost increases as the bond dimension increases, so
that efficient cost reduction techniques are highly demanded. We use noise
vectors for the low-rank approximation with the truncated SVD, by which the
truncation error is replaced with a statistical error due to noise, and an
error estimation could be improved. We test this method in the classical Ising
model and observe a better accuracy than TRG. We also discuss a cross
contamination issue in a multiple use of the same noise vectors, and to remove
this systematic error we consider position-dependent noise vectors.

##  1 Introduction

The tensor network is a powerful tool to study the many-body systems and
quantum field theory on the lattice. In a tensor network representation the
partition functions can be obtained by contracting the tensor products. One
efficient and simple algorithm for a tensor contraction is the tensor
renormalization group (TRG) approach [1]. Since TRG is a fully deterministic
algorithm, there is no sign problem and this method can also be applied to
systems with a complex action such as finite density QCD, the real-time
dynamics that are not easily accessible by the Monte Carlo methods.

TRG is a very efficient numerical method of a coarse graining, but contracting
the tensor products is still the most time-consuming part, whose computational
cost scales as e.g., $\mathcal{O}(D_{\rm cut}^{6})$ in the 2-D Ising model,
where $D_{\rm cut}$ is the truncated bond dimension used in the singular value
decomposition (SVD). The computational cost significantly increases as the
bond dimension increases, so that efficient cost reduction techniques are
highly demanded. In addition, there exists a systematic error, which arises
from the truncation in the SVD. Since an estimation of systematic errors
becomes difficult when the system becomes more complicated, it is important to
investigate a possibility of an alternative approach to improve the error
estimation that can be easily implemented on tensor networks. It is also known
that a short-distance effect remains after coarse graining [2], and the
resulting coarse-grained tensor does not obey a correct RG transformation. To
obtain a physically correct RG flow and to reduce the truncation error, other
renormalization methods have been proposed, such as Loop-TNR [3], TNR [4],
GILT [5]. Also several truncation techniques such as a randomized SVD [6], a
projective truncation method [7], have been proposed, which can reduce the
computation cost as $\mathcal{O}(D_{\rm cut}^{5})$. We note, however, that a
low-rank approximation of the tensor decompositions is commonly used in these
methods, so that there remains a systematic uncertainty due to the truncation
error.

In this work, we propose a new approach based on a stochastic calculation of
the low-rank approximation of a tensor decomposition. We use random noise
vectors in combining with the truncated SVD, in which the truncated bond
dimension is given as the sum of the number of singular values and the noise
vectors. The truncation error can be replaced with an error due to the noise,
and a stochastic determination of physical quantities is possible. A
feasibility study of this method is shown for the 2-D Ising model, where we
employ this new truncation method in the same coarse graining process as in
TRG. We study the noise structure dependence of the free energy and a noise
contamination issue which is found to be the only source of the systematic
uncertainty. We also discuss a possibility to completely remove the systematic
error by employing position-dependent noise vectors.

##  2 Rank truncation with random noise vector

Let us start by briefly explaining the TRG method for the 2-D Ising model. The
partition function $Z$ is given as a trace over the tensor indices in a
tensor-network representation,

$\displaystyle Z=\sum_{\\{\sigma\\}}e^{-\beta
H}=\sum_{ijk\cdots}T_{ijkl}T_{mni\ell}\cdots={\rm Tr}\left[\otimes T\right],$
(1)

where $T_{ijkl}$ is the initial tensor. We express the tensor in a matrix
representation, $T_{ijkl}=M_{ij;kl}\equiv M_{ab}$, and
$T_{ijkl}=\tilde{M}_{jk;li}\equiv\tilde{M}_{cd}$. Using the SVD, we represent
a matrix $M$ with a rank $R$ as follows,

$\displaystyle
M_{ab}=\sum_{s=1}^{R}\sqrt{\Lambda_{s}}u_{as}\sqrt{\Lambda_{s}}v_{sb},$ (2)

where $\Lambda_{s}$ is a singular value
($\Lambda_{1}\geq\Lambda_{2}\geq\Lambda_{3}\cdots$). Then the fourth-order
tensor is approximated by a product of two low-rank third-order tensors with a
rank $D\leq R$ via the SVD,

$\displaystyle
T_{ijkl}=M_{ab}\simeq\sum_{s=1}^{D}\sqrt{\Lambda_{s}}u_{as}\sqrt{\Lambda_{s}}v_{sb}=\sum_{s=1}^{D}S_{3as}S_{1sb},$
(3)

where $S_{3as}=\sqrt{\Lambda_{s}}u_{as}$, and
$S_{1sb}=\sqrt{\Lambda_{s}}v_{sb}$ are the third-order tensors. Similarly,
$S_{2,4}$ are also defined to approximate the other type of a tensor
decomposition, $\tilde{M}_{cd}\simeq\sum_{s=1}^{D}S_{4cs}S_{2sd}$. We then
consider a coarse graining by contracting all the internal indices of the four
third-order tensors of $S_{1,2,3,4}$. By repeating these processes $n$-times
we obtain a renormalized tensor
$T^{(n)}_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ as

$\displaystyle
T^{(n)}_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}=\sum_{ijkl}S^{(n-1)}_{1i^{\prime};ij}S^{(n-1)}_{2j^{\prime};kj}S^{(n-1)}_{3lk;k^{\prime}}S^{(n-1)}_{4il;l^{\prime}}\equiv{\rm
Tr}[S^{(n-1)}_{1i^{\prime}}S^{(n-1)}_{2j^{\prime}}S^{(n-1)}_{3k^{\prime}}S^{(n-1)}_{4l^{\prime}}],\
\ (n>0)$ (4)

where $T^{(0)}=T$ and $S^{(0)}_{1,2,3,4}=S_{1,2,3,4}$. The bond dimension of
the renormalized tensor $D_{\rm cut}$ is determined by the number of the
singular values which also determines the accuracy of the method. The dominant
computational cost for contractions scales as $\mathcal{O}(D_{\rm cut}^{6})$,
(the randomized TRG $\mathcal{O}(D_{\rm cut}^{5}))$.

As shown above, TRG and other TRG-inspired renormalization methods use a low-
rank approximation technique to reduce the computational cost, while it
becomes a source of the systematic uncertainty. In order to reduce the
systematic uncertainty with a limited number of bond dimensions, we propose a
different rank truncation method using random noises. We introduce
$D$-dimensional noise vectors
$\eta_{r}=(\eta_{1r},\eta_{2r},\cdots,\eta_{Dr})^{T}$ for $r=1,\cdots,N_{r}$,
where each vector component $\eta_{ir}$ has a random number. We define a
$D\times N_{r}$ matrix $\bm{\eta}=(\eta_{1},\cdots,\eta_{N_{r}})$, which
satisfies the completeness condition

$\displaystyle\frac{1}{N_{r}}{\bm{\eta}}\cdot{\bm{\eta}}^{\dagger}=\sum_{r=1}^{N_{r}}\frac{1}{N_{r}}\eta_{r}\cdot\eta_{r}^{\dagger}={\bm{1}}_{D\times
D}+\mathcal{O}(1/\sqrt{N_{r}}).$ (5)

This random noise vector can be used to decompose a matrix in combing with the
SVD. Substituting
${\bm{1}}\sim\frac{1}{N_{r}}{\bm{\eta}}\cdot{\bm{\eta}}^{\dagger}$ into Eq.
(3) with $D=R-D_{\rm svd}$ we approximate the matrix $M$ via the SVD as

$\displaystyle M_{ab}$ $\displaystyle\simeq\sum_{s=1}^{D_{\rm
svd}}(\sqrt{\Lambda_{s}}u_{as}\sqrt{\Lambda_{s}}v_{sb})+\sum_{s,t=D_{\rm
svd}+1}^{R}\sqrt{\Lambda_{s}}u_{as}\left(\frac{1}{N_{r}}{\bm{\eta}}\cdot{\bm{\eta}}^{\dagger}\right)_{st}\sqrt{\Lambda_{t}}v_{tb}.$
(6)

From this decomposition we suggest a modification of the third-order tensor to
include the noise vector part. For this purpose, we introduce the following
modified third-order tensors as a function of the noise vectors,

$\displaystyle\bar{S}_{3as}({\bm{\eta}})$
$\displaystyle\equiv\begin{cases}\sqrt{\Lambda_{s}}u_{as}&(1\leq s\leq D_{\rm
svd})\\\ \sum_{i=D_{\rm
svd}+1}^{R}\sqrt{\frac{\Lambda_{i}}{N_{r}}}u_{ai}\eta_{is}&(D_{\rm svd}+1\leq
s\leq D_{\rm svd}+N_{r})\end{cases},$
$\displaystyle\bar{S}_{1sb}({\bm{\eta}})$
$\displaystyle\equiv\begin{cases}\sqrt{\Lambda_{s}}v_{sb}&(1\leq s\leq D_{\rm
svd})\\\ \sum_{i=D_{\rm
svd}+1}^{R}\sqrt{\frac{\Lambda_{i}}{N_{r}}}\eta_{si}^{*}v_{ib}&(D_{\rm
svd}+1\leq s\leq D_{\rm svd}+N_{r})\end{cases},$ (7)

and $\bar{S}_{2,4}({\bm{\eta}})$ are also defined similarly. Using the
modified tensors $\bar{S}_{1,2,3,4}$ with Eq. (5) we immediately see
$\displaystyle M_{ab}=\lim_{N_{r}\to\infty}\sum_{s=1}^{D_{\rm
svd}+N_{r}}\bar{S}_{3as}\bar{S}_{1sb}$ for any value of $D_{\rm svd}$.
Therefore, we can stochastically calculate the matrix decomposition. Since the
modified third-order tensors include all the modes for the singular values,
there is only a statistical error due to noise vectors instead of the
truncation error, which ensures that the original matrix should be reproduced
within the statistical error. It should be emphasized that a similar
decomposition has already been employed in lattice QCD calculations for
solving the Dirac equation, where a dominant part of the Dirac propagator is
obtained using a low-mode part of the Dirac eigenvectors, and a high-mode part
is stochastically evaluated using the random noise vectors. The graphical
representation of the tensor decomposition is given in Fig. 1.

Figure 1:  Tensor decomposition and modified third order tensors in the random
noise vector method.

##  3 Hybrid stochastic method

While the truncation error can be replaced with the statistical error, the
bond dimension increases as $N_{r}$ increases and the low-rank approximation
of the matrix decomposition is not as good as the original TRG in the sense of
the Frobenius norm due to the theorem by Eckhart-Young-Mirsky. Instead of
increasing $N_{r}$, we adopt a statistical approach. Namely, we generate an
ensemble of random noise vectors ${\bm{\eta}}^{[\ell]}\ (\ell=1,\cdots,N)$
with $N$ being the total number of statistics, and $\ell$ labels a sample
number. Our strategy is that a matrix decomposition is approximately obtained
from an ensemble average of the modified third-order tensors with keeping a
lower value of $N_{r}$. In fact, in this hybrid stochastic method we can
decompose the matrix as follows,

$\displaystyle\displaystyle
M_{ab}=\frac{1}{N}\sum_{\ell=1}^{N}\sum_{s=1}^{D_{\rm
svd}+N_{r}}\bar{S}_{3as}({\bm{\eta}}^{[\ell]})\bar{S}_{1sb}({\bm{\eta}}^{[\ell]})+\mathcal{O}(1/\sqrt{N_{r}N}).$
(8)

Thus, an exact matrix decomposition can be obtained in the infinite
statistical limit ($N\to\infty$) while keeping $N_{r}$ finite. This property
is used to effectively reduce the bond dimensions. In the next section, we
discuss how to obtain a renormalized tensor using the above tensor
decompositions.

###  3.1 Coarse graining with common noise

In order to preserve the lattice isotropy which is important to efficiently
carry out multiple steps of RG transformation, we notice that the same noise
vector should be used for obtaining $\bar{S}_{1,3}$ or $\bar{S}_{2,4}$. Thus,
a simple way of coarse graining is to use a common noise vector for each site.
A renormalized tensor is defined as

$\displaystyle
T^{(1)[\ell]}_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}({\bm{\eta}}_{1},{\bm{\eta}}_{2},{\bm{\eta}}_{1},{\bm{\eta}}_{2})={\rm
Tr}[\bar{S}_{1i^{\prime}}({\bm{\eta}}^{[\ell]}_{1})\bar{S}_{2j^{\prime}}({\bm{\eta}}^{[\ell]}_{2})\bar{S}_{3k^{\prime}}({\bm{\eta}}^{[\ell]}_{1})\bar{S}_{4l^{\prime}}({\bm{\eta}}^{[\ell]}_{2})],$
(9)

where the tensor indices $i^{\prime}j^{\prime}k^{\prime}$, and $l^{\prime}$
run from 1 to $D_{\rm svd}+N_{r}$, and hence the bond dimension $D_{\rm cut}$
for a renormalized tensor $T^{(1)}$ is given as $D_{\rm cut}=D_{\rm
svd}+N_{r}$. Since the same renormalized tensor is obtained at each site, it
is sufficient to calculate a tensor contraction at a single plaquette for each
RG step, while ${\bm{\eta}}^{[\ell]}_{1}$ and ${\bm{\eta}}^{[\ell]}_{2}$ are
different in general (See Fig. 2). By iterating the RG process, we obtain a
sequence of renormalized tensors for each sample $\ell$,

$\displaystyle T^{(0)}\to T^{(1)[\ell]}\to T^{(2)[\ell]}\cdots\to
T^{(n)[\ell]},\ \ (\ell=1,\cdots,N),$ (10)

in which the ensemble of noise vectors should be generated for each RG step
independently. Thus, any physical observables are calculated as an ensemble
average in terms of $T^{(n)[\ell]}$. We note, however, that there exists an
unwanted additional systematic error that arises from a cross contamination
effect due to a multiple use of the same noise vectors. We will discuss the
effect later, and in order to avoid the additional systematic error we will
consider position-dependent noise vectors.

Figure 2:  RG transformation of the tensor network with common noises.

As a benchmark study of the common noise method, we calculate the free energy
density $f$ in the 2-D Ising model. We employ the $Z_{2}$ noise, i.e.,
$\eta_{ir}$ takes $\pm 1$ randomly. We calculate an averaged $f$ with $N$
statistics at the $n$-th RG step as

$\displaystyle
f=\frac{1}{2^{n}}\frac{1}{N}\sum_{\ell=1}^{N}(-T\log{Z(T^{(n)[l]})}),\quad
Z(T^{(n)[l]})={\rm
Tr}\left[T^{(n)[l]}\right]=\sum_{i,j=1}^{D_{\rm{cut}}}T^{(n)[l]}_{ijij}.$ (11)

The statistical error is estimated from the variance of the data, where we
discard some unphysical data ($Z(T^{(n)[l]})<0$) from the analysis. In the
left panel of Fig. 3, we show the result of the relative deviation from the
Onsager’s exact result, $\delta f=(f-f_{\rm ex})/f_{\rm ex}$ at the critical
temperature $T=T_{\rm c}$ for several values of noise vectors. In comparison,
the results obtained from the original TRG ($D_{\rm svd}=50,N_{r}=0$) are also
shown. We can observe a much better accuracy than the original TRG (more than
$\mathcal{O}(10^{-1}$) magnitude difference) by optimizing the numbers of
singular values and noise vectors. We also show a scattering plot of $\delta
f$ for each sample in the right panel of Fig. 3. The fluctuation is well
controlled and we find that all the samples except the unphysical samples have
better accuracy than the original TRG. We also observe a plateau (RG fixed-
point tensors) for $n\gtrsim 28$ as in the case of the original TRG, in which
the statistical error does not grow as increasing the RG steps. This plateau
also indicates that there exists a systematic deviation that can not be
reduced by increasing the statistics. We consider this systematic deviation as
a cross contamination effect due to a multiple use of the same noise vectors.
Since we use a common noise vector for each RG process, e.g. in
$\bar{S}_{1,3}$ (Eq.(9)), we have a contact term originated from e.g.
$\frac{1}{N_{r}^{2}}\sum_{r,s=1}^{N_{r}}\eta_{ir}\eta_{rj}^{*}\eta_{ks}\eta_{sl}^{*}$,
which causes a systematic deviation from the Kronecker deltas. In fact, the
existence of a contact term is inevitable to preserve the isotropy, so that we
could not completely remove all the systematic error in the common noise
method. We note, however, that the resulting error is much smaller than the
original TRG. It should be also noted that this noise method can be
straightforwardly implemented in other tensor networks with the same order of
the computational cost per sample.

Figure 3:  Benchmark results for the RG flow in the common noise method.
(Left) Relative deviation between our result and the exact result for the free
energy density $|\delta f|$ with several values of the noise dimensions at the
critical temperature $T=T_{c}$. The number of statistics is $N=200$ for all
the results and the bond dimension is fixed to $D_{\rm cut}=50$. The TRG
result ($D_{\rm svd}=50,N_{r}=0$) is also shown for comparison. (Right) The
sample dependence of $\delta f$ at $n=50$ (RG step). The dashed (red) line
represents the mean value.

###  3.2 Coarse graining with position-dependent noise

In order to completely remove the systematic error, we shall consider
position-dependent noise vectors. The renormalized tensor is defined as

$\displaystyle T^{(1)[\ell]}_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}={\rm
Tr}[\bar{S}_{1i^{\prime}}({\bm{\eta}}^{[\ell]}_{1})\bar{S}_{2j^{\prime}}({\bm{\eta}}^{[\ell]}_{2})\bar{S}_{3k^{\prime}}({\bm{\eta}}^{[\ell]}_{3})\bar{S}_{4l^{\prime}}({\bm{\eta}}^{[\ell]}_{4})],$
(12)

where ${\bm{\eta}}^{[\ell]}_{i}$ is independently generated at each site $i$.
We refer this method as position-dependent noise method. Since the isotropy is
completely broken, we have to calculate a renormalized tensor for each
plaquette, where we also have to take into account the boundary conditions in
a finite volume. Therefore, the system volume $V$ has to be fixed, resulting
in a restriction on the maximum number of the RG transformations. See Fig. 4
for multiple steps of RG transformations in this method. At each coarse-
graining step, we obtain a statistical ensemble of tensor configurations, so
that this approach is similar to a Monte Carlo method, where the probability
of a spin variable is given by the Boltzmann distribution, while in this case
the tensor variables are generated randomly by noise vectors in analogous to a
random walk. Thanks to this tensor configurations the long distance
correlation can also be taken into account automatically. If we repeat this RG
process until being a single tensor, the total computational cost will be
proportional to the volume. Therefore, this method has the disadvantage of
computational cost scaling as $\mathcal{O}(VD^{6})$ per configuration. It,
however, guarantees that a correct result should be obtained within the
statistical error, as we desired.

Figure 4:  RG transformation of the tensor network with position-dependent
noise vectors for $2\times 2$ lattice. In the initial lattice four different
noise vectors ${\bm{\eta}}_{1,2,3,4}$ are distributed to the four sites. The
tensor contraction for each plaquette is carried out. After the RG
transformation, there exists a $1\times 2$ lattice, where two different noise
vectors ${\bm{\eta}}_{5,6}$ are distributed. The tensor contraction for one
plaquette is carried out, and a renormalized tensor $T^{(2)}$ is obtained.

Fig. 5 shows our preliminary results for the position-dependent noise method.
Since an ensemble average of the partition function
$Z^{(n)}=\frac{1}{N}\sum_{\ell=1}^{N}T^{(n)[\ell]}$ has no systematic error,
we calculate the free energy density from the partition function as
$f=-\frac{T}{2^{n}}\log{Z^{(n)}}$. The results give a much more accurate
result than both the common noise and the original TRG methods for each RG
step. But the results largely fluctuate depending on RG steps and the
statistical error may not be correctly estimated for larger RG steps due to a
small number of statistics $(N=20)$. Since the fluctuation comes from
truncated components of the SVD, a more stable and accurate result should be
obtained by increasing $D_{\rm svd}$. A detailed study will be reported in
Ref. [8].

Figure 5:  Preliminary results for the position-dependent noise method in
comparison with the common noise method and the original TRG for $N=20$.

##  4 Summary

We have proposed new low-rank approximation methods for TRG by utilizing
random noise vectors in combining with the SVD. In this method we generate
tensor ensembles, and any physical observables are statistically calculated
from this ensembles. In the case of the common noise method, there exists a
small systematic error due to a multiple use of the same noise vectors. In
order to completely remove the truncation error, we consider position-
dependent noise vectors, by which a statistical ensemble for the tensor
configuration is generated. Since there is no systematic error, it is
guaranteed that the result should become exact in the infinite statistical
limit with finite $D_{\rm cut}$. We have tested these two methods in the 2-D
Ising model, where we obtain a better accuracy in both cases than the original
TRG with a small statistics. The computational costs for the common noise
method per sample is the same order of the original TRG. On the other hand,
the position-dependent noise method will cost as expensive as
$\mathcal{O}(VD_{\rm cut}^{6})$. It should be emphasized that our noise
methods can be implemented to other coarse graining processes of tensor
networks as long as the SVD is used. A more detailed study of the noise
methods are given in Ref. [8].

## Acknowledgments

We would like to thank Shinji Takeda for fruitful discussions. M. T. is
supported by U.S. DOE grant #DE-SC0010339. H. O. is supported in part by JSPS
KAKENHI Grant Numbers 17K14309, 18H03710, and 21K03554.

## References

  * [1] M. Levin and C. P. Nave, Phys. Rev. Lett. 99, no.12, 120601 (2007) doi:10.1103/PhysRevLett.99.120601 [arXiv:cond-mat/0611687 [cond-mat.stat-mech]].
  * [2] Z.-C. Gu and X.-G. Wen, Phys. Rev. B 80, 155131 (2009).
  * [3] S. Yang, Z.C. Gu, X.G. Wen, Phys. Rev. Lett. 118, 110504 (2017)
  * [4] G. Evenbly, G. Vidal, Phys. Rev. Lett. 115, 180405 (2015).
  * [5] M. Hauru, C. Delcamp and S. Mizera, Phys. Rev. B 97, 045111 (2018).
  * [6] Satoshi Morita and Ryo Igarashi and Hui-Hai Zhao and Naoki Kawashima, Physical Review E. 97 033310 (2018).
  * [7] Y. Nakamura, H. Oba and S. Takeda, Phys. Rev. B 99, no.15, 155101 (2019) doi:10.1103/PhysRevB.99.155101 [arXiv:1809.08030 [cond-mat.stat-mech]].
  * [8] H. O. et al., in preparation.


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

Emily M. Boudreaux Department of Physics and Astronomy, Dartmouth College,
Hanover, NH 03755, USA  Elisabeth R. Newton Department of Physics and
Astronomy, Dartmouth College, Hanover, NH 03755, USA  Nicholas Mondrik
Department of Physics , Harvard University , 17 Oxford Street Cambridge MA
02138, USA  David Charbonneau Center for Astrophysics, Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, USA  Jonathan Irwin  Center for
Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA

(Received 23 September, 2021; Revised 1 March, 2022; Accepted 9 March, 2022)

###### Abstract

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^{\prime}_{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.

M dwarf stars (982), Stellar Activity (1580)

##  1 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 $\alpha\Omega$ 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 (Thompson et al., 1996), 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., 2004, 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
$\alpha$ and $\Omega$ 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 number111The Reynolds
Number is the ratio of magnetic induction to magnetic diffusion; consequently,
a plasma with a larger magnetic Reynolds number will sustain a magnetic field
for a longer time than a plasma with a smaller magnetic Reynolds number..

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 ($\tau_{c}$) — dominates
the motion of a charge packet within a star. Therefore, the Rossby Number
($Ro=P_{rot}/\tau_{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 transition between the saturated and non-saturated
regions occurs at $Ro\sim 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-to-late 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^{\prime}_{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.

##  2 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. (2016, 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
$\sim 5$ up to $\sim 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.

##  3 Analysis

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

$S=\alpha\frac{f_{H}+f_{K}}{f_{V}+f_{R}}$ (1)

where $f_{H}$ and $f_{K}$ are the integrated flux over triangular passbands
with a full width at half maximum of $1.09\text{ \AA}$ centered at
$3968.47\text{ \AA}$ and $3933.66\text{ \AA}$, 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, $\alpha$ is a scaling factor
with $\alpha=2.4$.

Figure 1: Spectrum of 2MASS J06105288-4324178 overplotted with the S index
bandpasses. (top) V band and Ca II K emission line. (bottom) Ca II H emission
line and R band. Note that the rectangular and triangular regions denote both
the wavelength range of the band and the relative weight assigned to each
wavelength in the band while integrating.

Following the procedure outlined in Lovis et al. (2011) we use the mean flux
per wavelength interval, $\tilde{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 $\lambda_{b,i}$ and a red most
wavelength $\lambda_{r,i}$, $\tilde{f}_{i}$ is the summation of the product of
flux ($f$) and weight ($w_{i})$ over the passband.

$\tilde{f}_{i}=\frac{\sum_{l=\lambda_{b,i}}^{\lambda_{r,i}}f(l)w_{i}(l)}{\lambda_{r,i}-\lambda_{b,i}}$
(2)

where $w_{i}$ represents the triangular passband for $f_{H}$ & $f_{K}$ and the
tophat for $f_{V}$ & $f_{R}$.

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 $\alpha_{H}$,
defined as:

$\alpha_{H}=8\alpha\left(\frac{1.09\text{ \AA}}{20\text{ \AA}}\right)$ (3)

Therefore, S-indices are calculated here not based on the historical
definition given in Equation 1; rather, the slightly modified version:

$S=\alpha_{H}\frac{\tilde{f}_{H}+\tilde{f}_{K}}{\tilde{f}_{V}+\tilde{f}_{R}}$
(4)

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^{\prime}_{HK}$ index (Middelkoop,
1982) is a transformation of the S-index intended to remove the contribution
of the photosphere.

$R^{\prime}_{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 mid-late 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^{\prime}_{HK}$ is defined as

$R^{\prime}_{HK}=K\sigma^{-1}10^{-14}C_{cf}(S-S_{phot})$ (5)

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\times 10^{6}\text{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.
$\sigma$ is the Stephan-Boltzmann constant. If we define

$R_{phot}\equiv K\sigma^{-1}10^{-14}C_{cf}S_{phot}$ (6)

then we may write $R^{\prime}_{HK}$ as

$R^{\prime}_{HK}=K\sigma^{-1}10^{-14}C_{cf}S-R_{phot}.$ (7)

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^{\prime}_{HK}$ as the standard deviation of
a distribution of $R^{\prime}_{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${}^{\prime}_{HK}$ uncertainty varies drastically with signal-to-noise;
targets with signal-to-noise ratios $\sim 5$ have typical uncertainties of a
few percent whereas targets with signal-to-noise ratios $\sim 100$ have
typical uncertainties of a few tenth of a percent.

###  3.1 Rotation and Rossby Number

The goal of this work is to constrain the rotation activity relation;
therefore, in addition to the measured $R^{\prime}_{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, $\tau_{c}$, such that the Rossby Number,
$Ro=P_{rot}/\tau_{c}$ . To first order $\tau_{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
$\tau_{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:

$\log_{10}(\tau_{c})=(0.64\pm 0.12)+(0.25\pm 0.08)(V-K)$ (8)

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.

##  4 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\pm 2.1$ km s-1
(Kesseli et al., 2018) and is therefore quite rotationally broadened, which is
known to affect $R^{\prime}_{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\geq 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).

$\log(R^{\prime}_{HK})=\begin{cases}\log(R_{s})&Ro<Ro_{s}\\\
k\log(Ro)+\log(R_{s})-k\log(Ro_{s})&Ro\geq Ro_{s}\end{cases}$ (9)

$Ro_{s}$ is the Rossby number cutoff between the saturated and unsaturated
regime. $R_{s}$ is the maximum, saturated, value of $R^{\prime}_{HK}$ and $k$
is the index of the power law when $Ro\geq Ro_{s}$. Due to the issues
measuring $R^{\prime}_{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.

2MASS ID | Mass | $Ro$ | $\log(R^{\prime}_{HK})$ | $\log(R^{\prime}_{HK})_{err}$ | $V_{mag}$ | $V-K$ | prot | $r_{prot}$ | Estimate
---|---|---|---|---|---|---|---|---|---
| $\mathrm{M_{\odot}}$ |  |  |  | $\mathrm{mag}$ | $\mathrm{mag}$ | $\mathrm{d}$ |  |
06000351+0242236 | 0.24 | 0.020 | -4.5475 | 0.0021 | 11.31 | 5.268 | 1.809 | 2016ApJ…821…93N | False
02125458+0000167 | 0.27 | 0.048 | -4.6345 | 0.0014 | 13.58 | 5.412 | 4.732 | 2016ApJ…821…93N | False
01124752+0154395 | 0.28 | 0.026 | -4.4729 | 0.0017 | 14.009 | 5.240 | 2.346 | 2016ApJ…821…93N | False
10252645+0512391 | 0.11 | 0.000 | -4.9707 | 0.0380 | 18.11 | 7.322 | 0.102 | 2016ApJ…821…93N | False
05015746-0656459 | 0.17 | 0.873 | -5.0049 | 0.0028 | 12.2 | 5.464 | 88.500 | 2012AcA….62…67K | False
06022261-2019447 | 0.23 | 1.307 | -5.6980 | 0.0192 | 13.26 | 4.886 | 95.000 | This Work | False
06105288-4324178 | 0.30 | 0.705 | -5.2507 | 0.0139 | 12.28 | 4.968 | 53.736 | 2018AJ….156..217N | False
09442373-7358382 | 0.24 | 0.542 | -5.6026 | 0.0147 | 15.17 | 5.795 | 66.447 | 2018AJ….156..217N | False
14211512-0107199 | 0.24 | 1.160 | -5.5846 | 0.0125 | 13.12 | 5.027 | 91.426 | 2018AJ….156..217N | False
14294291-6240465 | 0.12 | 0.394 | -5.0053 | 0.0014 | 11.13 | 6.746 | 83.500 | 1998AJ….116..429B | False
16352464-2718533 | 0.23 | 1.423 | -5.5959 | 0.0108 | 14.18 | 5.182 | 122.656 | 2018AJ….156..217N | False
16570570-0420559 | 0.24 | 0.014 | -4.3071 | 0.0014 | 12.25 | 5.130 | 1.212 | 2012AcA….62…67K | False
02004725-1021209 | 0.34 | 0.188 | -4.7907 | 0.0026 | 14.118 | 5.026 | 14.793 | 2018AJ….156..217N | False
18494929-2350101 | 0.18 | 0.034 | -4.5243 | 0.0015 | 10.5 | 5.130 | 2.869 | 2007AcA….57..149K | False
20035892-0807472 | 0.33 | 0.946 | -5.6530 | 0.0077 | 13.54 | 5.254 | 84.991 | 2018AJ….156..217N | False
21390081-2409280 | 0.21 | 1.152 | -6.1949 | 0.0190 | 13.45 | 5.091 | 94.254 | 2018AJ….156..217N | False
23071524-2307533 | 0.30 | 0.720 | -5.2780 | 0.0077 | 13.587 | 4.849 | 51.204 | 2018AJ….156..217N | False
00094508-4201396 | 0.30 | 0.009 | -4.3392 | 0.0018 | 13.62 | 5.397 | 0.859 | 2018AJ….156..217N | False
00310412-7201061 | 0.31 | 0.906 | -5.3879 | 0.0074 | 13.69 | 5.245 | 80.969 | 2018AJ….156..217N | False
01040695-6522272 | 0.17 | 0.006 | -4.4889 | 0.0024 | 13.98 | 5.448 | 0.624 | 2018AJ….156..217N | False
02014384-1017295 | 0.19 | 0.034 | -4.5400 | 0.0022 | 14.473 | 5.284 | 3.152 | 2018AJ….156..217N | False
03100305-2341308 | 0.40 | 0.028 | -4.2336 | 0.0017 | 13.502 | 4.935 | 2.083 | 2018AJ….156..217N | False
03205178-6351524 | 0.33 | 1.029 | -5.6288 | 0.0096 | 13.433 | 5.238 | 91.622 | 2018AJ….156..217N | False
07401183-4257406 | 0.15 | 0.002 | -4.3365 | 0.0022 | 13.81 | 6.042 | 0.307 | 2018AJ….156..217N | False
08184619-4806172 | 0.37 | 0.021 | -4.2834 | 0.0025 | 14.37 | 5.019 | 1.653 | 2018AJ….156..217N | False
08443891-4805218 | 0.20 | 1.348 | -5.6682 | 0.0067 | 13.932 | 5.370 | 129.513 | 2018AJ….156..217N | False
09342791-2643267 | 0.19 | 0.007 | -4.3415 | 0.0025 | 13.992 | 5.373 | 0.694 | 2018AJ….156..217N | False
09524176-1536137 | 0.26 | 1.342 | -5.6319 | 0.0110 | 13.43 | 4.923 | 99.662 | 2018AJ….156..217N | False
11075025-3421003 | 0.25 | 0.068 | -4.2250 | 0.0032 | 15.04 | 5.633 | 7.611 | 2018AJ….156..217N | False
11575352-2349007 | 0.39 | 0.031 | -4.2952 | 0.0026 | 14.77 | 5.415 | 3.067 | 2018AJ….156..217N | False
12102834-1310234 | 0.36 | 0.435 | -4.6892 | 0.0029 | 13.83 | 5.418 | 42.985 | 2018AJ….156..217N | False
12440075-1110302 | 0.18 | 0.020 | -4.4053 | 0.0033 | 14.22 | 5.546 | 2.099 | 2018AJ….156..217N | False
13442092-2618350 | 0.35 | 2.032 | -5.9634 | 0.0253 | 13.253 | 4.968 | 154.885 | 2018AJ….156..217N | False
14253413-1148515 | 0.51 | 0.301 | -4.7641 | 0.0030 | 13.512 | 5.121 | 25.012 | 2018AJ….156..217N | False
14340491-1824106 | 0.38 | 0.271 | -4.6093 | 0.0038 | 14.346 | 5.638 | 30.396 | 2018AJ….156..217N | False
15154371-0725208 | 0.38 | 0.050 | -4.6214 | 0.0023 | 12.93 | 5.224 | 4.379 | 2018AJ….156..217N | False
15290145-0612461 | 0.46 | 0.095 | -4.2015 | 0.0017 | 14.011 | 5.230 | 8.434 | 2018AJ….156..217N | False
16204186-2005139 | 0.45 | 0.031 | -4.3900 | 0.0035 | 13.68 | 5.261 | 2.814 | 2018AJ….156..217N | False
16475517-6509116 | 0.17 | 0.889 | -4.8744 | 0.0045 | 13.98 | 5.101 | 73.142 | 2018AJ….156..217N | False
20091824-0113377 | 0.15 | 0.010 | -4.3772 | 0.0023 | 14.47 | 5.958 | 1.374 | 2018AJ….156..217N | False
20273733-5452592 | 0.35 | 1.520 | -5.9982 | 0.0181 | 13.18 | 5.259 | 136.924 | 2018AJ….156..217N | False
20444800-1453208 | 0.49 | 0.073 | -4.4912 | 0.0023 | 14.445 | 5.305 | 6.715 | 2018AJ….156..217N | False
15404341-5101357 | 0.10 | 0.318 | -5.0062 | 0.0081 | 15.26 | 7.317 | 93.702 | 2018AJ….156..217N | False
22480446-2422075 | 0.20 | 0.005 | -4.4123 | 0.0016 | 12.59 | 5.384 | 0.466 | 2013AJ….146..154M | False
06393742-2101333 | 0.26 | 0.952 | -5.2524 | 0.0069 | 12.77 | 5.120 | 79.152 | 2018AJ….156..217N | False
04130560+1514520 | 0.30 | 0.019 | -4.4775 | 0.0088 | 15.881 | 5.437 | 1.881 | 2016ApJ…818..46M | False
02411510-0432177 | 0.20 | 0.004 | -4.4272 | 0.0016 | 13.79 | 5.544 | 0.400 | 2020ApJ…905..107M | False
11381671-7721484 | 0.12 | 0.958 | -5.5015 | 0.0369 | 14.78 | 6.259 | 153.506 | This Work | True
12384914-3822527 | 0.15 | 2.527 | -6.0690 | 0.0156 | 12.75 | 5.364 | 241.913 | This Work | True
13464102-5830117 | 0.48 | 1.340 | -5.6977 | 0.0146 |  |  | 65.017 | This Work | True
15165576-0037116 | 0.31 | 0.157 | -4.0704 | 0.0024 | 14.469 | 5.364 | 15.028 | This Work | True
19204795-4533283 | 0.18 | 1.706 | -5.8392 | 0.0091 | 12.25 | 5.405 | 167.225 | This Work | True
21362532-4401005 | 0.20 | 1.886 | -5.8978 | 0.0168 | 14.14 | 5.610 | 207.983 | This Work | True
Table 1: Calculated Rossby Numbers and $R^{\prime}_{HK}$ values. All circular
data points in Figures 2 & 3 are present in this table. Masses are taken from
the MEarth database. A machine readable version of this table is available
Figure 2: Rotation activity relation from this work. The color axis gives each
stars mass. The dashed line is the best fit to our data set. Figure 3:
Rotation activity relation for both our work and Astudillo-Defru et al.
(2017). The dotted line is the best fit to the re-derived rotation-activity
relation from Astudillo-Defru et al. (2017). Note that targets from Astudillo-
Defru et al. (2017) are systematically higher than targets presented here as a
consequence of the range in mass probed by the samples. Figure 4: Derived
rotation-activity curves from this work, Astudillo-Defru et al. (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.

We find best fit parameters with one $\sigma$ errors:

  * •

$k=-1.347\pm 0.203$

  * •

$Ro_{s}=0.155\pm 0.045$

  * •

$\log(R_{s})=-4.436\pm 0.048$

A comparison of the rotation activity derived in this work to those from both
Astudillo-Defru et al. (2017) and Mamajek & Hillenbrand (2008) is presented in
Figure 4. For the 6 targets which do not have measured rotational periods we
include an estimate of $Ro$ and $p_{rot}$ in the machine readable version of
Table 1. The convective overturn timescale for one of these 6 targets (2MASS
J13464102-5830117) can not be inferred via Equation 8 as it lacks a V-K color
measurement. Instead, we infer $\tau_{c}$ via Wright et al. (2018) Equation 6
(this paper Equation 10) using mass. Similar to our manner of inferring
$\tau_{c}$ via color, when inferring $\tau_{c}$ via mass, we adopt the larger
of the two antisymmetric errors from Wright et al. (2018).

$\log_{10}(\tau_{c})=2.33\pm 0.06-1.5\pm 0.21\left(M/M_{\odot}\right)+0.31\pm
0.17\left(M/M_{\odot}\right)^{2}$ (10)

Note that $R^{\prime}_{HK}$ for one of six of these targets (2MASS
J15165576-0037116) is consistent to within 1$\sigma$ 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^{\prime}_{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 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\pm 0.117$, consistent with
0 to within 1$\sigma$. Moreover, all other parameter for this model are
consistent to within one $\sigma$ 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\sigma$
confidence level. Ultimately, we adopt the most standard activity
interpretation, a fully-saturated 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, rotation-activity 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 regime — to 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 x-ray 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
($\tau_{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 $\tau_{c}$ (Wright et al., 2018). Doing this for all
targets presented in Astudillo-Defru et al. (2017) Table 3 and fitting the
same piecewise power law as before, we find best fit parameters of
$Ro_{s}=0.17\pm 0.04$, $\log(R_{s})=-4.140\pm 0.067$, and $k=-1.43\pm 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 H-alpha; 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_{\odot}$) and fit the same function as above, we find that all
best fit parameters are consistent to within one sigma between the two data-
sets.

Figure 5: Distribution of masses between our sample and the sample presented
in Astudillo-Defru et al. (2017). Note how the two studies have approximately
the same sample sizes; however, our sample is more tightly concentrated at
lower masses \later spectral classes.

We also compare our best fit $Ro_{s}$ to both those derived in Newton et al.
(2017) using $H_{\alpha}$ 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.

##  5 Conclusions

In this work we have approximately doubled the number of M-dwarfs with both
empirically measured $R^{\prime}_{HK}$ with $M_{*}<0.5M_{\odot}$. This has
enabled us to more precisely constrain the rotation-activity relation. This
relationship is consistent with other measurements using $R^{\prime}_{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.

## References

  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astudillo-Defru et al. (2017) Astudillo-Defru, N., Delfosse, X., Bonfils, X., et al. 2017, A&A, 600, A13
  * Barnes et al. (2005) Barnes, J., Cameron, A. C., Donati, J.-F., et al. 2005, Monthly Notices of the Royal Astronomical Society: Letters, 357, L1
  * Barnes et al. (2004) Barnes, J., James, D., & Cameron, A. C. 2004, Monthly Notices of the Royal Astronomical Society, 352, 589
  * Baroch et al. (2020) Baroch, D., Morales, J. C., Ribas, I., et al. 2020, A&A, 641, A69
  * Bernstein et al. (2003) Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Proc. SPIE, ed. M. Iye & A. F. M. Moorwood, 1694–1704
  * Berta et al. (2012) Berta, Z. K., Irwin, J., Charbonneau, D., Burke, C. J., & Falco, E. E. 2012, AJ, 144, 145
  * Browning et al. (2010) Browning, M. K., Basri, G., Marcy, G. W., West, A. A., & Zhang, J. 2010, AJ, 139, 504
  * Cameron et al. (2017) Cameron, R. H., Dikpati, M., & Brandenburg, A. 2017, Space Science Reviews, 210, 367.  https://doi.org/10.1007/s11214-015-0230-3
  * Catalano & Marilli (1983) Catalano, S., & Marilli, E. 1983, A&A, 121, 190
  * Chabrier & Baraffe (1997) Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039
  * Charbonneau (2014) Charbonneau, P. 2014, ARA&A, 52, 251
  * Delfosse et al. (1998) Delfosse, X., Forveille, T., Perrier, C., & Mayor, M. 1998, A&A, 331, 581
  * Dikpati & Charbonneau (1999) Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508
  * Earl et al. (2021) Earl, N., Tollerud, E., Jones, C., et al. 2021, astropy/specutils: v1.3, vv1.3, Zenodo, doi:10.5281/zenodo.4987675.  https://doi.org/10.5281/zenodo.4987675
  * Gastine et al. (2014) Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76
  * Hall et al. (2007) Hall, J. C., Lockwood, G. W., & Skiff, B. A. 2007, The Astronomical Journal, 133, 862.  https://doi.org/10.1086%2F510356
  * Hawley et al. (1996) Hawley, S. L., Gizis, J. E., & Reid, I. N. 1996, AJ, 112, 2799
  * Hawley & Pettersen (1991) Hawley, S. L., & Pettersen, B. R. 1991, ApJ, 378, 725
  * Houdebine et al. (2017) Houdebine, E. R., Mullan, D. J., Bercu, B., Paletou, F., & Gebran, M. 2017, ApJ, 837, 96
  * Irwin et al. (2011) Irwin, J., Berta, Z. K., Burke, C. J., et al. 2011, ApJ, 727, 56
  * Kelson (2003) Kelson, D. D. 2003, PASP, 115, 688
  * Kelson et al. (2000) Kelson, D. D., Illingworth, G. D., van Dokkum, P. G., & Franx, M. 2000, ApJ, 531, 159
  * Kesseli et al. (2018) Kesseli, A. Y., Muirhead, P. S., Mann, A. W., & Mace, G. 2018, AJ, 155, 225
  * Lehtinen et al. (2020) Lehtinen, J. J., Käpylä, M. J., Olspert, N., & Spada, F. 2020, arXiv e-prints, arXiv:2007.00040
  * Lovis et al. (2011) Lovis, C., Ségransan, D., Mayor, M., et al. 2011, Astronomy & Astrophysics, 528, A112
  * Magaudda et al. (2020) Magaudda, E., Stelzer, B., Covey, K. R., et al. 2020, A&A, 638, A20
  * Mamajek & Hillenbrand (2008) Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264
  * Mann et al. (2016) Mann, A. W., Gaidos, E., Mace, G. N., et al. 2016, ApJ, 818, 46
  * Medina et al. (2020) Medina, A. A., Winters, J. G., Irwin, J. M., & Charbonneau, D. 2020, ApJ, 905, 107
  * Middelkoop (1982) Middelkoop, F. 1982, A&A, 107, 31
  * Newton et al. (2017) Newton, E. R., Irwin, J., Charbonneau, D., et al. 2017, The Astrophysical Journal, 834, 85
  * Newton et al. (2016) Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93
  * Newton et al. (2018) Newton, E. R., Mondrik, N., Irwin, J., Winters, J. G., & Charbonneau, D. 2018, AJ, 156, 217
  * Noyes et al. (1984) Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763
  * Núñez et al. (2015) Núñez, A., Agüeros, M. A., Covey, K. R., et al. 2015, ApJ, 809, 161
  * Nutzman & Charbonneau (2008) Nutzman, P., & Charbonneau, D. 2008, PASP, 120, 317
  * Owen & Mohanty (2016) Owen, J. E., & Mohanty, S. 2016, MNRAS, 459, 4088
  * Pallavicini et al. (1981) Pallavicini, R., Golub, L., Rosner, R., et al. 1981, ApJ, 248, 279
  * Parker (1955) Parker, E. N. 1955, ApJ, 122, 293
  * Parker (1975) —. 1975, ApJ, 198, 205
  * Perdelwitz et al. (2021) Perdelwitz, V., Mittag, M., Tal-Or, L., et al. 2021, A&A, 652, A116
  * Pizzolato et al. (2000) Pizzolato, N., Maggio, A., Micela, G., et al. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 198, Stellar Clusters and Associations: Convection, Rotation, and Dynamos, ed. R. Pallavicini, G. Micela, & S. Sciortino, 71
  * Rauscher & Marcy (2006) Rauscher, E., & Marcy, G. W. 2006, PASP, 118, 617
  * Reiners et al. (2014) Reiners, A., Schüssler, M., & Passegger, V. M. 2014, ApJ, 794, 144
  * Robertson et al. (2014) Robertson, P., Mahadevan, S., Endl, M., & Roy, A. 2014, Science, 345, 440
  * Rutten (1984) Rutten, R. G. M. 1984, A&A, 130, 353
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55.  https://doi.org/10.7717/peerj-cs.55
  * Schmidt et al. (2014) Schmidt, S. J., Prieto, J. L., Stanek, K. Z., et al. 2014, ApJ, 781, L24
  * Schrijver & Rutten (1987) Schrijver, C. J., & Rutten, R. G. M. 1987, A&A, 177, 143
  * Schröder et al. (2009) Schröder, C., Reiners, A., & Schmitt, J. H. M. M. 2009, A&A, 493, 1099
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1
  * Skumanich (1972) Skumanich, A. 1972, ApJ, 171, 565
  * Suárez Mascareño et al. (2016) Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12
  * Suárez Mascareño et al. (2015) Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., & Esposito, M. 2015, MNRAS, 452, 2745
  * Thompson et al. (1996) Thompson, M. J., Toomre, J., Anderson, E. R., et al. 1996, Science, 272, 1300
  * Tomczyk et al. (1996) Tomczyk, S., Schou, J., & Thompson, M. J. 1996, in American Astronomical Society Meeting Abstracts, Vol. 188, American Astronomical Society Meeting Abstracts #188, 69.03
  * Torres et al. (2006) Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695
  * Vanderburg et al. (2016) Vanderburg, A., Plavchan, P., Johnson, J. A., et al. 2016, Monthly Notices of the Royal Astronomical Society, 459, 3565
  * Vaughan et al. (1981) Vaughan, A. H., Baliunas, S. L., Middelkoop, F., et al. 1981, ApJ, 250, 276
  * Vaughan & Preston (1980) Vaughan, A. H., & Preston, G. W. 1980, Publications of the Astronomical Society of the Pacific, 92, 385
  * Vilhu (1984) Vilhu, O. 1984, A&A, 133, 117
  * West et al. (2015) West, A. A., Weisenburger, K. L., Irwin, J., et al. 2015, ApJ, 812, 3
  * Wilson (1963) Wilson, O. C. 1963, ApJ, 138, 832
  * Wilson (1968) —. 1968, ApJ, 153, 221
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, The Astrophysical Journal, 743, 48.  https://doi.org/10.1088%2F0004-637x%2F743%2F1%2F48
  * Wright et al. (2018) Wright, N. J., Newton, E. R., Williams, P. K. G., Drake, J. J., & Yadav, R. K. 2018, MNRAS, 479, 2351