Measuring the Spot Variability of T Tauri Stars Using Near-IR Atomic Fe and Molecular OH Lines

0000-0003-4247-1401]Shih-Yun Tang
Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA
Lowell Observatory, 1400 West Mars Hill Road, Flagstaff, AZ 86001, USA

0000-0002-8828-6386]Christopher M. Johns-Krull
Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA

0000-0001-7998-226X]L. Prato
Lowell Observatory, 1400 West Mars Hill Road, Flagstaff, AZ 86001, USA

0000-0002-0848-6960]Asa G. Stahl
Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA

As part of the Young Exoplanets Spectroscopic Survey (YESS), this study explores the spot variability of 13 T Tauri Stars (TTSs) in the near-infrared $H$ band, using spectra from the Immersion GRating INfrared Spectrometer (IGRINS). By analyzing effective temperature () sensitive lines of atomic Fe i at $\sim$1.56259  and $\sim$1.56362 , and molecular OH at $\sim$1.56310 and $\sim$1.56317 , we develop an empirical equivalent width ratio (EWR) relationship for  in the range of 3400–5000 K. This relationship allows for precise relative  estimates to within tens of Kelvin and demonstrates compatibility with solar metallicity target models. However, discrepancies between observational data and model predictions limit the extension of the –EWR relationship to a broader parameter space. Our study reveals that both classical and weak-line TTSs can exhibit  variations exceeding 150 K over a span of two years. The detection of a quarter-phase delay between the EWR and radial velocity phase curves in TTSs indicates spot-driven signals. A phase delay of 0.06 $\pm$ 0.13 for CI Tau, however, suggests additional dynamics, potentially caused by planetary interaction, inferred from a posited 1:1 commensurability between the rotation period and orbital period. Moreover, a positive correlation between  variation amplitude and stellar inclination angle supports the existence of high-latitude spots on TTSs, further enriching our understanding of stellar surface activity in young stars.

§ INTRODUCTION

To comprehend how planets come into being, studying them while they are still forming is crucial. Yet searches for planets around young active stars, like T Tauri Stars (TTSs), often face challenges related to stellar activity <cit.>. The more we know about these surface features, the better chance we have of finding planets around young active stars.

One way for spots to hinder radial velocity (RV) planet searches is by generating periodic signals that mimic the RV variations induced by a planet. For simplicity, we refer to these spot-induced apparent RV shifts as “spot-RV signals". As star spots are cooler areas on the stellar disk, with local convective flux suppressed by strong magnetic fields, their nonaxisymmetric distribution will cause observed spectral lines to become asymmetric. Because spots co-rotate with the star, the periodic variation in the distortion of the absorption lines can be mistaken for a periodic RV signal, which can then be interpreted as a planet-induced signal <cit.>.

Purely data-driven methods for modeling the spot-RV signal, like Gaussian Process Regression (GPR) <cit.>, have lately seen wide use <cit.>; however, prior knowledge of the stellar rotation period and spot(s) lifetime is essential for the GPR to trace the activity signal properly. Knowing these stellar properties is even more critical for applying GPR to young stellar systems where the spot-RV signals from the extreme stellar activity can be several times larger than those induced by planets <cit.>. Young systems in which a planet's orbital period () is approximately equivalent to the rotation period of the host star () [Lanza, 2022] are even more challenging to characterize. An approximate 1:1 ratio between : has been found in mature solar-type systems like CoRoT-4 <cit.> and $\tau$ Boo <cit.>, and may not be surprising to find in young stars [Dawson & Johnson, 2018].

A star's  can often be determined from analyzing its lightcurve; however, for disk-bearing classical TTS (CTTS), the rotation period is often more difficult to determine compared with that of a weak-line TTS (WTTS) as the CTTS' photometry can be contaminated by the accretion variability, variable extinction, and disk-scattered light <cit.>. An alternative way to measure the star's  is to trace the apparent effective temperature () variation caused by changes in the area of cool spot coverage on the surface of a rotating star <cit.>. This provides a powerful, photometry-independent tool for the measurement of stellar .

Line-depth ratios (LDRs) have been an established technique for estimating the  of a stellar object since <cit.>. By estimating the LDR of a pair of isolated lines (e.g., V i at 6251.83 Å and Fe i at 6252.57 Å), this technique is immune to spectral broadening such as that caused by stellar rotation and instrumental resolution. With improvements over time, the precision of LDR-based relative  measurements for main-sequence and giant stars has reached $\leq$10 K <cit.>. This method has also facilitated the measurement of rotationally-modulated temperature variations for mature main-sequence stars such as $\epsilon$ Eri <cit.> and $\sigma$ Dra <cit.>, and more active stars like VY Ari <cit.>.

As TTSs have typical spectral types of late K to M, observing them in the near-IR is more efficient given that the spectral energy distributions of late-type stars peak in this wavelength region. Studies of the LDR– relationship in the near-IR began about two decades after work in the optical region. For example, Fukue et al., 2015 used absorption lines in the $H$ band to determine the temperature for G- and K-type giants and supergiants, Taniguchi et al., 2018 used lines in $Y$ and $J$ bands to establish the LDR- for G- and M-type giants, Kovtyukh et al., 2023 used near-IR lines for obtaining Cepheids' , and Afşar et al., 2023 used atomic lines in $H$ and $K$ bands to establish LDR– for late-type stars with a range of gravity () and metallicity ([M/H]). The effect of metallicity and  on the near-IR LDR– relationship has also been studied [Jian et al., 2019, Jian et al., 2020].

The  sensitive OH lines at $\sim$1.56310 [All wavelengths used in this study are in vacuum.] and 1.56317  were first used by O'Neal & Neff, 1997 and O'Neal et al., 2001 to study starspots on active stars. They reported a nearly linear relationship between the OH total equivalent width and the , shown an increasing trend from 5000 K to 3000 K. These OH lines, along with the nearby Fe i lines at $\sim$1.56259  and 1.56362  were later shown by Prato et al., 2002 to display inverse depth growth with spectral type from G0 to M9, i.e., the line depths increase for OH lines from G0 to M9, but the line depths of Fe lines get shallower from G0 to M9. López-Valdivia et al., 2019 went on to use LDs (as opposed to LDRs) of the aforementioned atomic and molecular lines to measure  for 162 K- and M-type dwarfs with the Immersion GRating INfrared Spectrometer <cit.>. These authors determined the temperature of each target, which they call $T_{\rm spec}$, by matching the targets' LDs to those from the BT-Settl spectral model grid [Allard et al., 2011, Allard et al., 2012]. They then calibrated systematic offsets between the models and observations using the color–temperature relation from Mann et al., 2015, ultimately establishing a linear $T_{\rm spec}$–LDR relationship from about 3300 to 3900 K. This relationship is more broadly applicable than that between $T_{\rm spec}$ and LDs, because LDRs are less affected by spectral broadening effects such as those induced by the instrumental profile. However, López-Valdivia et al., 2019 did not discuss how the $T_{\rm spec}$–LDR relationship might possibly change with stellar parameters such as , the line-of-sight projected equatorial rotation velocity (), metallicity, and average surface magnetic field strength ($\bar B$)  all of which have a larger range of values among TTSs than among main-sequence dwarfs [Sokal et al., 2020, López-Valdivia et al., 2021].

In this work, we explore the use of the OH/Fe LDR and equivalent width ratio (EWR) methods for assessing the relative apparent  variability of TTSs. We further apply these techniques to estimate the  of selected TTSs through their periodic LDR and EWR variations. The following sections lay out the structure of our study: Section <ref> first describes our target samples, the TTS and the  calibration sources, following with the IGRINS observations and data reduction. Section <ref> introduces the relevant spectral lines and analysis of the LD and EW ratios, taking into account the impact of spectral broadening. Section <ref> establishes empirical –EWR relationships from IGRINS data, which we show are more robust to use compared to the –LDR relationship. Section <ref> compares the observation results to spectral models. In Section <ref> we interpret the  variability observed in the TTS sample, present our  estimates, including a potential 1:1 : commensurability in the CI Tau system, and explore the dependence of  variations on other known properties of these stars. We conclude with a summary in Section <ref>.

§ TARGET SAMPLES, OBSERVATION, AND DATA REDUCTION

This study uses two target samples: the science target TTSs and the  calibration targets. The  calibration targets are main-sequence stars less active than TTSs, meaning almost no  variation. The  calibration targets also have known physical parameters, such as , allowing us to establish an empirical –EWR relationship after measuring the EWR from their spectra. This relationship is then used to estimate the  of TTSs after measuring their EWR from their spectra. In the following sections, we describe the two target samples and then discuss the IGRINS observation and data reduction.

§.§ T Tauri Star Sample

The TTS sample was selected from targets in the Young Exoplanets Spectroscopic Survey (YESS) [Prato et al., 2008, Crockett et al., 2012, Johns-Krull et al., 2016, Tang et al., 2023], a long-term program to search for stellar <cit.> and substellar <cit.> companions to young active stars using optical and near-IR RVs. Of the $\sim$120 TTSs in the YESS target list, we focussed on the 13 with at least six IGRINS observations taken during a single season. This cutoff maximizes our sample size while also ensuring sufficient data in each season for a substantive analysis of  variability.

The 13 targets can be broadly divided into two categories: those with two or more seasons of at least six observations in each season, and those with only one such season. The former category includes AA Tau, CI Tau, DK Tau, Hubble 4, V827 Tau, V830 Tau, and V1075 Tau; these targets' data not only allow the study of short-term (i.e., within a rotation period)  variations but also track long-term (seasonal) changes. The latter category includes DH Tau, DM Tau, DS Tau, GI Tau, IQ Tau, and LkCa 15.

The observations spanned five observational seasons, from 2014 to 2019, with a typical cadence of one observation per night for one or more weeks during each season. Figure <ref> displays the number of IGRINS observations for each target during these periods, while basic information on the targets can be found in Table <ref>. In Table <ref>, columns (2) to (11) contain values from existing literature, whereas column (12) present measurements from this study. The results of the analysis from the time series data for these TTSs are detailed in Table <ref>, where column (1) gives the target name, columns (2) and (3) list the observation UT dates and Julian Dates (JDs), respectively. The median signal-to-noise ratio (S/N) of the spectra from the reduction pipeline <cit.> is provided in column (4). For our TTS spectra the typical median S/N is around 200. Finally, measurements from this study are documented in columns (5) through (9).

The number of IGRINS observations of each TTS target over the five observing seasons between 2014–2019.
(a) Targets with two or more seasons having at least six observations; (b) Targets with only one season of at least six observations.

r C c RRRRRR c CCC
Basic Information for the TTS Sample

Name          Hmag$^a$                 SpT$^b$
$\bar B$      $\sigma \bar B$
$\sigma$     avg. RV

4-5 8-9 11-13
(1)     (2)     (3)
(4)     (5)
(6)     (7)
(8)     (9)
(10)     (11)     (12)

AATau     8.5     M0.6     3751     171     3.9     0.3     2.16     0.41        12.5     2.3     17.7

CITau     8.4     K5.5     3951     94     3.8     0.2     1.95     0.31        12.5     1.9     17.1

DHTau     8.8     M2.3     3477     125     3.9     0.2     2.21     0.32        8.4     2.1     15.9

DKTau     7.8     K8.5     3809     183     4.0     0.3     2.55     0.64        17.6     3.3     15.9

DMTau     9.8     M3.0     3449     104     4.1     0.2     1.63     0.29        5.7     2.0     18.5

DSTau     8.6     M0.4     3879     138     3.9     0.2     2.09     0.41        13.4     2.3     15.9

GITau     8.4     M0.4     3689     186     3.8     0.3     1.92     0.44        12.1     2.6     17.1

IQTau     8.4     M1.1     3612     177     3.7     0.3     1.81     0.53        14.3     2.5     15.5

LkCa15     8.6     K5.5     4156     123     4.1     0.2     1.83     9.99        15.4     2.3     17.9

HubbleI4     7.6     K8.0     3806     82     3.7     0.1     2.63     0.27        16.8     1.7     17.0

V1075Tau     9.1     K6.0     4122     97     4.2     0.2     2.62     9.99        30.9     2.1     18.2

V827Tau     8.5     M2.0     3610     92     3.9     0.1     2.42     9.99        20.8     1.8     18.0

V830Tau     8.6     K7.5     3878     80     3.9     0.1     2.4     9.99        31.6     1.8     17.0

$^a$ 2MASS H magnitude [Skrutskie et al., 2006].
$^b$ Spectra type from Luhman et al., 2017.
Information in columns (4)–(11) is from López-Valdivia et al., 2021.
Column (12) shows results from this study using  [Stahl et al., 2021, Tang et al., 2021, Tang et al., 2023]. The average radial velocity (avg. RV) values are calculated using data in Table <ref>.

§.§ Effective Temperature Calibration Sample

In order to derive an empirical –EWR relationship, we first identified targets for the calibration of effective temperature. These sources are nearby main-sequence late K and M dwarfs with IGRINS observations that have also been studied by <cit.>, <cit.>, and are included in the Gaia FGK benchmark star sample <cit.>. In the following, we first describe how  and  were calculated for each source and then introduce the [Fe/H] values adopted or derived in these studies.

§.§.§ Effective Temperature & Surface Gravity

The  values from [Boyajian et al., 2012] and GBS [Heiter et al., 2015] were estimated from the bolometric flux and the angular diameter of each stellar object using the Stefan–Boltzmann law . Stellar angular diameter can be measured directly using long-baseline interferometry (LBI), and bolometric flux can be calculated from combined photometry and/or spectra. This method therefore provides by far the most accurate  estimates <cit.>. A physical stellar radius can be calculated with parallax and the angular diameter of the star. The stellar mass can be estimated from the evolutionary tracks (GBS), or via the mass–luminosity relationship [Boyajian et al., 2012]. In this study, we use radii and masses from GBS and [Boyajian et al., 2012] to derive  of the target from the fundamental relation $g=GM/R^2$, where $M$ is the stellar mass, $R$ the stellar radius, and $G$ the gravitational constant. These  estimates are in agreement with literature values, when available. Passegger et al., 2022 reports  values for 10 of the targets in our final  calibration sample (Section <ref>) and in all cases the  adopted in this study is within the reported uncertainties around the literature median value.

The  values from [Mann et al., 2015] were estimated by fitting observed spectra to a model grid based on the method described in Mann et al., 2013. As this approach was tuned to match the results of LBI measurements to spectral fits with careful selection of the atmospheric models and wavelength regions used, the  measurements of [Mann et al., 2015] are comparable to those of GBS and [Boyajian et al., 2012]. The mean difference in  between [Mann et al., 2015] and literature LBI measurements is only $\sim$20$ \pm$11 K. Using the Stefan–Boltzmann law with parallax and bolometric flux as described above, [Mann et al., 2015] also estimated the angular diameter and physical radius for all their targets. The mean difference in angular diameter between their measurements and those derived from LBI is $\sim$1.4% $\pm$0.7%. Like [Boyajian et al., 2012], [Mann et al., 2015] also adopted a mass–luminosity relationship to estimate the masses and thus the  for each target. Figure <ref>a shows good agreement in  for the 13 targets common to both [Mann et al., 2015] and [Boyajian et al., 2012].

The accuracy of  measurements via interferometry critically hinges on the the measurement of angular diameters and on the calculation of bolometric fluxes, with the primary uncertainties in the latter stemming from the establishment of photometric zero points, model atmosphere selection, and extinction estimation. Interferometric calibration further introduces discrepancies, especially given systematic differences at between different facilities and beam combiners [Casagrande et al., 2014, White et al., 2018]. A comprehensive analysis by Tayar et al., 2022 identified systematic uncertainty floors of 2.4% $\pm$ 0.6% for bolometric fluxes and 4% $\pm$ 1% for angular diameters, contributing to an overarching  measurement uncertainty minimum of 2.0% $\pm$ 0.5%. For a typical  of TTS, 3800 K, the minimum uncertainty is about 75 K.

For targets with angular diameters $\lesssim$ 1 mas, Casagrande et al., 2014 identified a systematic offset where interferometric  measurements were found to be about 100 K higher compared to those derived using the infrared flux method (IRFM). To assess the impact of such discrepancies on our study, we performed IRFM  calculations for targets in [Boyajian et al., 2012] and [Mann et al., 2015] with available bolometric flux data. Utilizing the BT-Settl spectral library [Allard et al., 2011] for model flux and the IRTF Spectral Library [Cushing et al., 2005, Rayner et al., 2009] for estimating monochromatic flux on Earth, we narrowed our analysis to four targets with [Mann et al., 2015]'s bolometric fluxes (GJ 411, GJ 581, GJ 806, and GJ 846) and two targets with data from [Boyajian et al., 2012] (GJ 411 and GJ 581). A range of discrepancies in  are observed based on the bolometric flux source: $-81$ K to 87 K for [Mann et al., 2015]'s fluxes and $-119$ K to 132 K for [Boyajian et al., 2012]'s, with an average absolute discrepancy of 7 K across all measurements. Our analysis of this limited sample does not reveal a significant  offset. Instead, the differences seem largely tied to the bolometric flux calculation method.

ccCR RRRRR
Measurements of the TTS Samples

Name      UT              JD$-$2450000
RV        $\sigma$RV
EWR       $\sigma$EWR

(yyyy-mm-dd)    (days)

(1)     (2)     (3)
(4)     (5)
(6)     (7)
(8)     (9)

AATau     2015-01-05     7028.647     155     17.37     0.07     0.680     0.017     3839

AATau     2016-11-14     7707.980     105     17.63     0.08     0.695     0.020     3827

AATau     2016-11-17     7710.916     66     17.87     0.16     0.793     0.028     3756

AATau     2016-11-22     7715.891     68     18.77     0.08     0.756     0.027     3782

AATau     2016-11-23     7716.824     117     17.95     0.09     0.748     0.019     3788

(This table is available in its entirety in machine-readable form.)

Comparison of (a)  and (b) [Fe/H] for the  calibration sample targets common to both Mann et al., 2015 and Boyajian et al., 2012. The mean difference between the two studies for each parameter is shown; X symbols denote targets without a reported [Fe/H] uncertainty in Boyajian et al., 2012.

§.§.§ Metallicity

The metallicity calculation for the GBS [Jofré et al., 2014] was based on high-resolution and high S/N ratio optical spectra (i.e., HARPS, NARVAL, and UVES). The reported final [Fe/H] values were computed based on the abundances of selected Fe$\textsc{i}$ lines with an absolute solar Fe abundance from Grevesse et al., 2007. During the fitting process, the  and  were fixed to values determined from fundamental relations <cit.>. In [Mann et al., 2015], [Fe/H] values were calculated using the relationship between the near-IR atomic lines' EWs and [Fe/H], established in Mann et al., 2013, Mann et al., 2014, which made use of 156 binary systems with both a solar-type star of known [Fe/H] and a late K- or M-dwarf companion <cit.>. Lastly, [Boyajian et al., 2012] reported metallicities derived from a variety of sources: Anderson & Francis, 2011, Neves et al., 2012, Rojas-Ayala et al., 2012. Figure <ref>b compares [Fe/H] between the 13 targets common to both [Mann et al., 2015] and [Boyajian et al., 2012]. Unlike the  values shown in Figure <ref>a, the [Fe/H] display a systematic offset and scatter of $\sim$0.1, but overall the metallicities from [Boyajian et al., 2012] are consistent with those of [Mann et al., 2015] within the uncertainties. Also, not all [Fe/H] values collected by [Boyajian et al., 2012] have reported uncertainties.

r CCCCCCCCc RR cCC CC
Calibration Sources

Mass     $\sigma$Mass
R     $\sigma$R
[Fe/H]      $\sigma$[Fe/H]
Ref.2$^b$      S/N$^c$
EWR       $\sigma$EWR

(1)     (2)     (3)
(4)     (5)    (6)     (7)
(8)     (9)     (10)     (11)
(12)     (13)     (14)     (15)
(16)     (17)

GJ 581     3418     80     0.29     0.04     0.30     0.02     4.94     0.08     1,2     -0.15     0.08     1     535     0.00     1.419     0.021

GJ 725A     3424     61     0.33     0.05     0.35     0.01     4.85     0.07     1,2     -0.23     0.08     1     545     1.00     1.566     0.021

GJ 687     3426     66     0.41     0.06     0.42     0.02     4.81     0.07     1,2     0.05     0.08     1     597     0.00     1.470     0.019

GJ 436     3447     80     0.46     0.06     0.45     0.03     4.79     0.08     1,2     0.01     0.08     1     767     0.00     1.258     0.015

GJ 615.2C     3454     63     0.42     0.04     0.44     0.02     4.78     0.06     1     -0.06     0.03     1     236         1.256     0.024

GJ 70     3458     60     0.40     0.04     0.41     0.02     4.80     0.06     1     -0.13     0.08     1     265         1.404     0.025

GJ 625     3475     60     0.32     0.03     0.33     0.01     4.89     0.05     1     -0.35     0.08     1     374     0.00     1.653     0.034

GJ 745B     3494     62     0.30     0.03     0.32     0.01     4.90     0.06     1     -0.35     0.08     1     327     0.00     1.593     0.033

GJ 745A     3500     60     0.30     0.03     0.31     0.01     4.92     0.06     1     -0.33     0.08     1     345     0.00     1.584     0.030

GJ 3195     3500     61     0.41     0.04     0.42     0.02     4.80     0.06     1     -0.12     0.08     1     235         1.502     0.029

PM I18007+2933     3509     61     0.46     0.05     0.46     0.02     4.76     0.06     1     -0.06     0.03     1     218         1.289     0.026

GJ 411     3513     61     0.39     0.06     0.39     0.01     4.84     0.06     1,2     -0.38     0.08     1     978     0.00     1.640     0.018

GJ 806     3542     61     0.43     0.04     0.44     0.02     4.78     0.06     1     -0.15     0.08     1     552     0.00     1.279     0.018

GJ 393     3548     60     0.43     0.04     0.42     0.02     4.82     0.06     1     -0.18     0.08     1     466     0.00     1.335     0.020

GJ 412A     3558     71     0.40     0.06     0.39     0.02     4.85     0.07     1,2     -0.37     0.08     1     589     0.00     1.540     0.022

GJ 752A     3558     60     0.47     0.05     0.47     0.02     4.76     0.05     1     0.10     0.08     1     551     0.20     1.099     0.015

GJ 15A     3585     61     0.41     0.06     0.39     0.01     4.88     0.06     1,2     -0.30     0.08     1     810         1.446     0.019

GJ 382     3623     60     0.53     0.05     0.52     0.02     4.72     0.05     1     0.13     0.08     1     414     1.50     0.897     0.015

GJ 526     3633     67     0.49     0.07     0.48     0.02     4.76     0.07     1,2     -0.31     0.08     1     514     0.00     1.489     0.022

GJ 87     3638     62     0.44     0.04     0.44     0.02     4.79     0.06     1     -0.36     0.08     1     434         1.648     0.028

GJ 908     3646     60     0.41     0.04     0.41     0.01     4.83     0.05     1     -0.45     0.08     1     435     0.00     1.913     0.035

GJ 3408B     3656     62     0.43     0.04     0.42     0.02     4.83     0.06     1     -0.26     0.03     1     374         1.369     0.025

GJ 686     3657     60     0.44     0.04     0.42     0.01     4.83     0.05     1     -0.25     0.08     1     539     0.00     1.343     0.020

GJ 176     3680     60     0.49     0.05     0.45     0.02     4.82     0.06     1     0.14     0.08     1     227     0.00     1.048     0.019

GJ 887     3682     92     0.51     0.07     0.47     0.02     4.80     0.08     1,2     -0.06     0.08     1     348         1.294     0.021

PM I02441+4913W     3685     60     0.52     0.05     0.50     0.02     4.76     0.05     1     0.06     0.03     1     353         0.958     0.016

GJ 134     3700     61     0.64     0.06     0.63     0.03     4.65     0.06     1     0.53     0.08     1     442     1.40     0.603     0.013

GJ 649     3700     60     0.53     0.05     0.51     0.02     4.75     0.05     1     0.03     0.08     1     329     0.30     0.905     0.016

GJ 505B     3709     60     0.54     0.05     0.54     0.02     4.71     0.05     1     -0.12     0.03     1     396         0.938     0.015

GJ 880     3716     61     0.57     0.08     0.55     0.02     4.71     0.06     1,2     0.21     0.08     1     993     1.30     0.751     0.012

GJ 514     3727     61     0.53     0.05     0.48     0.02     4.79     0.05     1     -0.09     0.08     1     468     0.00     0.998     0.015

GJ 809     3741     63     0.58     0.08     0.54     0.02     4.73     0.07     1,2     -0.06     0.08     1     804     0.00     0.785     0.013

GJ 212     3765     60     0.59     0.06     0.57     0.02     4.70     0.06     1     0.19     0.03     1     327     1.90     0.695     0.014

GJ 281     3771     60     0.63     0.06     0.63     0.03     4.64     0.06     1     0.12     0.08     1     468     1.70     0.523     0.013

GJ 96     3785     62     0.61     0.06     0.60     0.02     4.67     0.05     1     0.14     0.08     1     920         0.606     0.012

GJ 205     3801     60     0.62     0.09     0.58     0.02     4.71     0.06     1,2     0.49     0.08     1     612     2.10     0.571     0.013

GJ 685     3846     61     0.59     0.06     0.54     0.02     4.74     0.05     1     0.10     0.08     1     225     0.90     0.742     0.015

GJ 846     3848     60     0.59     0.06     0.55     0.02     4.73     0.05     1     0.02     0.08     1     396     1.10     0.676     0.013

GJ 338B     3867     37     0.60     0.06     0.57     0.01     4.71     0.05     2     -0.15     0.17     2,4     549     0.50     0.579     0.013

GJ 338A     3913     69     0.61     0.09     0.56     0.03     4.72     0.08     1,2     -0.01     0.08     1     617     1.00     0.559     0.013

GJ 208     3966     60     0.65     0.07     0.60     0.02     4.69     0.05     1     0.05     0.08     1     337     2.00     0.519     0.013

GJ 820B     4044     32     0.64     0.06     0.60     0.02     4.67     0.04     1,2,3     -0.22     0.08     1     552         0.542     0.013

GJ 169     4124     62     0.74     0.07     0.69     0.02     4.63     0.05     1     0.39     0.08     1     1085     1.30     0.352     0.012

GJ 673     4124     60     0.71     0.07     0.65     0.02     4.66     0.05     1     0.19     0.08     1     550         0.374     0.013

GJ 820A     4374     22     0.68     0.07     0.66     0.00     4.63     0.04     2,3     -0.33     0.38     3     779         0.328     0.012

GJ 702B     4393     149     0.70     0.07     0.67     0.01     4.63     0.05     2     0.06     0.04     2,5,6     624         0.237     0.013

GJ 570A     4507     58     0.74     0.07     0.74     0.02     4.57     0.05     2     0.05     0.04     2,6,7     698         0.172     0.013

GJ 892     4699     16     0.76     0.08     0.78     0.01     4.54     0.04     2     0.01     0.04     2,5,6     913         0.134     0.013

eps Vir     4983     61                     2.77     0.02     3     0.15     0.16     3     418         0.183     0.014

$^a$ Reference for  and  measurements:
(1) Mann et al., 2015, (2) Boyajian et al., 2012, and
(3) Gaia FGK benchmark stars <cit.>.
$^b$ Reference for [Fe/H] measurement: (1)–(3) the same as reference
for  and . (4) Rojas-Ayala et al., 2012, (5) Ramírez et al., 2012,
(6) Luck & Heiter, 2006, and (7) Ghezzi et al., 2010.
$^c$ S/N refers to the median S/N values between 1.561–1.565  from the
IGRINS pipeline package version 2.2.0 <cit.>.
Column (15) give  measurements from Reiners et al., 2022, and columns
(16) and (17) give equivalent width ratio (EWR) measurements and their uncertainties
from this study.
(This table is available in machine-readable form.)

§.§.§ Final Samples

Targets with  between 3200 and 5000 K in [Boyajian et al., 2012], [Mann et al., 2015], and GBS were cross-matched with The Raw & Reduced IGRINS Spectral Archive <cit.>[<https://igrinscontact.github.io>]. With a S/N cut of $>$100, we are left with 49 final  calibration targets with high-quality IGRINS spectra, of which 19 have measurements from [Boyajian et al., 2012], 43 have measurements from [Mann et al., 2015], and three have measurements from GBS.

For targets with multiple measurements in  and , we took the mean and adopted an uncertainty based on the quadratic sum of all uncertainties. For [Fe/H], we adopted measurements from [Mann et al., 2015] and GBS, if available. If not, we used the values in [Boyajian et al., 2012], except for the three targets that lack [Fe/H] uncertainties in [Boyajian et al., 2012]. For these targets (GJ 702B, GJ 570A, and GJ 892), we draw [Fe/H] values and uncertainties from Luck & Heiter, 2006, Ramírez et al., 2012, and Ghezzi et al., 2010.

Table <ref> presents the final parameters for the  calibration targets. Column (1) shows the target name, and (2) and (3) provides the  values and their uncertainties, respectively. Columns (4) to (9) give stellar mass, stellar radii ($R$), , and their respective uncertainties. The sources for literature values of mass and radii are cited in column (10). Columns (11) and (12) provide the [Fe/H] values and uncertainties, with their references listed in column (13). The median S/N values from the plp reduction pipeline are included in column (14). Column (15) shows  values from Reiners et al., 2022. Measurements conducted in this study are in columns (16) and (17). Figure <ref> displays the IGRINS $H$ band spectra centered on the Fe and OH lines, arranged from top to bottom by increasing . The values of , [Fe/H], and  for the targets, as listed in Table <ref>, are also illustrated in the figure.

IGRINS spectra for the  calibration sample. Spectra are sorted by , with high  on the top and low  at bottom. The , [Fe/H], and  values are from Table <ref>.

§.§ Observations and Data Reduction

All spectra used in this study were taken with IGRINS, a cross-dispersed echelle spectrograph which can simultaneously cover the entire $H$ (1.49–1.80 ) and $K$ bands (1.96–2.46 ) while delivering a resolution of $R\sim$45,000 with a fixed width slit (08). With no moving parts, IGRINS has been installed at a number of different sites, including the McDonald Observatory's 2.7 m Harlan J. Smith Telescope, the 4.3 m Lowell Discovery Telescope, and the 8-m Gemini South telescope. Our target spectra were observed at these three sites from July 2014 to April 2019. Observations were all taken with one or more A-B nodding sequence(s), and the reduction was carried out using the IGRINS pipeline package version 2.2.0 <cit.>.

We downloaded the reduced spectra with telluric line correction from RRISA. The telluric line correction was performed by dividing the science spectrum by an A0V spectrum observed shortly before or after the science target. (typically within an air-mass difference of 0.1). As RRISA also provides a cross-match (XMatch) catalog with SIMBAD that has been checked by hand to ensure the XMatch accuracy, we used the RVs given in the XMatch catalog for the  calibration sample. For the TTS sample, we calculated RVs for each epoch using [<https://github.com/shihyuntang/igrins_rv>] [Stahl et al., 2021, Tang et al., 2021, Tang et al., 2023] because RV variations in these young active stars can be as high as several  <cit.>. The average RV for each target TTS, determined from , is given in column (12) of Table <ref>.

Spectra showing the $H$ band region centered around the Fe and OH lines, all normalized as per Section <ref>.
The three sets of colored lines on the top are model spectra (see Section <ref>) of various  with  = 4.0, and $R$ = 45,000; black lines at the bottom are the IGRINS spectra of all  calibration sources.
The locations of several strong lines (e.g., CN i, Fe i, OH i, Y i, and Ni i) are marked and a legend is shown for these in the upper left corner. Relatively strong, unlabelled lines that appear in sources with  $\lesssim$ 3600 K are mostly water lines. Vertical black lines mark the wavelength regions for the EW calculation. The two main Fe i lines discussed in Section <ref>, LFe i at $\sim$1.56259  and RFe i at $\sim$1.56362  are highlighted. Line information is taken from the VALD3 line list [Ryabchikova et al., 2015].

§ SPECTRAL ANALYSIS

§.§ The OH and Fe Region

The OH lines at $\sim$1.56310  and $\sim$1.56317  and the Fe i lines at $\sim$1.56259  (the left-hand side Fe i, LFe i) and $\sim$1.56362  (the right-hand side Fe i, RFe i, Figure <ref>) in the near-IR $H$ band are advantageous for LD and EW analysis for a number of reasons. The lines are close in wavelength, so they will be minimally affected by imperfect spectral normalization, and, for CTTSs, the wavelength dependent veiling. Additionally, the OH and Fe i lines have contrasting responses to changes in : LD increases with  for the Fe i lines from $\sim$3200 K to $\sim$5000 K, but for the OH lines, LD decreases with increasing  from $\sim$3400 K to $\sim$5000 K (Fig. <ref>).

The telluric contamination in this wavelength region is also minimal. The main source of telluric absorption is water vapor, which is only present at substantial levels during winter and at lower altitude observing sites (in this case, the McDonald Observatory and the Lowell Discovery Telescope). Telluric water absorption only affects the LFe i line and the shorter wavelength of the OH lines ($\sim$1.56310 ). The typical flux difference in this wavelength region, when choosing an A0V standard star with an air-mass difference of 0.1, is less than 0.3%, which is about a change of 0.005 in EWR (see Section <ref>).

By contrast, the effects of line blending are both significant and unavoidable given the large  values of our TTS sample (up to about $30$ , Table <ref>). To show the level of significance of blending, we display model spectra with various values of  and  in colored lines in Figure <ref>. The line blending effect can be seen, for example, in the LFe i line. It is blended with another Fe i line and a CN i line. Another example is the blending of the RFe i with two other Fe i lines, a CN i line, and a Ni i line. Some lines can be distinguished in the low  spectra but start to blend in with each other at a higher ; for example, the two separate Fe i lines within 1.5635–1.5636  are distinguishable when =0 but blend together at =10 .

§.§ Equivalent Width and Line Depth Calculation

To ensure LD and EW measurements are consistent between synthetic models and IGRINS spectra, as well as between stars with different , we develop a standardized normalization approach: all spectra are normalized to the 95th percentile value of the flux measured between 1.561–1.565  (Fig. <ref>). This wavelength region was chosen to ensure that spectra maintain a similar continuum level at $\sim$1.5612 and $\sim$1.5646  regardless of the rotational broadening exhibit by our targets, even with  as high as 30 .

To streamline the measurement process and minimize the blending effect discussed in Section <ref>, we calculated LD and EW using wavelength regions which include lines that might be blended at higher . The EW is measured from 1.56230–1.56290  for the LFe i line, between 1.56290–1.56335  for the OH lines, and from 1.56335–1.56395  for the RFe i line (these regions are demarcated by vertical black lines in Figure <ref>). The final LD and EW for Fe is the mean of the values measured in the two Fe i regions. The OH lines form a close doublet that is effectively unresolved at the IGRINS resolution of 45,000.

We used the  package  [Earl et al., 2022] to measure the EW of lines in the three wavelength regions. Uncertainties were calculated using the Monte Carlo (MC) method with 5000 spectra generated based on the measured S/N. We tested the robustness of this procedure by measuring the EWR in spectra of a known quiescent star, GJ 281, which has been used as a standard because of its low RV variability <cit.>. Our EWR measurements of GJ281 show a scatter of $\sigma_{\rm EWR}\sim$ 0.012 over 61 nights of observation spanning a period of 4 years (Fig. <ref> a). This scatter is about 2.8 times larger than the median uncertainty of individual EWR values calculated based on the S/N from the plp pipeline. To check for systematic errors in the EWR, we assumed GJ 281 has negligible intrinsic variation between observations and we used the average scatter in the residual between adjacent spectra, sorted by S/N, to infer a more realistic S/N. We find these residuals to be systematically larger than implied by the plp S/N estimate, and this yields scaling factors for the plp spectral uncertainty: 0.75 for plp S/N below 225, and 0.58 for those above. The revised S/N values lead to a median uncertainty of 0.0051, which is still significantly less (by a factor of 2.4) than the observed $\sigma_{\rm EWR}$ of GJ 281. Consequently, we define an additional systematic uncertainty, $\sigma_{\rm sys} = \sqrt{0.012^2 - 0.0051^2} \sim 0.011$, which is incorporated quadratically into all final EWR measurement uncertainties. Possible sources of this systematic uncertainty include inaccurate continuum normalization and the influence of spot coverage <cit.>. The EWR uncertainty introduced by adopting different A0V standards during the telluric correction (Section <ref>), could be one of these sources. Using the –EWR relationship established for solar metallicity targets (Section <ref>, Figure <ref>), we convert GJ 281's EWRs to , as shown in Figure <ref>b. The resultant  scatter, $\sigma_{T_{\rm eff}}\sim$ 12 K, indicates our measurement precision for  around 3950 K (an EWR of approximately 0.53). This analysis employs results from Section <ref>, discussed later, to maintain the logical flow of our findings and discussions.

The (a) EWR and (b)  time series for GJ 281 spanning about 4 years. The  values in (b) are calculated from the EWRs in (a) using the “solar” relationship given in Figure <ref>. The small scatter in EWR, $\sigma_{\rm EWR}\sim$ 0.012, and in the , $\sigma_{T_{\rm eff}}\sim$ 12 K, demonstrate the precision of our technique at $\sim$3950 K (EWR $\sim$ 0.53).

There are three common techniques to measure the LD: fit a Gaussian function to the main profile of the line <cit.>, fit a parabolic curve to pixels around the lowest flux <cit.>, and calculate the difference between the average flux near the minimum flux point and the continuum <cit.>. We tested all three methods in the three wavelength regions used for EW measurement (Fig. <ref>) by (a) fitting a Gaussian function within each region, (b) fitting a second-order polynomial to the 7 pixels ($\sim$2 resolution elements) around the minimum flux in each region, and (c) calculating the mean flux of the 3 pixels around the minimum flux in each region.

Because the spectra in our TTS sample exhibit generally lower S/N than the  calibration sample, we smoothed the spectra using a Gaussian (standard normal) kernel with $\sigma =$ 1 pixel before performing the LD measurements. This smoothing process helps minimize the effect of outlier flux measurements in the lower S/N spectra.

§.§ Spectral Broadening Effect

In order to be robust and comprehensive, an empirical –EWR (and/or LDR) relationship must not significantly depend on other factors such as , spectral resolution ($R$), or the surface average magnetic field strength ($\bar B$). We study the effect of these parameters on LDRs and EWRs using synthetic spectra with  = 4000 K,  = 4.0, and S/N = 150. Figure <ref> illustrates the impact on EWR and LDR of (a) , (b) $R$, and (c) magnetic field strength. Each LDR measurement method (Gaussian fit, polynomial fit, minimum) is displayed in a different color; note the high resolution of IGRINS spectra make the minimum (yellow) and polynomial (blue) results often indistinguishable.

Simulated EWR and LDR as a function of (a) , (b) spectral resolution ($R$), and (c) magnetic field strength.
Model spectra are generated with  = 4000 K,  = 4.0, [Fe/H] = 0.0, and S/N = 150. A  of 10  is adopted for model spectra used in (b) and (c), and an $R$ of 45,000 is adopted for model spectra used in (a) and (c).
Red lines show LDR calculated by fitting a Gaussian, orange lines show LDR calculated by adopting the mean flux of 3 pixels around the minimum flux, blue lines are LDR calculated by fitting a second-order polynomial to the 7 pixels around the absorption lines' minimum flux, and the black lines show the EWR (see Section <ref>).
In all panels, the EWR lines are shifted by the amount given in the legend for a better comparison. Also, in all panels, the orange and the blue lines are on top of each others. The EWR method is less affected by all three spectral broadening effects.

The –EWR plot for the  calibration sample.
(a) sample color coded by [Fe/H], and (b) sample color coded by . The empirical relationship of –EWR is fitted using a natural log function. The solid black line shows the fitting result for solar-type ($-0.15<$[Fe/H]$<0.15$) stars, highlighted with plotted outlined with black circles, and the dashed line shows the fitting result of all stars in the sample. The best-fit results and their associated uncertainties are given in panel (a). In (b), an outlier with a low  of 2.77 is marked with a star symbol.

First, we observe that LDRs clearly vary with  regardless of measurement method (Fig. <ref>a) because of the blending of lines. For example, the dip in LDR around a  of 4–10  with the minimum and polynomial fit methods is primarily the result of blending of the OH doublet. The maximum variation in the LDR from  of 0 to 30  is about 0.21, which is equivalent to a temperature difference of about 230 K at 4000 K if we adopt the solar –EWR relationship derived in the next section (Section <ref>, and Figure <ref>). The EWR measurements, on the other hand, only yield a maximum variation of about 0.02, equivalent to a $\sim$20 K difference at 4000 K.

To account for variations in the $R$ across the IGRINS detector <cit.>, and to generalize this method for use with other spectrographs, Figure <ref>b illustrates the relationship of LDR and EWR with $R$, ranging from 10,000 to 120,000. Notably, there is a marked decrease in LDR for $R$ values below approximately 30,000. This decline is attributed to the increasingly dominant influence of spectral resolution blending close lines, which becomes similarly significant to a  of 10 . In general, the impact of $R$ on the ratio is relatively minor (with a maximum LDR variation of 0.07) in comparison to the effect of .

The $B$ strength has the largest impact on both the LDR and EWR measurements. Given a typical range of $\bar B$ values from $\sim$1–3 kG for TTSs <cit.>, the effect on the ratio measurements are $\lesssim$ 0.04 for EWR and $\lesssim$ 0.13 for LDR. A difference of 0.04 in EWR corresponds to an approximate $\Delta$ of $\sim$15 K at an EWR of 1.5 (equivalent to $\sim$3400 K) and $\sim$72 K at an EWR of 0.3 (equivalent to $\sim$4300 K), according to the solar –EWR relationship. Unlike  and $R$, the response of different atomic lines to the magnetic field strength is unique to that particular line. Hence, the findings presented here are specific to the spectral region examined in this study. Ultimately, we adopt EWR, not LDR, for deriving the  relationship as it demonstrates reduced sensitivity to spectral broadening effects.

§.§ –EWR Empirical Relationship

Figure <ref> shows EWR versus  for our calibration sources. In Figure <ref>a, the sample is color-coded by [Fe/H], and in <ref>b, it is color-coded by . The tight –EWR trend at high  is weakened at temperatures cooler than $\sim$3700 K because of the presence of sub-solar metallicities and higher  targets (we later show in Section <ref> and Figure <ref> that the scatter in low temperature is mainly caused by metallicity differences). At EWR of $\sim$1.5, the  difference between a solar metallicity target and a sub-solar ([Fe/H]$\sim-0.4$) metallicity target can be as large as 200K. We therefore fit two empirical –EWR relationships: one including the entire sample, and the other using only a “solar” metallicity ($-0.15<$ [Fe/H] $<0.15$) sub-sample. Both –EWR relationships were fit with a natural log function ( $=a + b\ln{\rm EWR}$) using the orthogonal distance regression ( package ) to consider uncertainties in both  and EWR. Best-fit coefficients are given in Figure <ref>a. These empirical –EWR relationships are best characterized for EWR from about 0.1 ( $\sim 5000$ K), before it asymptotes, to about 1.5, where we start to lose constraints from the observed samples ( $\sim 3400$ K).

Comparison plots of Fe and OH equivalent widths (EW) from model spectra versus IGRINS observational data for the  calibration targets. The top row (panels a, b, and c) show the Fe EW outcomes, while the bottom row (panels d, e, and f) details the OH EW results. Column one illustrates the EW derived from the PHOENIX ACES spectral grid, column two from synmast using the VALD3 line list without water lines, and column three incorporates the same line list with additional BT2 water lines [Barber et al., 2006] in synmast. Data points are color-coded according to  as indicated by the color bar on the right. The dashed line across each panel represents the 1:1 correspondence line, serving as a reference for perfect agreement between modeled and observed EW values.

§ MODEL COMPARISON

To obtain the most precise  estimates possible, it would be ideal to customize the –EWR relationship for different , [Fe/H], and magnetic field strengths using models calibrated with observations.
We compare our observational results against theoretical models to explore the possibility that spectral models might provide a useful extension of our observed sample in this way.

First, we compare our empirically-derived relationship with measurements based off of the publicly available pre-computed high-resolution model spectral library PHOENIX-ACES-AGSS-COND-2011 <cit.>[<http://phoenix.astro.physik.uni-goettingen.de>]. We interpolate across the spectral grid to get model spectra for each of our  calibration stars with the adopted parameters from Table <ref>. The adoption of the $\alpha$ element (O, Ne, Mg, Si, S, Ar, Ca, and Ti) enhanced spectra was based on the following rule:
[$\alpha$/Fe] = 0.4 for $-5.0$ $\leq$ [M/H] $\leq-1.0$;
[$\alpha$/Fe] = $-0.4 \times$ [M/H] for $-1.0 \leq$ [M/H] $\leq$ 0.0; and [$\alpha$/Fe] = 0.0 for [M/H] $\geq$ 0.0 [Gustafsson et al., 2008]. The same procedures introduced in Section <ref> were then applied to get the model EW for PHOENIX-ACES spec. The first column of Figure <ref> compares the EW results derived from the PHOENIX-ACES spectra and IGRINS data for Fe (Fig. <ref>a) and OH (Fig. <ref>d). Deviation from the slope $=$ 1 line is evident, especially for Fe. This mismatch indicates that pre-computed spectral models are not a useful extension of observations when it comes to calibrating a –EWR relationship. We therefore explore a number of means to improve on these spectral models, relying on spectral synthesis to generate our own model spectra for which we can specify the spectral line list and make other detailed changes.

We utilize the synmast spectral synthesis code [Piskunov, 1999, Kochukhov, 2007, Kochukhov et al., 2010] because it allows us to incorporate the effects of $B$ fields. Given that TTSs can exhibit surface magnetic field strengths several kG higher than those found in M dwarfs, accounting for the influence of $B$ fields on the –EWR relationship could be crucial. To this end, we employ a line list from the VALD3 database [Ryabchikova et al., 2015][<http://vald.astro.uu.se>]. Additionally, a comparison between the results of using different model atmosphere profiles, i.e., BT-NextGen [Allard et al., 2011, Allard et al., 2012], PHOENIX-ACES <cit.>, and MARCS [Gustafsson et al., 2008] is shown in Appendix <ref>. Because the difference is minimal, we opt only to show results using the BT-NextGen atmospheres with ASGG2009 [Asplund et al., 2009] solar abundance.

Spectral tuning results using synmast, incorporating a tuned VALD3 line list, BT2 water lines, and adopted the BT-NextGen atmospheric profile. The tuning exercise involved a simultaneous fitting process applied to solar spectra (top row) and spectra from three M-dwarfs (bottom row): GJ 205 (first column), GJ 880 (second column), and GJ 338B (third column). Each panel displays the observed spectrum as a solid black line with the model spectra in color. The model spectra generated with the un-tuned line list are depicted by faded color lines, and the model spectra generated with tuned line list are in dashed color lines. Residuals between the observations and models with tuned line list are indicated by grey dashed lines, offset 0.55 for visibility. The $\chi^2$ statistic of fit is also shown. Key spectral lines, including CN i, Fe i, OH i, Y i, and Ni i, are identified and labeled within each panel. The six lines that underwent specific tuning are highlighted by inverted triangles.

Similar to Figure <ref> but showing synmast model results with VALD3 tuned line list, BT2 water lines, and adopted the BT-NextGen atmosphere profile. The simultaneous fit to solar and the three M-dwarfs are shown in the first column for GJ 205, the second for GJ 880, and the third for GJ 338B. Each of the M-dwarfs used in the tuning process is highlighted in each panel. The top row, panels (a), (b), and (c), show Fe results, and the bottom row, panels (d), (e), and (f) show OH results.

The middle column of Figure <ref> displays these results versus our data-derived  relationships for Fe EW, Figure <ref>b, and OH EW, Figure <ref>e. The synmast spectra produce a better match than the PHOENIX-ACES spec for Fe EWs, displaying a strong correlation with the data-derived results albeit with a slight offset. The synmast spectra also show a marked improvement for OH EWs among stars warmer than $\sim$3800 K. The more divergent results from synmast spectra cooler than $\sim$3800 K may be the result of the lack of water lines in VALD3 at these temperatures, and/or the relative insensitivity of the OH line depth at lower temperatures in synmast.

To further improve on these results, we experimented with adding a water line list to our synmast model generation. We tested two different water ($^1{\rm H}_2$$^{16}{\rm O}$) line lists, the BT2 [Barber et al., 2006] and the POKAZATEL [Polyansky et al., 2018] from the ExoMol database  [Tennyson & Yurchenko, 2012][<https://www.exomol.com>]. The difference between the resulting EWR by adding BT2 versus the POKAZATEL water lines is small, as shown in Appendix <ref> Figure <ref>; thus, here we only show results using the BT2 line list. The improvement in the model EWR after adding water lines is significant, as shown in the third column (panels c and f) of Figure <ref>. The extra opacity from water lines in cooler stars results in larger Fe and OH EWR values, making them move upward on the plot closer to the slope = 1 line. Nevertheless, a larger deviation from the unity line remains for Fe EWR in Figure <ref>c around targets with $\sim$ 4000 K. One reason for this deviation is that spectral lines from atomic and molecular databases are often not precise enough to match with observational data <cit.>. Below, we show our results on manually adjusting the log$gf$ and the van der Waals constant for six dominant lines.

§.§ VALD3 Line Tuning

It is well known that spectral line data collected in various databases are often not accurate enough to properly reproduce stellar spectra <cit.>. Investigators often choose to tune line parameters such as the oscillator strength (in this case log$gf$, where $g$ is the statistical weight and $f$ is the oscillator strength) and the Van der Waals (VdW) damping parameter by fitting the spectrum of the quiet Sun. In our case, however, simply adjusting the log$gf$ values and VdW damping parameters of the two dominant Fe i lines based on the NSO disk center solar atlas [Wallace et al., 1996] results in systematically weakened Fe i lines for stars with  $\leq$ 4000 K. This suggests a degree of degeneracy between the log$gf$ values and VdW damping parameters that cannot be resolved by fitting the solar spectrum alone. A simultaneous fit incorporating a cooler star is necessary to disentangle these effects. We opt to calibrate our line parameters through a joint fit of both solar and cool ( $\sim$ 3800 K) star spectra, given our focus on stars substantially cooler than the Sun, which exhibit distinct features like OH lines absent in the solar spectrum (Figure <ref>). The cool stellar spectra used were observed with IGRINS under similar S/N conditions, ensuring that variations in instrument, resolution, and S/N do not skew the tuning process.

We first created a line list using the VALD3 line database [Ryabchikova et al., 2015] by performing two “Extract Stellar" queries, one at  = 5775 K and one at  = 3800 K, and combining the two returned lists. Spectra generated using the original VALD3 line list are shown in Figure <ref> as faded color lines. The line tuning process is done by simultaneously fitting the Solar spectrum (here we are using the spectrum of the asteroid Ceres from IGRINS) and a M-dwarf spectrum in the  calibration sample. When performing this tuning, it is key to use stars with accurately known stellar parameters given that parameters such as , gravity, and metallicity affect line strength. While the Solar values are quite reliable, there is more uncertainty in these parameters for M stars. We therefore choose 3 stars from Table 2 with interferometrically determined  values and with metallicities from Mann et al., 2015, based on relations calibrated to FGK stars where metallicity determinations are more reliable. The log$g$ values are again based on the interferometric radii combined with a mass-luminosity relationship [Boyajian et al., 2012]. Passegger et al., 2022 collected log$g$ values from the literature for all of these stars and they show generally excellent agreement, giving confidence in the values we adopt.  Finally, we also select stars with minimal contamination from water lines since the line data for water has significant uncertainties. Using the spectra of GJ 205, GJ 880, and GJ 338B, we performed the line tuning three times using the Sun and each of these stars in turn. The locations of these M-dwarfs in EW– space are highlighted in Figure <ref>c and Figure <ref>f. We used the parameters in Table <ref> for the M stars when doing the tuning. For the Solar spectrum, we adopt  = 5775 K,  = 4.44, and [Fe/H] = 0.0.

We performed the line tuning using a non-linear least squares fitting procedure based on the Marquardt method [Bevington & Robinson, 1992] to simultaneously fit synthesized spectra to the observed spectrum of the Sun (Ceres) and the cooler star of choice. While the observed S/N for each spectrum varies somewhat, they are all above 100. We arbitrarily set the S/N for each observation to 100 in order to equally weight the fitting to each spectrum. The free parameters of the fit are the log$gf$ and VdW broadening terms of lines that showed significant mismatch between the initially synthesized spectra and the observations (Figure <ref>). For weaker lines, only the log$gf$ value is fit. Model spectra were calculated with the synthesis code synmast without water lines included, and the BT-NextGen models were interpolated to the specific stellar parameters (, , [Fe/H]) of each star using the PySME code [Wehrhahn et al., 2023].

§.§ Line tuning Results

Figure <ref> shows the line tuning results on the three M-dwarfs, and Figure <ref> shows the resulting models EW compared to the observation for the entire sample. Line tuning using the sun and GJ 205 seems to give the lowest residual between the model and observed spectra as shown in the first column of Figure <ref> (panels a and d). However, this set of tuned lines gives the most significant deviation of Fe EW between the models and the observations of the whole  calibration target set (Fig. <ref>a). GJ 338B, on the other hand, shows an inverse result. While lines tuned with the sun and GJ 338B give the largest residual between the model and the observed spectra for these stars (Fig. <ref>c and Fig. <ref>f), Fe EW comparison between the models and the observations have the tightest relationship for cool stars  $\lesssim$ 4400 K (Fig. <ref>c). The OH EWs (Fig. <ref>f) show a relatively low-scatter correspondence. Nevertheless, this tuning systematically underestimates the EW of the model for the hotter stars. Of the three M-dwarfs, the result for GJ 880 lies between GJ 205's and GJ 338B's results, thus this is the one we choose to use for the further analysis.

Line tuning effectively brings the model EW of selected targets closer to the equality line, as evidenced by the adjustments for GJ 205 and GJ 880 in Figure <ref>a and <ref>b, respectively. Nonetheless, this method does not address deviations for the entire sample. Additionally, even targets selected for tuning, such as GJ 338B, may still deviate from the equality line. This issue arises during simultaneous fitting with the Sun, which is presumed to have accurately determined physical parameters, limiting the enhancement of the VdW damping strength. The spectral fitting results for the Sun, when simultaneously fitting with GJ 338B, illustrated in Figure <ref>c, already shows excessively deep Fe cores, thus preventing the fitting routine from further increasing the strength of the VdW value.

As mentioned above, adopting an accurate set of stellar parameters is essential in the line tuning. However, this is challenging for the cool M-dwarfs, especially with respect to the [Fe/H] values. Passegger et al., 2022 studied metallicities of 18 M-dwarfs determined with different techniques, and all three of the M-dwarfs we used in the tuning are included in their sample. Passegger et al., 2022 collected literature [Fe/H] values for their entire sample and listed them in their Table A.1, where they also provided literature median values and associated uncertainties. Although the [Fe/H] values adopted for the three targets in our tuning process are all within the uncertainties in the literature median from Passegger et al., 2022, the range in the literature [Fe/H] values can be as large as 0.69. For example, for GJ 205, the lowest [Fe/H] is $+0.0 \pm 0.09$ [Maldonado et al., 2015] and the highest [Fe/H] is $+0.69 \pm 0.10$ [Terrien et al., 2015].

We evaluate our line tuning process on GJ 205 with three different [Fe/H] values (0.20, 0.39, and 0.49) from the range of estimates in Passegger et al., 2022. Using 0.49 from [Mann et al., 2015] as our initial (default) value yielded the smallest discrepancy between the model and observed spectra when fitting the Sun and GJ 205 simultaneously (Fig. <ref>). Using the other two [Fe/H] values did not produce simultaneous fits as good as that when we assume [Fe/H] = 0.49. Line tuning critically relies on accurate stellar parameters when making adjustments to database-provided line parameters to enhance the observational fit. However, if a target's stellar parameters are poorly defined, biases can arise, limiting line tuning to merely compensating for these uncertainties. As a result, while line tuning may align the EW of a specifically tuned target with the line of equality (Fig. <ref>), it cannot rectify overall sample deviations. Our experiments highlight the need for precisely defined physical parameters in calibration targets for  and improvements in spectral modeling for cooler M-dwarfs.

The effective temperature () versus equivalent width ratio (EWR) for model grids and IGRINS observational data. The model grids were produced using synmast with a VALD3 line list calibrated on solar and GJ 880 spectra, incorporating BT2 water lines within the BT-NextGen atmospheric profile. Panel (a) shows the effect of metallicity [Fe/H], (b) the influence of surface gravity (), and (c) magnetic field strength impact. Models are shown as thin lines, observational data as circles, and a grey line denotes the empirical fit for solar metallicity ($-0.15<$[Fe/H]$<0.15$), detailed in Section <ref> and Figure <ref>. Panels (a) and (b) assume 0 $B$ field strength, with (a) set at =4.5 and (b) using solar metallicity. Panel (c) is also at =4.5 and solar metallicity.

§.§ EWR Comparison of Model and Observation

Figure <ref> is similar to Figure <ref> showing  versus EWR, but now overplotted with model results generated using atomic lines tuned on GJ 880, the BT2 water list, and BT-NextGen atmospheres. By employing models, we are able to investigate the EW across a continuous range of , offering deeper understanding of the –EWR relationship and its dependency on parameters such as , [Fe/H], and the $B$ field. This approach overcomes the constraints of our limited sample.

In Figure <ref>a, we see that [Fe/H] has a more significant effect on the OH/Fe EWR toward cooler  compared to the impact of  in Figure <ref>b and $B$ field in Figure <ref>c. This larger difference is expected as our EWR includes Fe lines. The more intriguing result occurs at sub-solar metallicities below $\sim$3700 K where the models converge toward a single EWR value. The backtracking in sub-solar metallicity models below $\sim$3700 K occurs as Fe lines effectively vanish, leaving the models to reflect predominantly water opacities. For , the effect is indistinguishable between  of 3.5 to 5.0 above $\sim$3900 K. Models only start to be slightly separated below $\sim$3900 K to $\sim$3200 K. The small separation between different  models implies the  estimated using the –EWR relationship is less susceptible to the accuracy of  values. Lastly, we can see that there is almost no difference in models as a function of B field strength (Fig. <ref>c).

In Figure <ref>a and <ref>b we can also see that models with  = 4.5 and solar metallicity (solid black thin lines) match well with the empirical –EWR fit for solar metallicity targets (data points with open black circles, also see Section <ref>). Still, below $\sim$3400 K, the model deviates from the empirical fit, turning almost straight down to trace what would be a non-unique –EWR relationship. Taken altogether, these results indicate a robust empirical –EWR relationship for cool stars and TTSs in particular, so long as  is not below 3400 K, in agreement with our observations in Section <ref>.

§ DISCUSSION

§.§ Effective Temperature Variability in T Tauri Stars

To study the  variability of the 13 TTSs in our sample, we first measure their time series EWR using the procedures described in Section <ref>. We then use the empirical –EWR relationship in section <ref> to obtain their time series . As TTSs mostly have solar metallicities, e.g., D'Orazi et al., 2011, we use the “solar” empirical –EWR relationship. With time series  measurements, we study the overall  variability, likely caused by the variation in spot(s) coverage, of the TTSs (nine CTTSs and four WTTSs), then explore the seasonal  variations of the seven TTSs (three CTTSs and four WTTSs) that were most densely observed.

§.§.§ Overall Variability

The box plots in Figure <ref> display the  distribution of each of our TTS targets, with the extent of temperature variability illustrated by the size of the boxes, from the first quartile (Q1, 25th percentile) to the third quartile (Q3, 75th percentile). The two whiskers (vertical bars) correspond to the minimum and maximum values of the , such that the range the whiskers span gives the total amplitude of the  variation. Targets exhibiting significant  variability (amplitude $\gtrsim$150 K), include V1075 Tau ($\sim$150 K), DK Tau ($\sim$180 K), CI Tau ($\sim$240 K), and LkCa 15 ($\sim$210 K). These variations are all more than 3.5 times the measured  sensitivity ($\Delta$) of each target listed in Table <ref>. Targets with smaller  amplitudes ($\lesssim$60 K) include Hubble 4 ($\sim$65 K), DS Tau ($\sim$80 K), and DM Tau ($\sim$47 K). For these sources, the  variation is at least 2.5 times greater than the respective $\Delta$.

c CCCCC
TTS EWR Temperature Sensitivity

Name      avg. EWR
avg. $\sigma$EWR      $\Delta$

(1)     (2)     (3)
(4)     (5)

AA Tau     0.728     3804     0.020     30

CI Tau     0.512     3994     0.016     34

DH Tau     0.859     3714     0.018     23

DK Tau     0.659     3858     0.019     31

DM Tau     0.988     3637     0.016     18

DS Tau     0.604     3903     0.017     30

GI Tau     0.732     3801     0.019     28

Hubble 4     0.637     3875     0.014     24

IQ Tau     0.755     3783     0.018     26

LkCa 15     0.370     4170     0.016     47

V1075 Tau     0.383     4151     0.014     40

V827 Tau     0.763     3777     0.014     20

V830 Tau     0.544     3960     0.014     28

Column 2 displays each TTS's average EWR from all observations, while column 3 lists their corresponding  (K), estimated via the solar –EWR relationship shown in Figure <ref>. Column 4 lists the average uncertainties in EWR. By applying the average EWR from column 2 and its uncertainties from column 4 to the –EWR relationship, we calculate the systematic scatter in , the $\Delta$ (K), presented in column 5. This scatter is interpreted as the sensitivity of the  measurement for each TTS.

Box plot showing the  “variations” for the TTS sample, with CTTSs in red and WTTSs in blue. The boxes span the range between the first quartile (Q1, 25th percentile) and the third quartile (Q3, 75th percentile); the median is indicated with a dark colored horizontal line. The whiskers (upper and lower bars) show the maximum and minimum values. Black squares show  values estimated by López-Valdivia et al., 2021 using MCMC fitting to IGRINS spectra, and the errorbars show the “uncertainties” in the measurements.

No apparent trend links  variability to whether a star is classified as a CTTS or a WTTS. This may suggest similar spot distributions, sizes, and spot-to-photosphere temperature contrasts between both classes of stars. Moreover, as shown in Figure <ref>, where targets are ordered by their increasing median  for CTTSs and WTTSs separately, a positive correlation between  and the activity level, i.e., the  variation amplitude, is observed. To quantify this result, we calculate the Pearson correlation coefficient and Spearman’s rank correlation coefficient. We use the inter-quartile range (IQR, Q3$-$Q1) instead of the  variation amplitude to mitigate the effects of outliers. Both tests yield a correlation of 0.6 with p-values of 0.02, indicating a significant linear relationship. This result might suggest higher surface activity in TTSs with higher average . However, further observations and model testing are needed.

López-Valdivia et al., 2021 analyzed the same TTSs as shown in Figure <ref>, using the identical IGRINS dataset. However, they employed a forward modeling approach to determine  for each TTS, averaging between 1 and 10 observational epochs per target, with uncertainties assessed via MCMC fitting. This sets the stage for an insightful comparison between  “uncertainty” derived from forward modeling and the “variability” we identify through EWR measurements. In Figure <ref>, we show López-Valdivia et al., 2021's model fitting results (represented as black squares with error bars) against our findings. Notably, the  “uncertainty” reported by López-Valdivia et al., 2021 often matches or surpasses the “variability” detected in our study. Additionally, López-Valdivia et al., 2021's average  values are consistently lower than our median values derived through the EWR method. This systematic lower  is underscored by a significant negative correlation, indicated by a Pearson coefficient of $r=-0.75$ and a p-value of 0.003, suggesting a larger discrepancy at lower  values. This trend could be explained by the fact that López-Valdivia et al., 2021's  estimates are based off of measurements in the $K$ band instead of the $H$ band. At these longer wavelengths, the  estimates would be more sensitive to cool stellar spots [Gully-Santiago et al., 2017]: as the temperature of the photosphere decreases, the contrast between the spot and the photosphere increases more in the $K$ band than in the $H$ band. This would produce systematically lower temperature results, potentially explaining the  negative correlation observed.

Box plots showing the seasonal  “variation” of seven TTSs which have at least six observations in two or more seasons. Three observing seasons from 2016–2019 are shown. Red boxes are CTTSs, and blue boxes are WTTSs. The black bar in the lower left of each panel shows the  sensitivity from Table <ref>, column (5). The smaller the black bar's length compared to the variability, the more robustly the changes are characterized. The numbers at the top of each box show the total observations for each season.

§.§.§ Seasonal Variability

In this section, we delve deeper into the seasonal variability of seven TTSs that feature two or more seasons with at least six observations each. Figure <ref> illustrates the  distributions of these TTSs across three observation seasons (2016–2019). By focusing on the median values and the size of the IQR to minimize the impact of outliers, seasonal  variations become apparent for targets like AA Tau, CI Tau, V827 Tau, and V830 Tau. Interpreting these fluctuations is challenging without model simulations, as both spot size and temperature contrast with the stellar photosphere likely influence the apparent . However, assuming minor temperature contrast changes between seasons, an increase in overall , such as observed for V827 Tau, may suggest a reduction in spot sizes. Additionally, the IQR's “size” serves as an indicator of changes in a star's activity level across seasons. An increase in either the number or size of spots, along with greater temperature contrast with the photosphere, suggests heightened stellar activity, reflected by a larger IQR. Consequently, an expanding IQR, as seen with AA Tau, V830 Tau, and the last two seasons in V827 Tau, could signal rising stellar activity.

In addition to illuminating the overall  variability amplitudes, examining trends in EWR when folded to the target's  can provide a more insightful view as now we can examined the shape of the EWR curve. Figure <ref> displays, for each target, a Lomb-Scargle periodogram <cit.> analysis of RVs (red), LDRs (blue), and EWRs (black) in an upper panel. The lower panels show the EWRs and RVs phase folded at . We leave the discussion of the periodograms to the next section (Section <ref>), and only focus on the shapes of the folded EWR and RV curves here.

Figure 14 shows that a strong periodic signal is present for all targets [Hubble 4 is excluded from Figure <ref> because of the complexities arising from its binary system nature [Rizzuto et al., 2020]. Despite detecting a periodic signal when phase-folding the LDR and EWR to the reported  of 1.5459 days [Carvalho et al., 2021], a significant phase shift was observed, consistent with findings from Carvalho et al., 2021's optical RV data. Additionally, our observations were conducted during the periastron phase of the binary orbit. Without well-defined orbital parameters, detrending the long-term RV variations attributable to the binary system is challenging. V827 Tau is also reported as a binary system, but we include it in our analysis because only one prominent signal is present in its EWR and RV curves, which we expect to be solely from the primary star. This expectation is supported by a primary-to-secondary flux ratio of about 1.7 in the $H$ band, indicating that we mainly observe the primary star [Kraus et al., 2011].]. Moreover, color-coding each observation by time often reveals changes in the EWR and RV variations between seasons, i.e., on an annual time scale. For example, the EWRs of V827 Tau peak more sharply during the winter of 2016 (blue), around phase 1.25, than in the winters of 2017 (cyan) and 2018 (red). Its RV data, meanwhile, appears to trace a higher amplitude trend during the 2017 (cyan) season than in either 2016 (blue) or 2018 (red). The characteristics of the EWR and RV curves could offer insights into the size and configurations of stellar spots through model simulations. However, accurately interpreting these features is challenging without such modeling. It is plausible that spots appearing as patches could cause significant scatter in the EWR curve, whereas a larger, more uniform spot might result in a more concentrated EWR curve, as observed in V1075 Tau. Specifically, V1075 Tau exhibited a wide scatter in EWRs during the winter of 2017 (depicted in cyan), in contrast to a more defined pattern in the winter of 2018 (shown in red). Continued investigation with simulations will be crucial for a deeper understanding of these dynamics but is beyond the scope of this paper.

The potential of using the EWR and RV curves to study spot lifetimes is promising, as significant changes in the size of the IQR are observed from season to season in our TTS sample. Additionally, the shape of the EWR and RV folded curves varies dramatically from season to season in several of our TTSs. By observing over a longer period, we should gain a clearer understanding of the lifetimes of spots on these young active stars.

Periodograms, as well as EWRs and RVs folded at the rotation period, for our TTS targets with at least six observations in two or more seasons (Hubble 4 is not shown given its complex nature, see Section <ref>). Each data point is plotted twice as we are showing two phases. In the left column are the WTTSs (V827 Tau, V830 Tau, and V1075 Tau) and in the right are the CTTSs (AA Tau, DK Tau, and CI Tau). Periodograms are calculated from RVs (red), LDRs (blue), and EWRs (black). The strongest periodogram peak for all three datasets is marked with a solid triangle, and the second strongest peak is marked with an open triangle if significant. The adopted rotation period from the EWR data is displayed on the upper right of each panel. For comparison, rotation periods derived from lightcurve analysis, e.g., K2 [Rebull et al., 2020] and AAVSO [Watson et al., 2006], are printed as text and also indicated with purple carets above the periodogram. The 1% false alarm probabilities (FAP) estimated using the bootstrap method are indicated with horizontal dotted lines. For each target, the middle and bottom panels correspond to the folded EWR and RV curves (to the EWR period), respectively; data points are color-coded by date (scale bars on right). The $\Delta$ phase values indicate the phase delay between the EWR and RV curves, as detailed in Section <ref>. The sine waves used to calculate $\Delta$phases are shown in lines with matching color to the adopted season (Section <ref>). $^*$ Final results for CI Tau are presented in Figure <ref>.

§.§ A New Way to Estimate Stellar Rotation Period?

In Figure <ref>, for the WTTSs, all datasets (EWRs, LDRs, and RVs) display the strongest periodicity at the same period determined by lightcurve analysis to within about 0.02 d. This suggests that, even with our limited samples, power spectrum analysis of EWRs and LDRs have similar sensitivity to rotation periods as photometry for WTTSs <cit.>. We show the adopted  from the EWR data on the upper right of each panel. The quoted period uncertainties are calculated using MC sampling of the data, with a simulated scatter equivalent to the standard deviation of the data around a sinusoidal fit. This effectively takes into account the influence of astrophysical variability on the rotation period determination process. Moreover, to avoid signals that are too different from the adopted , we only search for periods within $\pm$1 day of the adopted  for each sample in the MC process.

The periodograms for the CTTSs are much noisier than those of the WTTSs. Greater discrepancies between the periodicities can be seen in EWR, LDR, RV, and the  determined from photometric studies. Strong peaks around 1 day are discarded as they often resemble the nightly observation cadence. To fortify the results on CTTSs, we perform another period-searching algorithm that is designed to process non-sinusoidal and irregularly spaced time-series data, the phase dispersion minimization (PDM) technique [Stellingwerf, 1978], on the EWR data. Using the same MC method to estimate the uncertainties, we find the following  values: for AA Tau, the LSP method yields 8.086 $\pm$ 0.245 d, while the PDM method yields 8.085 $\pm$ 0.292 d; for DK Tau, the LSP method gives 8.150 $\pm$ 0.025 d, and the PDM method gives 8.151 $\pm$ 0.041 d; for CI Tau, the LSP method results in 9.324 d, and the PDM method results in 9.017 d. For CI Tau, the strongest signal (9.324 d) from the LSP does not reveal any periodicity as shown in Figure <ref>. Thus, the uncertainties of its EWR  are not reported.

For both AA Tau and DK Tau, the  values determined from the LSP and PDM methods are nearly identical, with only a 0.001 d difference. However, both methods provide slightly different uncertainties. We choose to adopt the higher values to be conservative.

In Figure <ref>, DK Tau's EWR, LDR, and RV data are all periodic at $\sim$8.15 $\pm$ 0.15 d, about 0.31 days longer than the period from the K2 data analysis [Rebull et al., 2020]. However, our EWR  of 8.15 $\pm$ 0.041 days agrees well with AAVSO's result of  $\sim$ 8.18 days [Watson et al., 2006]. For AA Tau, the EWRs are periodic at $\sim$ 8.086 d, the LDRs at $\sim$ 8.077 d, and the RVs at $\sim$ 8.288 d. These periods are, on average, $\sim$0.1 days different from the  determined from AAVSO photometric data <cit.>, which is within the uncertainty of EWR's , 0.292 d.

In contrast to the LSP of other targets shown in Figure <ref>, the periodogram of CI Tau shows less agreement between the signals dominating the RV, LDR, and EWR data. The strong periodicity of LDR and EWR near 9.324 $\pm$ 0.03 d is about 0.3 d different then the period determined by K2 data analysis $\sim$9.034 d [Rebull et al., 2020]. Folding the EWR and RV data at 9.324 d in the bottom two panels does not show any coherent signal. In search for an alternative signal, we pick next strongest peaks around 9.3 days for EWR, LDR, and RV, marked with a y symbol in Figure <ref>a, and list their corresponding periods in the upper right corner of the plot (RV at 8.994 days, LDR at 8.984 days, and EWR at 9.014 days). Interestingly, these periods are all very close to the  determined from the PDM above, 9.017 d. Figures <ref>b, <ref>c, and <ref>d illustrate the EWR, LDR, and RV data, respectively, folded to the average period of these four signals, 9.002 days. This new period caused the folded curves in Figures <ref>b–d to exhibit a more pronounced sine wave-like pattern, indicating a more uniform trend across different seasons compared to CI Tau's results in Figure <ref>. The reason that the best period that produces a coherent signal in all the data does not appear as the strongest peak in the LSP may be attributed to the scatter observed around the best-fit sine curve for each season in Figures <ref>b–d. Additionally, the mean values of EWR, LDR, and RV are shifted from season to season. Together, these factors contribute to the LSP's decreased sensitivity in accurately detecting the correct period, and indicate the preferability of the PDM approach under these conditions.

A re-analysis of CI Tau structured similarly to Figure <ref>. The y symbols in  the periodogram (panel a) identify strong peaks around 9.3 d that produce more coherent signals in the folded EWR, LDR, and RV data. The color-coded periods of the selected peaks are given just outside the upper right of panel (a). The two strong periods identified by Biddle et al., 2018 using K2 lightcurve analysis are shown in purple carets above the periodogram. Panels (b), (c), and (d) show EWR, LDR, and RV data folded to 9.002 d and color-coded by date (scale bar on right). A coherent sinewave-like modulation is clearer for CI Tau when folded to 9.002 d, as opposed to the 9.324 d period used in Figure <ref>. The magnitude of the phase delay ($\Delta$Phase) between the RV and EWR curves for four seasons (S1: 2014–2015, S3: 2016–2017, S4: 2017–2018, and S5: 2018–2019) are shown in panels (b) and (c) and the fitted sine curves for each season are plotted, color-coded to match the $\Delta$Phase text.

§.§ Phase Delay Between EWR and RV Curves

If RV modulation originates from spot(s), we would expect to see a phase delay between the RV and the LDR/EWR curve. This is because the RV extrema occur when spot coverage on the stellar surface is at the most asymmetrical phases (when spot(s) are closest to each of the limbs). As for the LDR and EWR, in the case of OH/Fe, the signal's maximum occurs during the greater cool spot(s) coverage phase, and the minimum occurs during the lower cool spot(s) coverage phase. Here, we define a positive delay when the RV curve leads the EWR curve. Depending on the exact spot configuration, the phase delay between the RV and LDR/EWR curves can vary, but is expected to be a quarter-phase for a single dominant cool spot. As long as spots dominate the RV signal, the phase delay should not be near zero.

We fit a sine function to the EWR and RV curves to measure this phase offset (Fig. <ref>). To obtain the best fitting results for each target, we only fit data from the same observation season for which we have dense sampling. For WTTSs, we fit data from 2016–2019; for DK Tau, we fit from 2017–2019, and AA Tau we only fit for the last season of data, 2018–2019. The phase differences between the EWR and RV curves are calculated with uncertainties estimated via the MC method. The weighted results among seasons are given in Figure <ref> as $\Delta$Phase. The $\Delta$Phase, excluding V830 Tau, ranges from 0.20 to 0.32 with a typical uncertainty of 0.1. The average of these targets' $\Delta$Phase is $\sim$0.28, close to a quarter-phase delay. Although V830 Tau also shows a significant phase delay between the RV and the EWR curves, its $\Delta$Phase, about 0.68, is almost half a phase larger than the targets mentioned before. One possibility for this result is for V830 Tau to be dominated by cool regions and have one significant relatively hot spot. Such a configuration would produce RV and EWR curves with the phase delay as defined here equal to $\sim 0.75$. Future work on modeling can help confirm this hypothesis and can also help explore more potential scenarios.

§.§ Possible Commensurability of : = 1:1 for CI Tau

CI Tau's $\Delta$Phase significantly deviates from the expected quarter-phase delay, showing almost no delay during Season 3 (S3, 2016–2017) and Season 5 (S5, 2018–2019), as shown in Figure <ref>. Season 4 (S4, 2017–2018) exhibits a marginally larger $\Delta$Phase in LDR, approximately 0.1, yet it still lacks significance relative to a quarter-phase delay. Notably, for both LDR and EWR, Season 1 (S1, 2014–2015) has the greatest $\Delta$Phase, around 0.4. Although the modulation of EWR's S1 data might not be convincing because of the small variation magnitude, the LDR data surprisingly follows the best-fit sine curve. This result in S1 may be attributed to the significantly lower number of data points, only half as many as in S3 and S5, but if it is real, this change in the value of $\Delta$Phase can indicate a migrating of spot(s).

The large uncertainties in the $\Delta$Phase value mainly reflect the RV curves' small amplitude. With each season's data points tightly following the best-fit sine curve for both the EWR and LDR data, we consider the near-zero phase delay in CI Tau to be genuine. The weighted average $\Delta$Phase is 0.29 $\pm$ 0.07, and by excluding S1, we get 0.06 $\pm$ 0.13. Therefore, CI Tau stands out from other targets shown in Figure <ref> by demonstrating a near zero, non-quarter phase delay.

This lack of a quarter-phase delay between the EWR and the RV curve might indicate that CI Tau's RV signal is not solely caused by spot(s). A candidate hot Jupiter (HJ) companion to CI Tau with an  of $\sim$9 d, was identified by the YESS team based on optical and near-IR RVs [Johns-Krull et al., 2016]. The $\sim$11.6 $M_{\rm Jup}$ HJ candidate was later confirmed by Flagg et al., 2019 through a direct detection of CO in the planet's atmosphere, modulated by a period of $\sim$8.9891 d. Moreover, this 9 d signal was also present in the K2 lightcurve, and Biddle et al., 2018 associates this signal with the planet-disk interaction. However, the optical spectropolarimetric analysis of Donati et al., 2020 led the authors to suggest that the 9 d RV modulation was instead associated with stellar activity.

A recent study by the GRAVITY team [Gravity Collaboration et al., 2023] using the Very Large Telescope Interferometer (VLTI) provided an updated look at the CI Tau system on inner sub-au scales. They concluded that a massive inner planet orbiting at $\sim$0.08 au can explain the extension of the inner disk edge compared to the dust sublimation radius. This orbital radius of 0.08 au is approximately the stellar corotation radius of CI Tau, when estimated with a stellar rotation period of 9 d and a stellar mass of 0.9 $M_{\odot}$ <cit.>. Additionally, Kozdon et al., 2023's spectroscopic analysis of data taken at the NASA IRTF revealed a nine-day period in the accretion variability measured with the hydrogen Pf$\beta$ line, hinting at companion-driven accretion dynamics. The discrepancy between the inner disk's inclination ($\sim$70 degrees) and the outer disk's ($\sim$49 degrees) further supports the existence of a close-in massive planet [Gravity Collaboration et al., 2023, Clarke et al., 2018]. As the measurement of EWR is not sensitive
to accretion, our independent determination of CI Tau's rotation period at approximately 9.002 days using EWR (Section <ref>),  suggests that CI Tau and CI Tau b might be in a : = 1:1 commensurability, a result of the in situ planet formation and/or stalling at the corotation radius during disk migration [Dawson & Johnson, 2018].

Manick et al., 2024 recently identified a potential planet of 3.6 $\pm$ 0.3 Jupiter masses orbiting CI Tau with a  of about 25.2 d, a finding consistently observed across all their photometric and spectroscopic data. Concurrently, Donati et al., 2024 observed a similar periodic signal, but determined it to be 23.86 d, noting it was distinct in CO lines but absent in atomic lines. Consequently, Donati et al., 2024 suggested this signal likely corresponds to a non-axisymmetric structure within CI Tau's inner disk, rather than a planetary signal. Our periodogram analysis of CI Tau, shown in Figure <ref>a, aligns with these findings; while we do not detect a signal around 25.2 d, we do identify a peak around 23.86 d in the LDR and EWR data, although the false alarm probability is 76%. Interestingly, Biddle et al., 2021 had previously noted a signal around 24 d, specifically at 24.44 d, suggesting an earlier observation of this phenomenon. The absence of an RV signal at 23.86 d in our data merits further investigation, especially considering that  measurements heavily depend on CO lines in the K band. However, exploring this result falls outside of the scope of this paper.

§.§ Dependence of the Effective Temperature Variation

To explore whether the magnitude of the  variations correlates with parameters such as , average surface magnetic field strength ($\bar B$), $K$ band veiling (r$_{\rm K}$), and , we plot each target's median absolute deviation (MAD) of their EWR-derived  against these four properties (Fig. <ref>). Values of , $\bar B$, r$_{\rm K}$, and  are from López-Valdivia et al., 2021 (Table <ref>). The MAD of the  is meant to represent the level of the average stellar activity. Therefore, for targets with multi-season observation, the  MAD shown is the average of each season's values.

We use a Kendall rank correlation coefficient to test the dependence between parameters, as we have no assumptions about their underlying relationships. We calculate the $\tau$ and p-values for the full sample (13 TTSs) as well as for just the CTTSs. For a sample size of 13 (full sample), the critical $\tau_c$ is 0.359 at $\alpha=0.05$, and for a sample size of 9 (CTTSs only) $\tau_c$ is 0.5 at $\alpha=0.05$. As shown in Figure <ref>, only  for the CTTS sample has a $\tau$ value, 0.54, larger than the critical value of 0.5, i.e., that shows a significant correlation at the 95% confidence level; this is also evident from the corresponding p-value of 0.05. However, given that  is a function of both the stellar inclination angle and the stellar rotation velocity, interpreting this correlation becomes challenging. This complexity further motivates our investigation into the correlation between the MAD of  and the inclination angle, as detailed below.

The median absolute deviation (MAD) of  for the TTS sample as a function of (a) , (b) average surface magnetic field strength ($\bar B$), (c) $K$ band veiling (r$_{\rm K}$), and (d) .
Red corresponds to CTTSs, blue to WTTSs.
The Kendall rank correlation coefficient ($\tau$) and associated p-values

are given as determined for all objects (black) and for CTTSs alone (red). Measurement uncertainties are shown for , $\bar B$, r$_{\rm K}$, and . The y-axis errorbars show the standard deviations of the  variation amplitudes from different seasons, so not all targets have y-axis errorbars.

Figure <ref> shows EWR 's MAD against the stellar inclination angle ($i$). For WTTSs, we obtain $i$ from the literature, if possible, such as for V830 Tau from Donati et al., 2015 and for Hubble 4 from Carvalho et al., 2021. Both studies estimate the $i$ from  using estimates of stellar radius ($R_\ast$) and measurements of stellar rotation period ($P_{\rm rot}$). For V827 Tau, we calculated an $i \sim$ 41.1 deg from  = 20.8  (Table <ref>), $P_{\rm rot} = 3.758$ d [Rebull et al., 2020], and $R_\ast = 2.35$ $R_\odot$ [Gangi et al., 2022]. We leave V1075 Tau out of Figure <ref> as we cannot find a dedicated study which determines a reliable $R_\ast$ for this star.

For CTTSs, we assume that the circumstellar disk is perpendicular to the stellar rotation axis and adopt the disk inclination angle as the stellar inclination angle. Inner disk $i$ is preferred, if available, as it is more representative of the stellar rotation $i$ than outer disk $i$ [Nelissen et al., 2023]. Disk inclination angles collected from the literature are listed in Table <ref>, and we take averages for targets with multiple measurements.

The $\tau$ value for the full sample of 12 targets between the 's MAD to the stellar $i$ is 0.351 as shown in Figure <ref>. This value indicates a correlation between the two parameters at a 90% confidence level as $\tau_c = 0.303$ at $\alpha=0.1$ with a sample size of 12. Although the sample size is small, this relationship, if real, can be explained by TTSs hosting high latitude spot(s) <cit.>. With spot(s) at high latitude, the changes in the spot coverage area is minimized when the star is nearly pole on (small $i$), such that less temperature variability is expected to be observed. For larger $i$, we are closer to facing the equator of the TTS, and so spot coverage area will vary more over the visible hemisphere to produce a larger observed temperature variability.

The median absolute deviation (MAD) of  for all our TTS samples as a function of stellar inclination angle ($i$). For CTTSs (red symbols), the $i$ values are estimated either from the outer circumstellar disk's inclination angle (circles) or the inner disk's inclination angle (squares). Values for $i$ are estimated from , stellar radius, and stellar rotation period (triangles, see Section <ref>). The Kendall rank correlation coefficient ($\tau$) and the p-values are given for the full sample (black), and only for the CTTSs alone (red).

c L l
TTS Sample's Inclination Angles

Target            Incl.

3cOuter Disk Incl.

AA Tau           59                          Francis & van der Marel, 2020

59.1 ±0.3                Loomis et al., 2017

CI Tau           50.0 ±0.3                 Long et al., 2019

49.24                       Clarke et al., 2018

DH Tau           16.9 ^+2.0_-2.2         Long et al., 2019

48.4 ^+1.4_-1.5         Sheehan et al., 2019

21.4 ^+4.2_-3.0         Rota et al., 2022

DK Tau           12.8 ^+2.5_-2.8         Long et al., 2019

21.1 ^+3.3_-2.7         Rota et al., 2022

DM Tau           35                          Francis & van der Marel, 2020

35.2 ±0.7                Kudo et al., 2018

36.0 ^+0.12_-0.09       Flaherty et al., 2020

DS Tau           65.2 ±0.3                Long et al., 2019

70 ±3                    Long et al., 2019

GI Tau           43.8 ±1.1                Long et al., 2019

IQ Tau           62.1 ±0.5                Long et al., 2019

LkCa 15          55                          Francis & van der Marel, 2020

50.16 ±0.03              Facchini et al., 2020

43.95 ^+2.39_-2.06      Bohn et al., 2022

3cInner Disk Incl.

AA Tau           55 ±25               Francis & van der Marel, 2020

CI Tau           71 ±1               Gravity Collaboration et al., 2023

DK Tau           58 ^+18_-11         Nelissen et al., 2023

DM Tau           41 ±4               Francis & van der Marel, 2020

LkCa 15          61.02 ^+18.59_-20.80      Bohn et al., 2022

3cStellar Incl.

Hubble 4         11.8                    Carvalho et al., 2021

V827 Tau         39.9                    This study

V830 Tau         55 ±10               Donati et al., 2015

§ SUMMARY

As part of the YESS project, we used the atomic Fe i and molecular OH lines in the near-infrared $H$ band to investigate the spot variability of T Tauri Stars (TTSs). We demonstrate a positive correlation between the line depth of the OH lines at $\sim$1.56310  and $\sim$1.56317  as the effective temperature () drops from $\sim$5000 K to $\sim$3400 K, and a negative correlation between the Fe i lines at $\sim$1.56259  and $\sim$1.56362  as  decreases. These temperature sensitivities make the lines ideal for probing the  variability in TTSs within the 3400–5000 K range. Using data from 49 calibration targets observed with the Immersion GRating INfrared Spectrometer (IGRINS), we establish an empirical –equivalent width ratio (EWR) relationship. Key findings and contributions of this study include:

* We observe that line depth ratios (LDRs), heavily influenced by blending in our spectral analysis region, are susceptible to various broadening effects, including those resulting from rotational velocity (), spectral resolution, and magnetic field strength. In contrast, these influences have little to no impact on our EWR measurements.

* We verify the reliability of our EWR methodology through observations of the quiescent M-dwarf GJ 281. Across 61 nights spanning four years, the star exhibits a scatter in EWR of $\sigma_{\rm EWR} \approx 0.012$, equivalent to $\sigma_{T_{\rm eff}} \approx$ 12 K.

* A comparison between our empirical –EWR relationship and our model grid, particularly for solar metallicity targets in the $\sim$3400 K to 5000 K range, shows general agreement. However, there are noticeable discrepancies in equivalent width (EW) measurements between observed and modeled spectra. In order to extend the applicability of the EWR technique, refined physical parameters, such as metallicity, for the  calibration targets are needed, as well as improvements in spectral modeling and molecular line lists for cooler M-dwarfs.

* Our analysis of  variability in a sample of 13 TTSs reveals no noticeable trend distinguishing Weak-line TTSs and Classical TTSs, with both groups exhibiting substantial  variations ($\gtrsim$150 K) within a two-year observational period. This might indicate a common spot structure between these two types of TTSs. Furthermore, we discover a strong positive correlation between the median  of the TTSs to their  variation magnitudes, hinting a higher surface activity with higher average .

* We observe a quarter-phase delay between most of our targets' radial velocity (RV) and EWR modulations when folding them to their rotation periods. V830 Tau and CI Tau stand out for having non-quarter-phase delay results. The clean and stable $0.68$ phase delay seen in V830 Tau, over a period of 4 years, suggests the potential involvement of hot spot(s). With a phase delay of 0.06 $\pm$ 0.13, the CTTS CI Tau suggesting additional influences on its RVs and/or EWRs, possibly from a planetary companion, hinting at a 1:1 commensurability between CI Tau's rotation period and the orbital period of its Hot Jupiter companion.

* We report a positive correlation between  variability amplitude and stellar inclination angle with 90% confidence, in agreement with previous studies that have identified high-latitude spots on TTSs.

We are thankful for the comments and suggestions from an anonymous referee to help improve the quality of this paper. We also thank the technical and logistical staff at McDonald and Lowell Observatories for their excellent support of the Immersion Grating Infrared Spectrograph (IGRINS) installations, software, and observation program described here. In particular, D. Doss, C. Gibson, J. Kuehne, K. Meyer, B. Hardesty, F. Cornelius, M. Sweaton, J. Gehring, S. Zoonematkermani, E. Dunham, S. Levine, H. Roe, W. DeGroff, G. Jacoby, T. Pugh, A. Hayslip, and H. Larson. Partial support for this work was provided by NASA Exoplanet Research Program grant 80-NSSC19K-0289 (to LP). CMJ acknowledges partial support for this work through grants to Rice University provided by NASA (award 80-NSSC18K-0828) and the NSF (awards AST-2009197 and AST-1461918).

We are grateful for the generous donations of John and Ginger Giovale, the BF Foundation, and others which made the IGRINS-LDT program possible. Additional funding for IGRINS at the LDT was provided by the Mt. Cuba Astronomical Foundation and the Orr Family Foundation. IGRINS was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the US National Science Foundation under grant AST-1229522 and AST-1702267, of the University of Texas at Austin, and of the Korean GMT Project of KASI.

This work benefited from the language editing assistance of GPT-4, a language model developed by OpenAI in San Francisco, CA, USA <http://openai.com>. This work made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. This study also made use of the SIMBAD database and the VizieR catalogue access tool, both operated at CDS, Strasbourg, France. Some observations were obtained at the Lowell Discovery Telescope (LDT) at Lowell Observatory. Lowell is a private, non-profit institution dedicated to astrophysical research and public appreciation of astronomy and operates the LDT in partnership with Boston University, the University of Maryland, the University of Toledo, Northern Arizona University and Yale University. We have used IGRINS archival data older than the 2 year proprietary period.

LDT (IGRINS), Smith (IGRINS)

[Astropy Collaboration et al., 2013, Astropy Collaboration et al., 2018, Astropy Collaboration et al., 2022],
[Hunter, 2007],
, [Harris et al., 2020]
[Stahl et al., 2021, Tang et al., 2021],
[STScI Development Team, 2013]

§ ATMOSPHERE PROFILE COMPARISON

This appendix section describes the effects of employing different atmospheric profiles to determine equivalent widths (EWs) of the atomic Fei and molecular OH line region using the synmast code. The results presented here use the modified VALD3 line data tuned on GJ 205 (see Section <ref>) and supplemented with BT2 water lines. We investigate three distinct atmospheric profiles: BT-NextGen [Allard et al., 2011, Allard et al., 2012], PHOENIX-ACES[Atmospheric profiles retrieved from <ftp://phoenix.astro.physik.uni-goettingen.de/AtmosFITS/> under the v2.0 directory.], and MARCS [Gustafsson et al., 2008][Accessed via <https://marcs.astro.uu.se>]. The elemental abundances for the BT-NextGen and PHOENIX-ACES profiles are based on the ASGG2009 standard [Asplund et al., 2009], whereas the MARCS profiles conform to GAS2007 [Grevesse et al., 2007]. For all considered models, the scaling of $\alpha$-element abundances is consistent with the methodologies described in Section <ref>.

As depicted in Figure <ref>, the comparison of EW measurements across different atmospheric profiles reveals minimal variance for both the Fe and OH line regions. The choice of atmospheric model exhibits negligible influence on the precision of derived physical parameters for the targets. This consistency underscores the robustness of the employed spectroscopic techniques in yielding reliable stellar parameters irrespective of the chosen atmospheric profile.

Comparison of Fe and OH equivalent widths from the synmast model spectra using different atmospheric profiles and IGRINS data for the  calibration targets. Results derived from the PHOENIX ACES atmospheric profile are presented in the first column, while those obtained from the BT-NextGen and MARCS profiles are displayed in the second and third columns, respectively. Panels (a), (b), and (c) in the top row illustrate the results for Fe EW, and panels (d), (e), and (f) in the bottom row detail the OH EW findings. The color coding of data points corresponds to the effective temperature () scale, as indicated by the color bar on the right. The dashed line across each panel denotes the line of unity.

§ WATER LINE LIST COMPARISON

Figure <ref> illustrates the negligible impact on the equivalent width (EW) measurements of Fe (a) and OH (b) regions from employing different water ($^1{\rm H}_2$$^{16}{\rm O}$) line lists —BT2 [Barber et al., 2006] and POKAZATEL [Polyansky et al., 2018]— within the synmast spectral synthesis code. The differences, as shown, are minimal, with only a marginal shift in the Fe EW at an overall 6% difference.

POKAZATEL [Polyansky et al., 2018] and BT2 [Barber et al., 2006] water line list comparisons for Fe EW in (a) and OH EW in (b). All models are generated using synmast with an unmodified VALD3 line list and adopting the BT-NextGen atmosphere profile. Circles show results generated based on the  calibration targets' parameters, color-coded by their .

[Afşar et al., 2023]
Afşar, M., Bozkurt, Z., Topcu, G. B., et al. 2023, The Astrophysical
Journal, 949, 86, 10.3847/1538-4357/acc946

[Allard et al., 2011]
Allard, F., Homeier, D., & Freytag, B. 2011, 448, 91.

[Allard et al., 2012]
Allard, F., Homeier, D., Freytag, B., & Sharp, C. M. 2012, 57, 3,

[Anderson & Francis, 2011]
Anderson, E., & Francis, C. 2011, VizieR Online Data Catalog, V/137.

[Asplund et al., 2009]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of
Astronomy and Astrophysics, 47, 481,

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
Astronomy and Astrophysics, 558, A33, 10.1051/0004-6361/201322068

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, The Astronomical Journal, 156, 123, 10.3847/1538-3881/aabc4f

[Astropy Collaboration et al., 2022]
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, The
Astrophysical Journal, 935, 167, 10.3847/1538-4357/ac7c74

[Baines et al., 2021]
Baines, E. K., Thomas Armstrong, J., Clark, J. H., et al. 2021, The
Astronomical Journal, 162, 198, 10.3847/1538-3881/ac2431

[Baliunas et al., 1997]
Baliunas, S. L., Henry, G. W., Donahue, R. A., Fekel, F. C., & Soon, W. H.
1997, The Astrophysical Journal, 474, L119, 10.1086/310442

[Barber et al., 2006]
Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, Monthly
Notices of the Royal Astronomical Society, 368, 1087,

[Benatti et al., 2020]
Benatti, S., Damasso, M., Desidera, S., et al. 2020, A&A, 639, A50,

[Bevington & Robinson, 1992]
Bevington, P. R., & Robinson, D. K. 1992, Data Reduction and Error Analysis
for the Physical Sciences.

[Biddle et al., 2018]
Biddle, L. I., Johns-Krull, C. M., Llama, J., Prato, L., & Skiff, B. A.
2018, ApJ, 853, L34, 10.3847/2041-8213/aaa897

[Biddle et al., 2021]
Biddle, L. I., Llama, J., Cameron, A., et al. 2021, The Astrophysical
Journal, 906, 113, 10.3847/1538-4357/abc889

[Bohn et al., 2022]
Bohn, A. J., Benisty, M., Perraut, K., et al. 2022, Astronomy and
Astrophysics, 658, A183, 10.1051/0004-6361/202142070

[Bonomo et al., 2017]
Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, Astronomy and
Astrophysics, 602, A107, 10.1051/0004-6361/201629882

[Borsa et al., 2015]
Borsa, F., Scandariato, G., Rainer, M., et al. 2015, Astronomy and
Astrophysics, 578, A64, 10.1051/0004-6361/201525741

[Boyajian et al., 2012]
Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, The
Astrophysical Journal, 757, 112, 10.1088/0004-637X/757/2/112

[Boyajian et al., 2013]
—. 2013, The Astrophysical Journal, 771, 40,

[Brogi et al., 2012]
Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502,

[Cale et al., 2021]
Cale, B. L., Reefe, M., Plavchan, P., et al. 2021, AJ, 162, 295,

[Carroll et al., 2012]
Carroll, T. A., Strassmeier, K. G., Rice, J. B., & Künstler, A. 2012,
Astronomy and Astrophysics, 548, A95, 10.1051/0004-6361/201220215

[Carvalho et al., 2021]
Carvalho, A., Johns-Krull, C. M., Prato, L., & Anderson, J. 2021, ApJ, 910,
33, 10.3847/1538-4357/abe237

[Casagrande et al., 2014]
Casagrande, L., Portinari, L., Glass, I. S., et al. 2014, Monthly Notices of
the Royal Astronomical Society, 439, 2060, 10.1093/mnras/stu089

[Catalano et al., 2002]
Catalano, S., Biazzo, K., Frasca, A., & Marilli, E. 2002, Astronomy and
Astrophysics, 394, 1009,

[Clarke et al., 2018]
Clarke, C. J., Tazzari, M., Juhasz, A., et al. 2018, ApJ, 866, L6,

[Claytor et al., 2024]
Claytor, Z. R., van Saders, J. L., Cao, L., et al. 2024, The Astrophysical
Journal, 962, 47, 10.3847/1538-4357/ad159a

[Crockett et al., 2012]
Crockett, C. J., Mahmud, N. I., Prato, L., et al. 2012, ApJ, 761, 164,

[Cushing et al., 2005]
Cushing, M. C., Rayner, J. T., & Vacca, W. D. 2005, The Astrophysical Journal,
623, 1115, 10.1086/428040

[Dawson & Johnson, 2018]
Dawson, R. I., & Johnson, J. A. 2018, Annu. Rev. Astron. Astrophys., 56, 175,

[Dittmann et al., 2009]
Dittmann, J. A., Close, L. M., Green, E. M., & Fenwick, M. 2009, The
Astrophysical Journal, 701, 756, 10.1088/0004-637X/701/1/756

[Donati et al., 2014]
Donati, J. F., Hébrard, E., Hussain, G., et al. 2014, Monthly Notices of
the Royal Astronomical Society, 444, 3220, 10.1093/mnras/stu1679

[Donati et al., 2015]
Donati, J.-F., Hébrard, E., Hussain, G. A. J., et al. 2015, Mon. Not. R.
Astron. Soc., 453, 3707, 10.1093/mnras/stv1837

[Donati et al., 2017]
Donati, J.-F., Yu, L., Moutou, C., et al. 2017, Mon. Not. R. Astron. Soc.,
465, 3343, 10.1093/mnras/stw2904

[Donati et al., 2020]
Donati, J.-F., Bouvier, J., Alencar, S. H., et al. 2020, MNRAS, 491, 5660,

[Donati et al., 2024]
Donati, J.-F., Finociety, B., Cristofari, P. I., et al. 2024, The Classical
T Tauri Star CI Tau Observed with SPIRou: Magnetospheric
Accretion and Planetary Formation,  arXiv.

[D'Orazi et al., 2011]
D'Orazi, V., Biazzo, K., & Randich, S. 2011, Astronomy and Astrophysics, 526,
A103, 10.1051/0004-6361/201015616

[Earl et al., 2022]
Earl, N., Tollerud, E., Jones, C., et al. 2022, Astropy/Specutils:
V1.7.0, Zenodo, 10.5281/zenodo.6207491

[Endl et al., 2003]
Endl, M., Cochran, W. D., Tull, R. G., & MacQueen, P. J. 2003, The
Astronomical Journal, 126, 3099, 10.1086/379137

[Facchini et al., 2020]
Facchini, S., Benisty, M., Bae, J., et al. 2020, Astronomy and Astrophysics,
639, A121, 10.1051/0004-6361/202038027

[Flagg et al., 2019]
Flagg, L., Johns-Krull, C. M., Nofi, L., et al. 2019, ApJ, 878, L37,

[Flaherty et al., 2020]
Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, The Astrophysical
Journal, 895, 109, 10.3847/1538-4357/ab8cc5

[Flores et al., 2019]
Flores, C., Connelley, M. S., Reipurth, B., & Boogert, A. 2019, The
Astrophysical Journal, 882, 75, 10.3847/1538-4357/ab35d4

[Francis & van der Marel, 2020]
Francis, L., & van der Marel, N. 2020, The Astrophysical Journal, 892, 111,

[Fukue et al., 2015]
Fukue, K., Matsunaga, N., Yamamoto, R., et al. 2015, ApJ, 812, 64,

[Gangi et al., 2022]
Gangi, M., Antoniucci, S., Biazzo, K., et al. 2022, Astronomy and
Astrophysics, 667, A124, 10.1051/0004-6361/202244042

[Garnett, 2023]
Garnett, R. 2023, Bayesian Optimization (Cambridge University Press),

[Ghezzi et al., 2010]
Ghezzi, L., Cunha, K., Smith, V. V., & de la Reza, R. 2010, The
Astrophysical Journal, 724, 154, 10.1088/0004-637X/724/1/154

[Gravity Collaboration et al., 2023]
Gravity Collaboration, Soulain, A., Perraut, K., et al. 2023, Astronomy and
Astrophysics, 674, A203, 10.1051/0004-6361/202346446

[Gray, 1989]
Gray, D. F. 1989, The Astrophysical Journal, 347, 1021, 10.1086/168192

[Gray, 1994]
—. 1994, Publications of the Astronomical Society of the Pacific, 106, 1248,

[Gray & Baliunas, 1995]
Gray, D. F., & Baliunas, S. L. 1995, The Astrophysical Journal, 441, 436,

[Gray et al., 1992]
Gray, D. F., Baliunas, S. L., Lockwood, G. W., & Skiff, B. A. 1992, The
Astrophysical Journal, 400, 681, 10.1086/172030

[Gray & Brown, 2001]
Gray, D. F., & Brown, K. 2001, Publications of the Astronomical Society of the
Pacific, 113, 723, 10.1086/320811

[Gray & Johanson, 1991]
Gray, D. F., & Johanson, H. L. 1991, PASP, 103, 439, 10.1086/132839

[Grevesse et al., 2007]
Grevesse, N., Asplund, M., & Sauval, A. J. 2007, Space Science Reviews, 130,
105, 10.1007/s11214-007-9173-7

[Gully-Santiago et al., 2017]
Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836,
200, 10.3847/1538-4357/836/2/200

[Gustafsson et al., 2008]
Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, Astronomy and
Astrophysics, 486, 951, 10.1051/0004-6361:200809724

[Harris et al., 2020]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature,
585, 357, 10.1038/s41586-020-2649-2

[Haywood et al., 2014]
Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, Monthly Notices
of the Royal Astronomical Society, 443, 2517, 10.1093/mnras/stu1320

[Heiter et al., 2015]
Heiter, U., Jofré, P., Gustafsson, B., et al. 2015, Astronomy and
Astrophysics, 582, A49, 10.1051/0004-6361/201526319

[Hunter, 2007]
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90, 10.1109/MCSE.2007.55

[Husser et al., 2013]
Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553,
A6, 10.1051/0004-6361/201219058

[Jian et al., 2019]
Jian, M., Matsunaga, N., & Fukue, K. 2019, Monthly Notices of the Royal
Astronomical Society, 485, 1310, 10.1093/mnras/stz237

[Jian et al., 2020]
Jian, M., Taniguchi, D., Matsunaga, N., et al. 2020, Monthly Notices of the
Royal Astronomical Society, 494, 1724, 10.1093/mnras/staa834

[Jofré et al., 2014]
Jofré, P., Heiter, U., Soubiran, C., et al. 2014, Astronomy and
Astrophysics, 564, A133, 10.1051/0004-6361/201322440

[Johns-Krull, 2007]
Johns-Krull, C. M. 2007, ApJ, 664, 975, 10.1086/519017

[Johns-Krull et al., 2004]
Johns-Krull, C. M., Valenti, J. A., & Saar, S. H. 2004, ApJ, 617, 1204,

[Johns-Krull et al., 2016]
Johns-Krull, C. M., McLane, J. N., Prato, L., et al. 2016, ApJ, 826, 206,

[Kochukhov, 2007]
Kochukhov, O. 2007, arXiv:astro-ph/0701084.

[Kochukhov et al., 2010]
Kochukhov, O., Makaganiuk, V., & Piskunov, N. 2010, Astronomy and
Astrophysics, 524, A5, 10.1051/0004-6361/201015429

[Kovtyukh et al., 2023]
Kovtyukh, V., Lemasle, B., Nardetto, N., et al. 2023, Monthly Notices of the
Royal Astronomical Society, 523, 5047, 10.1093/mnras/stad1708

[Kozdon et al., 2023]
Kozdon, J., Brittain, S. D., Fung, J., et al. 2023, The Astronomical Journal,
166, 119, 10.3847/1538-3881/ace903

[Kraus et al., 2011]
Kraus, A. L., Ireland, M. J., Martinache, F., & Hillenbrand, L. A. 2011, ApJ,
731, 8, 10.1088/0004-637X/731/1/8

[Kudo et al., 2018]
Kudo, T., Hashimoto, J., Muto, T., et al. 2018, The Astrophysical Journal,
868, L5, 10.3847/2041-8213/aaeb1c

[Lanza, 2022]
Lanza, A. F. 2022, Astronomy and Astrophysics, 658, A195,

[Lanza et al., 2009]
Lanza, A. F., Aigrain, S., Messina, S., et al. 2009, A&A, 506, 255,

[Lee et al., 2017]
Lee, J.-J., Gullikson, K., & Kaplan, K. 2017, Igrins/Plp 2.2.0, Zenodo,

[Lomb, 1976]
Lomb, N. R. 1976, Astrophysics and Space Science, 39, 447,

[Long et al., 2019]
Long, F., Herczeg, G. J., Harsono, D., et al. 2019, The Astrophysical
Journal, 882, 49, 10.3847/1538-4357/ab2d2d

[Loomis et al., 2017]
Loomis, R. A., Öberg, K. I., Andrews, S. M., & MacGregor, M. A. 2017, The
Astrophysical Journal, 840, 23, 10.3847/1538-4357/aa6c63

[López-Valdivia et al., 2019]
López-Valdivia, R., Mace, G. N., Sokal, K. R., et al. 2019, The
Astrophysical Journal, 879, 105, 10.3847/1538-4357/ab2129

[López-Valdivia et al., 2021]
López-Valdivia, R., Sokal, K. R., Mace, G. N., et al. 2021, ApJ, 921,
53, 10.3847/1538-4357/ac1a7b

[Luck & Heiter, 2006]
Luck, R. E., & Heiter, U. 2006, The Astronomical Journal, 131, 3069,

[Luhman et al., 2017]
Luhman, K. L., Mamajek, E. E., Shukla, S. J., & Loutrel, N. P. 2017, The
Astronomical Journal, 153, 46, 10.3847/1538-3881/153/1/46

[Mace et al., 2016]
Mace, G., Kim, H., Jaffe, D. T., et al. 2016, in SPIE Astronomical
Telescopes + Instrumentation, ed. C. J. Evans, L. Simard, & H. Takami,
Edinburgh, United Kingdom, 99080C, 10.1117/12.2232780

[Maldonado et al., 2015]
Maldonado, J., Affer, L., Micela, G., et al. 2015, Astronomy and
Astrophysics, 577, A132, 10.1051/0004-6361/201525797

[Manick et al., 2024]
Manick, R., Sousa, A. P., Bouvier, J., et al. 2024, Long Period Modulation of
the Classical T Tauri Star CI Tau: Evidence for an Eccentric Close-in
Massive Planet at 0.17 Au,  arXiv.

[Mann et al., 2013]
Mann, A. W., Brewer, J. M., Gaidos, E., Lépine, S., & Hilton, E. J.
2013a, The Astronomical Journal, 145, 52,

[Mann et al., 2014]
Mann, A. W., Deacon, N. R., Gaidos, E., et al. 2014, The Astronomical
Journal, 147, 160, 10.1088/0004-6256/147/6/160

[Mann et al., 2015]
Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015,
The Astrophysical Journal, 804, 64, 10.1088/0004-637X/804/1/64

[Mann et al., 2013]
Mann, A. W., Gaidos, E., & Ansdell, M. 2013b, The Astrophysical
Journal, 779, 188, 10.1088/0004-637X/779/2/188

[Nelissen et al., 2023]
Nelissen, M., McGinnis, P., Folsom, C. P., et al. 2023, Astronomy and
Astrophysics, 670, A165, 10.1051/0004-6361/202245194

[Neves et al., 2012]
Neves, V., Bonfils, X., Santos, N. C., et al. 2012, Astronomy and
Astrophysics, 538, A25, 10.1051/0004-6361/201118115

[O'Neal & Neff, 1997]
O'Neal, D., & Neff, J. E. 1997, The Astronomical Journal, 113, 1129,

[O'Neal et al., 2001]
O'Neal, D., Neff, J. E., Saar, S. H., & Mines, J. K. 2001, The Astronomical
Journal, 122, 1954, 10.1086/323093

[Park et al., 2014]
Park, C., Jaffe, D. T., Yuk, I.-S., et al. 2014, in SPIE Astronomical
Telescopes + Instrumentation, ed. S. K. Ramsay, I. S. McLean, &
H. Takami, Montréal, Quebec, Canada, 91471D, 10.1117/12.2056431

[Parks et al., 2014]
Parks, J. R., Plavchan, P., White, R. J., & Gee, A. H. 2014, The Astrophysical
Journal Supplement Series, 211, 3, 10.1088/0067-0049/211/1/3

[Passegger et al., 2022]
Passegger, V. M., Bello-García, A., Ordieres-Meré, J., et al.
2022, Astronomy and Astrophysics, 658, A194,

[Piskunov, 1999]
Piskunov, N. 1999, 243, 515, 10.1007/978-94-015-9329-8_45

[Polyansky et al., 2018]
Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, Monthly Notices
of the Royal Astronomical Society, 480, 2597, 10.1093/mnras/sty1877

[Prato et al., 2008]
Prato, L., Huerta, M., Johns-Krull, C. M., et al. 2008, ApJ, 687, L103,

[Prato et al., 2002]
Prato, L., Simon, M., Mazeh, T., et al. 2002, ApJ, 569, 863,

[Ramírez et al., 2012]
Ramírez, I., Fish, J. R., Lambert, D. L., & Allende Prieto, C. 2012, The
Astrophysical Journal, 756, 46, 10.1088/0004-637X/756/1/46

[Rampalli et al., 2021]
Rampalli, R., Agüeros, M. A., Curtis, J. L., et al. 2021, The
Astrophysical Journal, 921, 167, 10.3847/1538-4357/ac0c1e

[Rayner et al., 2009]
Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, The Astrophysical Journal
Supplement Series, 185, 289, 10.1088/0067-0049/185/2/289

[Rebull et al., 2020]
Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2020, The Astronomical
Journal, 159, 273, 10.3847/1538-3881/ab893c

[Reiners et al., 2022]
Reiners, A., Shulyak, D., Käpylä, P. J., et al. 2022, Astronomy and
Astrophysics, 662, A41, 10.1051/0004-6361/202243251

[Rizzuto et al., 2020]
Rizzuto, A. C., Dupuy, T. J., Ireland, M. J., & Kraus, A. L. 2020, The
Astrophysical Journal, 889, 175, 10.3847/1538-4357/ab5aed

[Rojas-Ayala et al., 2012]
Rojas-Ayala, B., Covey, K. R., Muirhead, P. S., & Lloyd, J. P. 2012, The
Astrophysical Journal, 748, 93, 10.1088/0004-637X/748/2/93

[Rota et al., 2022]
Rota, A. A., Manara, C. F., Miotello, A., et al. 2022, Astronomy and
Astrophysics, 662, A121, 10.1051/0004-6361/202141035

[Ryabchikova et al., 2015]
Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr., 90,
054005, 10.1088/0031-8949/90/5/054005

[Saar & Donahue, 1997]
Saar, S. H., & Donahue, R. A. 1997, The Astrophysical Journal, 485, 319,

[Sawczynec et al., 2022]
Sawczynec, E., Mace, G., Gully-Santiago, M., & Jaffe, D. 2022, 54, 203.06.

[Sawczynec et al., 2023]
—. 2023, 55, 207.14.

[Scargle, 1982]
Scargle, J. D. 1982, The Astrophysical Journal, 263, 835,

[Sheehan et al., 2019]
Sheehan, P. D., Wu, Y.-L., Eisner, J. A., & Tobin, J. J. 2019, The
Astrophysical Journal, 874, 136, 10.3847/1538-4357/ab09f9

[Sikora et al., 2023]
Sikora, J., Rowe, J., Barat, S., et al. 2023, The Astronomical Journal, 165,
250, 10.3847/1538-3881/acc865

[Skrutskie et al., 2006]
Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, The Astronomical
Journal, 131, 1163, 10.1086/498708

[Sokal et al., 2020]
Sokal, K. R., Johns-Krull, C. M., Mace, G. N., et al. 2020, ApJ, 888, 116,

[Stahl et al., 2021]
Stahl, A. G., Tang, S.-Y., Johns-Krull, C. M., et al. 2021, AJ, 161, 283,

[Stellingwerf, 1978]
Stellingwerf, R. F. 1978, The Astrophysical Journal, 224, 953,

[STScI Development Team, 2013]
STScI Development Team. 2013, Astrophysics Source Code Library,

[Tang et al., 2021]
Tang, S.-Y., Stahl, A., Johns-Krull, C., Prato, L., & Llama, J. 2021, JOSS,
6, 3095, 10.21105/joss.03095

[Tang et al., 2023]
Tang, S.-Y., Stahl, A. G., Prato, L., et al. 2023, The Astrophysical Journal,
950, 92, 10.3847/1538-4357/acc58b

[Taniguchi et al., 2018]
Taniguchi, D., Matsunaga, N., Kobayashi, N., et al. 2018, Monthly Notices of
the Royal Astronomical Society, 473, 4993, 10.1093/mnras/stx2691

[Tayar et al., 2022]
Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2022, The
Astrophysical Journal, 927, 31, 10.3847/1538-4357/ac4bbc

[Tennyson & Yurchenko, 2012]
Tennyson, J., & Yurchenko, S. N. 2012, Monthly Notices of the Royal
Astronomical Society, 425, 21, 10.1111/j.1365-2966.2012.21440.x

[Terrien et al., 2015]
Terrien, R. C., Mahadevan, S., Bender, C. F., Deshpande, R., & Robertson, P.
2015, The Astrophysical Journal, 802, L10,

[Tran et al., 2023]
Tran, Q. H., Bedell, M., Foreman-Mackey, D., & Luger, R. 2023, The
Astrophysical Journal, 950, 162, 10.3847/1538-4357/acd05c

[Valenti & Fischer, 2005]
Valenti, J. A., & Fischer, D. A. 2005, The Astrophysical Journal Supplement
Series, 159, 141, 10.1086/430500

[Valenti & Piskunov, 1996]
Valenti, J. A., & Piskunov, N. 1996, Astronomy and Astrophysics Supplement
Series, 118, 595.
<https://ui.adsabs.harvard.edu/abs/1996A   AS..118..595V>

[VanderPlas, 2018]
VanderPlas, J. T. 2018, The Astrophysical Journal Supplement Series, 236, 16,

[von Braun et al., 2011]
von Braun, K., Boyajian, T. S., Kane, S. R., et al. 2011, The Astrophysical
Journal, 729, L26, 10.1088/2041-8205/729/2/L26

[Wallace et al., 1996]
Wallace, L., Livingston, W., Hinkle, K., & Bernath, P. 1996, The Astrophysical
Journal Supplement Series, 106, 165, 10.1086/192333

[Watson et al., 2006]
Watson, C. L., Henden, A. A., & Price, A. 2006, Society for Astronomical
Sciences Annual Symposium, 25, 47.

[Wehrhahn et al., 2023]
Wehrhahn, A., Piskunov, N., & Ryabchikova, T. 2023, Astronomy and
Astrophysics, 671, A171, 10.1051/0004-6361/202244482

[White et al., 2018]
White, T. R., Huber, D., Mann, A. W., et al. 2018, Monthly Notices of the
Royal Astronomical Society, 477, 4403, 10.1093/mnras/sty898

[Yu et al., 2019]
Yu, L., Donati, J. F., Grankin, K., et al. 2019, Monthly Notices of the Royal
Astronomical Society, 489, 5556, 10.1093/mnras/stz2481

[Yuk et al., 2010]
Yuk, I.-S., Jaffe, D. T., Barnes, S., et al. 2010, in SPIE Astronomical
Telescopes + Instrumentation, ed. I. S. McLean, S. K. Ramsay, &
H. Takami, San Diego, California, USA, 77351M, 10.1117/12.856864