Rovibrational Spectroscopy of CI Tau—Evidence of a Multicomponent Eccentric Disk Induced by a Planet

The protoplanetary disks of T Tauri stars often host substructures such as rings/gaps (e.g., Long et al. 2018), spiral arms (e.g., Phuong et al. 2020), and large-scale asymmetries (e.g., Dyda et al. 2015; see Bae et al. 2022 for a review). These features may be produced by processes such as condensation fronts or snow lines (e.g., Pinilla et al. 2017), radiation pressure (e.g., Bi & Fung 2022), gravitational instabilities (e.g., Dong et al. 2018), and dynamical interactions with planets (Paardekooper et al. 2022, and references therein). Among these possibilities, planet–disk interactions have generated the most interest because of the connections with exoplanet discoveries.

Exoplanets are now known to be ubiquitous (e.g., Batalha 2014; van der Marel & Mulders 2021; Christiansen 2022). The majority of detected exoplanets were found by transit and radial velocity surveys that target mature stars well beyond the stage when gas-rich, planet-building disks are present. On the other hand, detections of planets accompanying T Tauri stars (TTSs) and other young stellar objects (YSOs) that still have their protoplanetary disks are quite rare. This disparity arises because YSOs exhibit substantial stellar activity that, alongside an optically thick disk, obscures and dampens planetary signatures. Techniques that characterize planets through their imprints on the protoplanetary disks can supplement traditional planet detection techniques.

PDS 70 b and c are two known protoplanets around a TTS (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019), while AB Aur b (Currie et al. 2022) and CI Tau b (Johns-Krull et al. 2016; Clarke et al. 2018) are two of the next likely candidates. Among these systems, CI Tau b is unique because it is a potential hot Jupiter that has been detected by Doppler monitoring whereas PDS 70 b and c and AB Aur b were detected by direct imaging at large separations. This allows for CI Tau b to have a well-constrained mass and orbital parameters (Johns-Krull et al. 2016; Flagg et al. 2019). As such, one should be able to apply theories of planet–disk interactions to observations of the inner disk of CI Tau, and directly test the theoretical predictions. For example, simulations indicate that a massive companion should induce an eccentricity in the disk (Kley & Dirksen 2006; Teyssandier & Ogilvie 2017) and that the disk drives an eccentricity on the planet (Papaloizou et al. 2001; Duffell & Chiang 2015; Rosotti et al. 2017; Ragusa et al. 2018; Muley et al. 2019).

CI Tau is a ∼2 Myr old TTS of spectral type K5.5 located at a distance of d = 160 ± 10 pc (Gaia Collaboration et al. 2023) with a mass of M_* = 1.02 ± 0.001 M_☉ (Law et al. 2022). Its hot Jupiter companion, CI Tau b (Johns-Krull et al. 2016), has a mass of 11.6 ± 2.8 M_J (Flagg et al. 2019) and an eccentricity of e = 0.25 ± 0.16 (Flagg et al. 2019). CI Tau b's orbital period of P_orb = 9 ± 0.5 days (Biddle et al. 2021) translates to a semimajor axis of a = 0.085 au. The system parameters utilized in our equations are presented in Table 1.

Massive companions are expected to carve out deep and wide gaps, which are identifiable as a lack of near-IR (NIR) emission in the system's spectral energy distribution (SED; Gonzalez et al. 2012). This is not the case for CI Tau as it has a substantial NIR excess (McClure et al. 2013). However, Muley & Dong (2021) found that a gap carved by a companion like CI Tau b can potentially be indistinguishable from an undepleted disk—at least at the NIR end of the SED. The planet's existence was directly questioned by Donati et al. (2020), who suggested that the signatures that have been attributed to the planet can instead be replicated by stellar activity. However, Biddle et al. (2018) found separate signatures for the planet's orbital period and the star's rotational period. They argued that both signatures cannot be fully attributed to the star alone.

Massive companions may, to some degree, mediate stellar accretion (Fouchet et al. 2010). Teyssandier & Lai (2020) found that such modulation can even be driven to the point of matching the companion's orbital period (see also Artymowicz & Lubow 1996; Muñoz & Lai 2016). Such periodic accretion has been observed in binary systems before, e.g., in photometric monitoring of the young stars DQ Tau (Tofflemire et al. 2017a) and TWA 3A (Tofflemire et al. 2017b). We seek to capture similar behavior in the Pfβ line, which has been calibrated against accretion luminosity (Salyk et al. 2013), originating from CI Tau.

We utilize high-resolution spectroscopy for detailed analysis of the emission line profiles originating from CI Tau to characterize its inner protoplanetary disk. Specifically, we study the rovibrational NIR emission of ¹²CO, which classical TTS systems like CI Tau commonly emit (Najita et al. 2003). The distribution of the emitting gas can be determined from the line profiles given a known stellar mass and disk inclination. Asymmetric features can allude to the presence of disk substructures such as disk winds (Pontoppidan et al. 2011), circumplanetary disks (CPDs; Brittain et al. 2019), and disk eccentricities (Liskowsky et al. 2012). Also, the disk's atmospheric temperature can be estimated from the relative transition strengths.

The following section (Section 2) provides a brief overview of the data collection routine and reduction process. Section 3 describes the asymmetries in the emission line profiles, and the variabilities are discussed in Section 3.1. CI Tau's average 2022 profile is analyzed and compared to simulated models in Section 3.2, first fitted to a circular disk model (Section 3.2.1) and then to a two-component disk model (Section 3.2.2). Results from the profile fitting are compared to a hydrodynamic simulation in Section 3.3. Section 4 addresses lingering questions and how our results contribute to the debate around CI Tau's hot Jupiter companion. Lastly, Section 5 highlights our main results.

The spectra were collected with the iSHELL cross-dispersion echelle spectrograph at the NASA Infrared Telescope Facility (IRTF, Rayner et al. 2022) for 11 of the 12 epochs in this study. Two of the epochs, from 2018 October 3 and 2019 January 3, were first included in Banzatti et al. (2022), and the rest, obtained during the period of 2022 January 21–29, are presented in this work for the first time. In addition, archival data were downloaded for one epoch, 2008 October 10, from an older survey with NIRSPEC (McLean et al. 1998) at the W. M. Keck Observatory. See Table 2 for an observation log.

In 2022, we observed CI Tau for nine consecutive nights to cover one full orbital period of the companion CI Tau b. The spectra were acquired with an ABBA nodding pattern that allowed for the images to be combined in an A–B–B+A pattern. This removes sky emission to first order. The position angle of the slit on the sky was along the semimajor axis of the system's outer disk. For CI Tau we took 5 s exposures with 12 coadds, while the telluric standard was observed with 5 s exposures with 5 coadds. For the iSHELL observations, the slit width of 0375 provided a spectral resolution of ∼92,000 (Banzatti et al. 2022). The average during the 2022 observation run was 084. The NIRSPEC observations were acquired with a 0432 slit width that provided a spectral resolution of ∼25,000.

The iSHELL data were reduced using the SpeXtool5 pipeline (Vacca et al. 2003; Cushing et al. 2004), with which the calibration frames (sky frame and master flat) are prepared, and the spectral orders are then straightened, extracted, and wavelength-calibrated using a sky model. The orders are then stitched together to construct the full spectrum. The telluric corrections are carried out by ratioing the wavelength-calibrated science and standard stars. This cancels out the sky as well as stellar contributions. Utilizing a model generated by the Sky Synthesis Program (Kunde & Maguire 1974), the same process is followed to reduce the NIRSPEC data with a custom routine described in Brittain et al. (2003).

The average M-band spectrum from the 2022 observation run is presented in Figure 1. It covers the ¹²CO v = 1–0 transitions from R(15)–P(42) as well as the hydrogen Pfβ transition. The low-J lines where J ≤ 20 (J being the rotational quantum number of the lower state) and ${\tilde{\nu }}_{J}\geqslant 2060\,{\mathrm{cm}}^{-1}$ (where ${\tilde{\nu }}_{J}$ is the central wavenumber of a transition) have a narrow absorption component superimposed on the broad emission features. This is an indication of foreground CO at a lower temperature than the emitting gas of the same species.

Zoom In Zoom Out Reset image size

The average emission line profiles of the ¹²CO low-J (J ≤ 20) lines, high-J (J > 20) lines, and the hydrogen Pfβ line are presented in Figure 2. The epochs start at the bottom and progress upwards (the 2022 epochs begin with #4) along with the 2022 average (red) plotted over each for comparison. To construct the average profile, each individual profile is adjusted so that their centers are aligned with each other in Doppler-shifted velocity space (typically less than 0.5 km s⁻¹). This is done because the wavelength calibrations may not be the same on different nights. The profiles are then shifted again to correct for the Doppler motion on their respective dates to center the profiles at V = 0 km s⁻¹. Finally, the average profile is calculated while excluding the profiles that have significant corruption.

Zoom In Zoom Out Reset image size

The average stacked line profiles for the ¹²CO low-J lines are presented in the left column of Figure 2. They possess an absorption component and a telluric contribution that we remove during the reduction pipeline—this is seen as a gap in the profiles. Because of these external features the majority of the profile analysis is done on the high-J lines (center of Figure 2).

Emission from a CPD would appear in the profiles as an asymmetric feature that oscillates over an orbital period with Keplerian velocities. The high-J line profiles (center of Figure 2) displayed no clear asymmetric feature that varied over the nine days. The lack of a variable feature places a limit on the emission from a CPD associated with CI Tau b and can be used to constrain the CPD radius.

To estimate the lower limit of our ability to detect CPDs, we estimate the minimum radius (${R}_{\mathrm{CPD}}^{\mathrm{limit}}$) a CPD must have to be detected in our data. A best-fit temperature profile of the inner disk (Section 3.2.1) indicates the temperature at CI Tau b's location is about 2300 K. This temperature is assumed to correspond to that of the CPD for simplicity. Using the wavenumber of the center of the M band (2127.7 cm⁻¹), the CPD's flux density is ${B}_{J}({\tilde{\nu }}_{J},T)=4.1\times {10}^{4}$ erg s⁻¹ cm⁻² sr⁻¹ cm, where ${B}_{J}({\tilde{\nu }}_{J},T)$ is the Planck function.

The intensity of the CPD, I_CPD, and B_J (${\tilde{\nu }}_{J}$, T) are related by ${I}_{\mathrm{CPD}}={\rm{\Omega }}{B}_{J}({\tilde{\nu }}_{J},T)$, where Ω is the solid angle ${\rm{\Omega }}=\pi {({R}_{\mathrm{CPD}}/d)}^{2}$ and d is the distance to the system (see Table 1). We set I_CPD to our detection limit of 5σ_F, where σ_F = 4.1 × 10⁻¹⁵ erg s⁻¹ cm⁻² cm is the standard deviation of the spectrum noise. Finally, ${R}_{\mathrm{CPD}}^{\mathrm{limit}}$ can be estimated as

$$\begin{eqnarray}&&{R}_{\mathrm{CPD}}^{\mathrm{limit}}\approx d\sqrt{\displaystyle \frac{5{\sigma }_{F}}{\pi {B}_{J}({\tilde{\nu }}_{J},T)}},\end{eqnarray} \tag{ 1 }$$

(see also Pineda et al. 2019). We find ${R}_{\mathrm{CPD}}^{\mathrm{limit}}=0.013$ au, or about 27 Jupiter radii.

A CPD encompasses a fraction of its Hill radius (R_H) and may be as large as 0.5R_H (Machida et al. 2008). Using the system parameters from Table 1, and a = 0.085 au, we obtain

$$\begin{eqnarray}&&{R}_{{\rm{CPD}}}\approx 0.5{R}_{{\rm{H}}}=0.5a(1-e){\left(\displaystyle \frac{{M}_{{\rm{P}}}}{3{M}_{\ast }}\right)}^{1/3}\approx 0.005\,{\rm{au}}.\end{eqnarray} \tag{ 2 }$$

Since this value falls below our detection limit ${R}_{\mathrm{CPD}}^{\mathrm{limit}}$, we are not able to confirm whether the CPD around CI Tau b is present or not.

Although short-term variability associated with a CPD was not observed in the 2022 epochs (#4–12), there is a clear overall difference when compared to the earlier epochs (#1–3). This is further discussed in Section 4.

The hydrogen Pfβ profiles are presented in the right column of Figure 2 with the 2022 average overplotted. Since this is an individual line at 2148.7 cm⁻¹, no averaging was performed. Because hydrogen probes regions much closer to the star the maximum velocities of these regions are much higher than those observed in the CO. These profiles also do not exhibit the double-peaked structure that is often seen from disks. This transition suffers from a degree of blending because it occurs between the ¹²CO R(0) and R(1) lines. Because of this we limit the hydrogen profiles to velocities from V = −250 km s⁻¹ to +200 km s⁻¹.

Similar to the CO profiles, the earlier epochs of the hydrogen profiles differ from the 2022 ones. Unlike CO, though, they do display some variability through the consecutive 2022 epochs. Since hydrogen elucidates stellar behavior, this may be an indication of the companion modulating stellar accretion rates (see Section 3.1).

3.1. Variability

The equivalent widths of the average line profiles of the CO transitions are presented in Table 3 and plotted as red points in Figure 3. The 2022 average value is 21.8 km s⁻¹ and is represented as a horizontal red dashed line. During the 2022 observation run, fluctuations about the average value are generally consistent with the measurement uncertainties of about 0.4 km s⁻¹; as a result, we have not detected any significant variability over the nine-day period. We do note that epoch #7 has a notably lower value of 19.46 km s⁻¹ while epoch #11 has a higher value of 24.30 km s⁻¹. We consider them outliers but they may be physically meaningful if future observations find similar patterns. If there is a pattern over the course of an orbit then it is not clear.

Zoom In Zoom Out Reset image size

Table 3. Equivalent Widths and Accretion Rates

Epoch	EW(CO)	EW(Pfβ)	log($\dot{M}/{M}_{\odot }\,{{\rm{yr}}}^{-1}$)
#	(km s−1)	(cm−1)
1	26.22 ± 0.22	0.29 ± 0.14	−7.76 ± 0.19
2	29.96 ± 0.30	0.33 ± 0.07	−7.71 ± 0.08
3	27.56 ± 0.19	0.43 ± 0.04	−7.61 ± 0.04
4	21.80 ± 0.75	0.71 ± 0.05	−7.41 ± 0.03
5	20.75 ± 0.61	0.60 ± 0.06	−7.47 ± 0.04
6	20.78 ± 0.27	0.54 ± 0.01	−7.52 ± 0.01
7	19.46 ± 0.42	0.55 ± 0.02	−7.51 ± 0.01
8	21.11 ± 0.56	0.46 ± 0.02	−7.58 ± 0.02
9	22.41 ± 0.39	0.46 ± 0.02	−7.58 ± 0.02
10	21.57 ± 0.35	0.51 ± 0.01	−7.54 ± 0.01
11	24.30 ± 0.26	0.53 ± 0.01	−7.52 ± 0.01
12	21.45 ± 0.31	0.62 ± 0.02	−7.46 ± 0.01

Download table as:  ASCIITypeset image

Overall, the 2022 measurements (#4–12) are markedly lower than those seen in earlier epochs (#1–3). This indicates changes in the disk over timescales much longer than the nine-day period of the 2022 run. This is discussed further in Section 4.

Also tabulated in Table 3 are the equivalent widths of the hydrogen Pfβ line alongside derived accretion rates (described later in this section). The equivalent widths are plotted in Figure 3 as blue points with the same units as CO. The hydrogen equivalent widths have a 2022 average of 73.7 km s⁻¹ and it is plotted as a horizontal dashed blue line. The equivalent widths do display a noticeable amount of variation over the nine-day period. Starting with the higher value at epoch #4 with 99.2 km s⁻¹, the equivalent width then gradually dips to 64.3 km s⁻¹ by epochs #8/#9 and then rises again to 86.6 km s⁻¹ at epoch #12. The small uncertainties in our measurements lend confidence to this trend being physical.

The 2022 measurements are higher than the earlier measurements in general. Because the 2022 measurements are themselves variable and the earlier epochs have large uncertainties, it is inconclusive whether the prior epochs truly differ from the 2022 values. This is discussed further in Section 4.

Since we lack an absolute flux calibration, it remains uncertain during the 2022 observation run how much variability is attributable to the star itself. The stellar variability over the course of a day can be gauged by studying previous surveys. The AllWISE Multiepoch Photometry Database provides 4.6 μm photometry with 21 hr of coverage on 2010 February 24 and 24 hr on 2010 September 4. The variation during those nights averaged 7% and that is smaller than what is captured in the Pfβ line. Thus, fluctuations caused by the star are likely insufficient to explain the variability we observe in the Pfβ line. Nonetheless, a photometrically calibrated spectroscopic study would help clarify the nature of the observed variability.

The luminosity of the Pfβ line has been correlated to UV-excess-derived accretion luminosities and thus can be used as a reliable tracer of stellar accretion rates (Salyk et al. 2013). The hydrogen equivalent widths are converted to their corresponding line luminosities utilizing the ALLWISE photometric data available and the Gaia DR3 distance. We do not include any extinction correction because it is expected to be negligible. While an M-band extinction coefficient, A _M , for CI Tau is not available, one may estimate it using the J-band extinction coefficient, A_J , of 0.51 reported by Kenyon & Hartmann (1990). Assuming grains to have properties like those in the interstellar medium, we have A _M /A_J = 0.095 (Mathis 1990), which translates to A _M = 0.048, or a flux correction of only about 4%.

The Pfβ line luminosity (L_Pfβ ) is converted to an accretion luminosity (L_acc) using the relation given by Salyk et al. (2013):

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{acc}}}{{L}_{\odot }}\right)=0.91\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{Pf}\beta }}{{L}_{\odot }}\right)+3.29,\end{eqnarray} \tag{ 3 }$$

where we have elected to drop the uncertainties on the numerical constants as they are largely systematic and only serve to shift all data points equally. It is more meaningful to look at the uncertainties due to the measured equivalent widths alone.

The accretion luminosities are converted to stellar accretion rates ($\dot{M}$) by equating L_acc to the release of gravitational energy from R_cor, the corotation radius in the disk, to R_*, the host star's radius:

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{R}_{* }{L}_{\mathrm{acc}}}{{{GM}}_{* }}{(1-{R}_{* }/{R}_{\mathrm{cor}})}^{-1}.\end{eqnarray} \tag{ 4 }$$

The factor ${(1-{R}_{* }/{R}_{\mathrm{cor}})}^{-1}$ reduces to a factor of 1.25 when R_cor is assumed to be 5R_* (Hartmann et al. 1998). We find an average accretion rate of $\dot{M}$ = (3.1 ± 2.2) × 10⁻⁸ M_☉ yr⁻¹, which agrees with the previously reported value of $\dot{M}$ = (2.5 ± 1.8) × 10⁻⁸ M_☉ yr⁻¹ (Donati et al. 2020).

The accretion rate varied by 35% over the nine-day orbital period. Comparing to theoretical hydrodynamic models by Teyssandier & Lai (2020), this level of variability can be caused by a companion with a mass of 9.4 M_J and an orbital eccentricity of about 0.05. Although that eccentricity is lower than that reported for CI Tau (Flagg et al. 2019), it is not beyond the realm of possibility given the large uncertainty in the reported eccentricity. Also, one should bear in mind that Teyssandier & Lai (2020) studied the disk accretion rate, which may have quantitative differences from the stellar accretion rate. While the companion may drive strong modulation in the flow of gas near its orbit, for that gas to travel further inward and ultimately accrete onto the star, it must still be subjected to some disk transport mechanisms. Transport by turbulence, for example, is diffusive and would likely smooth out the accretion flow. Detailed modeling (e.g., Fouchet et al. 2010) is needed to verify whether the modulation we have detected here is consistent with the measured parameters of CI Tau b.

We estimate an upper limit to CI Tau b's accretion rate (${\dot{M}}_{{\rm{p}}}^{\mathrm{limit}}$) by applying Equation (4) where the factor ${(1-{R}_{* }/{R}_{\mathrm{cor}})}^{-1}$ is now 1 by assuming, for simplicity, that the gas falls onto the planet from a distance much larger than the planet's radius. We can write

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{p}}}^{\mathrm{limit}}=\displaystyle \frac{{R}_{{\rm{p}}}}{{{GM}}_{{\rm{p}}}}{L}_{\mathrm{acc}}^{\mathrm{limit}},\end{eqnarray} \tag{ 5 }$$

where we choose R_p to be 2 R_J. According to Spiegel & Burrows (2012, 2013), a 2 Myr old giant companion with a mass of 10 M_J will have a radius of 1.94 R_J or 1.14 R_J depending on how it formed. Since CI Tau b has a slightly larger mass we used a value of 2 R_J for simplicity.

${L}_{\mathrm{acc}}^{\mathrm{limit}}$ in Equation (5) is related to our 5σ_F detection limit (see the sixth paragraph in Section 3), or

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}^{\mathrm{limit}}\approx 5{\sigma }_{{\rm{F}}}4\pi {d}^{2}{\rm{\Delta }}V,\end{eqnarray} \tag{ 6 }$$

where ΔV is the width of the planet's hydrogen line. Assuming that the speed of the planet's magnetospheric accretion flow is of the order of the freefall speed $\sqrt{{{GM}}_{{\rm{p}}}/{R}_{{\rm{p}}}}\approx 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and the line width is roughly two times that, we get ΔV ≈ 200 km s⁻¹, or about 1.3 cm⁻¹ in wavenumber units. Finally, plugging these values into Equations (5) and (6), the upper limit for CI Tau b's accretion rate is about 1 × 10⁻¹¹ M_☉ yr⁻¹, or 1 × 10⁻⁸ M_J yr⁻¹.

3.2. Emission Line Analysis

3.2.1. Circular Disk Model

We begin our analysis of the emission lines by assuming a circular disk. The spectrum (Figure 1) is fitted with a simple two-dimensional slab model that treats the emission as if it arises from a geometrically thin disk. A radial grid is defined from r = [r_in, r_out], which translates to the Keplerian velocity

$$\begin{eqnarray}&&{V}_{{\rm{K}}}(r)=\sqrt{\displaystyle \frac{{{GM}}_{* }}{r}}.\end{eqnarray} \tag{ 7 }$$

We choose to make the step sizes in the radial grid correspond with a step size of 1 km s⁻¹ in Keplerian velocities. In this way the grid has a finer resolution closer to the star, where the disk is more luminous. At each radius, the annulus is further divided by an angular array bounded from θ = [0, 2π]. θ is defined as the angle from the line of sight that has been projected onto the plane of the disk. Similarly, we made the angular step sizes correspond to a change in the projected velocities of 1 km s⁻¹. The projected velocities are

$$\begin{eqnarray}&&{V}_{{\rm{p}}}(r,\theta )={V}_{{\rm{K}}}(r)\sin \theta \sin i,\end{eqnarray} \tag{ 8 }$$

where i is the system's inclination. The projected velocities are then made the center of normalized Gaussian profiles, which are represented as

$$\begin{eqnarray}&&G(V-{V}_{{\rm{p}}})=\displaystyle \frac{1}{b\sqrt{\pi }}\exp \left(\displaystyle \frac{-{\left(V-{V}_{P}\right)}^{2}}{{b}^{2}}\right).\end{eqnarray} \tag{ 9 }$$

Here the Gaussian line width, b, is given the value 8.8 km s⁻¹ to emulate a realistic level of blending between discrete projected velocities.

Each annulus corresponds to a radially dependent temperature and surface number density. Assuming basic power laws, these are expressed as

$$\begin{eqnarray}&&T(r)={T}_{0}{\left(\displaystyle \frac{r}{{r}_{{\rm{in}}}}\right)}^{\alpha }\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&N(r)={N}_{0}{\left(\displaystyle \frac{r}{{r}_{\mathrm{in}}}\right)}^{\beta }.\end{eqnarray} \tag{ 11 }$$

Assuming thermodynamic equilibrium, the populations of the rovibrational states are

$$\begin{eqnarray}&&{N}_{J}=\displaystyle \frac{{{Ng}}_{J}}{Q(T)}\exp \left(\displaystyle \frac{-{E}_{J}}{{k}_{{\rm{B}}}{\text{}}T}\right),\end{eqnarray} \tag{ 12 }$$

where E_J and g_J are the energy and degeneracy of state J, respectively. Q(T) is the partition function as a function of temperature.

The optical depths, τ_J , are

$$\begin{eqnarray}&&{\tau }_{J}=\displaystyle \frac{{N}_{J}}{8b{\pi }^{3/2}}\displaystyle \frac{{g}_{J}{A}_{J}}{{g}_{J}^{{\rm{{\prime} }}}{\tilde{\nu }}_{J}^{3}},\end{eqnarray} \tag{ 13 }$$

where A_J and ${\tilde{\nu }}_{J}$ are the Einstein A coefficient and the central wavenumber of the transition. ${g}_{J}^{{\prime} }$ is the degeneracy of the adjacent state.

The flux densities are

$$\begin{eqnarray}&&{F}_{J}=(1-{e}^{-{\tau }_{J}})\,\displaystyle \frac{{\tilde{\nu }}_{J}}{c}\,{B}_{J}({\tilde{\nu }}_{J},T),\end{eqnarray} \tag{ 14 }$$

where c is the speed of light in a vacuum, and ${B}_{J}({\tilde{\nu }}_{J},T)$ is the Planck function expressed in units of wavenumbers (cm⁻¹). The Planck function is evaluated at the central wavenumber of each transition, ${\tilde{\nu }}_{J}$, and at the local temperature, T. The factor ${\tilde{\nu }}_{J}/{\rm{c}}$ converts the units to Doppler-shifted velocities (km s⁻¹).

Finally, the line profile, I_J (V), is obtained by convolving the emission from every grid point with G(V − V_p) and integrating over the full disk, or

$$\begin{eqnarray}&&{I}_{J}(V)={f}_{\mathrm{norm}}{\int }_{0}^{2\pi }{\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}{F}_{J}(r)\,G(V-{V}_{{\rm{p}}}(r,\theta ))\,r\,{dr}\,d\theta ,\end{eqnarray} \tag{ 15 }$$

where f_norm normalizes I_J (V) to our observed line profiles. The profiles of the high-J lines are then averaged together and compared to CI Tau's profile (top left panel of Figure 4). We perform a χ² optimization over six parameters: T₀ and α in the temperature profile (Equation (10)), N₀ and β in the surface density profile (Equation (11)), r_in, and r_out. Our best fit has a reduced χ² value of 3.3.

Zoom In Zoom Out Reset image size

The best-fit parameters are presented in Table 4. We find that the temperature and surface density at the companion's location (r = 0.085 au) are 2297 K and 2.9 × 10¹⁹ cm⁻², respectively. The radial extent of the disk goes from r_in = 0.048 ± 0.001 au to r_out = 0.898 ± 0.001 au.

Table 4. Best-fit Parameters

Model	Parameter	Best Fit
Circular	rin	0.048 ± 0.001 au
	rout	0.898 ± 0.001 au
	T0	2770 ± 150 K
	α	−0.33 ± 0.05
	N0	9.0 ± 2.0 × 1019 cm−2
	β	−2.0 ± 0.1
Multicomponent	rin	0.052 ± 0.002 au
	rout	1.53 ± 0.16 au
	rb	0.14 ± 0.08 au
	ωin	<40°
	ein	0.056 ± 0.015
	eout	0.048 ± 0.008
	Δ	0.009 ± 0.003 au
	ιout	3.11 ± 0.04

Note. See Sections 3.2.1 and 3.2.2 for descriptions of the circular disk model and the multicomponent disk model, respectively.

Download table as:  ASCIITypeset image

In the top left panel of Figure 4 we see that the circular disk model does broadly describe the observed line structures. However, when comparing just the red side of the profiles (V > 0 km s⁻¹), we see velocities where the model overshoots the data (V ∼ +50 km s⁻¹) but then undershoots it (V ∼ +80 km s⁻¹). This is peculiar because it cannot be matched even if we assign an eccentricity to the disk. Applying an eccentricity shifts the profile in one direction, effectively raising/lowering the entirety of one side to create a typical asymmetric line profile, but it cannot simultaneously increase and decrease the fluxes at different velocities on the same side. This discrepancy motivates us to propose a more complex disk model in the following section.

3.2.2. Multicomponent Disk Model

CI Tau's profile is better modeled when its disk is instead composed of multiple components each having their own eccentricities and arguments of periapsis. This is motivated by the fact that the line core appears blueshifted while the wings are redshifted (Figure 4). If CI Tau's disk is divided into inner and outer components, the outer component can account for the blueshifted core while the inner one accounts for the redshifted wings of CI Tau's profile.

In this two-component model, we introduce a break radius, r_b, that separates the inner and outer components. Within this radius the inner component will have an argument of periapsis ω_in that is different from that of the outer component ω_out. In our model, ω_in is defined to be the angle between the inner disk's semimajor axis and our line of sight projected onto the plane of the disk. The disks are required to be oppositely oriented to some degree in order to replicate CI Tau's profile (Figure 4) so, to reduce the number of free parameters, we chose to fix them to be antiparallel (ω_out = ω_in + π). The validity of this assumption is discussed in Section 3.3.

The two components are also given their own eccentricities, e_in and e_out, which are treated as constants throughout their respective component. We assign a width parameter Δ around the break radius where neither disk will be contributing flux. This may capture any potential disk gap that might be present.

Unlike in Section 3.2.1 where we fitted the overall spectrum with a temperature and surface density profile (Equations (10) and (11)) as a means to compute the disk's intensity (Equation (15)), here we directly fit the intensity as a function of radius I(r). Since we are now fitting only to the high-J lines there is not enough information to be able to constrain the temperature and surface density. I(r) is fitted as a broken power law that follows one exponent for the inner component ι_in and another one for the outer component ι_out. The broken power law is defined as follows:

$$\begin{eqnarray}&&I(r)={I}_{0}{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{-{\iota }_{{\rm{in}}}}{\left(\displaystyle \frac{1}{2}\left(1+{\left(\displaystyle \frac{r}{{r}_{{\rm{b}}}}\right)}^{1/\delta }\right)\right)}^{({\iota }_{{\rm{in}}}-{\iota }_{{\rm{out}}})\delta },\end{eqnarray} \tag{ 16 }$$

where I₀ is a normalization constant. We fix ι_in = 0 because the inner disk is likely narrow and it helps simplify our fit. δ determines the "smoothness" of transition between power laws. We set it to be $\delta =\mathrm{log}(1+{\rm{\Delta }}/{r}_{{\rm{b}}})$ in order to emulate a smooth transition that occurs over the distance Δ. In the intensity profile, ι_out is the only parameter that we vary.

Since our model now considers eccentricity, it is simpler to operate in semimajor axis space, where the distance to the star is now dependent on the semimajor axis, eccentricity, and azimuthal angle:

$$\begin{eqnarray}&&r=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos (\theta -\omega )},\end{eqnarray} \tag{ 17 }$$

where the variables e and ω are e_in/e_out and ω_in/ω_out for the inner/outer components, respectively. As with the circular disk model, θ is the angle between the grid element and the line of sight in the disk's plane.

The projected velocities also change depending on the eccentricity and argument of periapsis of the annulus. The new expression for the projected velocities is

$$\begin{eqnarray}&&{V}_{{\rm{P}}}(a,\theta )=\sqrt{\displaystyle \frac{{{GM}}_{* }}{a(1-{e}^{2})}}(\sin \theta +e\sin \omega )\sin i.\end{eqnarray} \tag{ 18 }$$

I(r, θ) is calculated for every grid element and transformed into a line profile in velocity space by multiplying it with the Doppler-shifted normalized Gaussian profile described in Equation (9). This is done over the extent of the disk, summed together, and normalized to CI Tau's profile by adjusting I₀ in Equation (16).

Our χ² optimization finds a best fit that has a reduced χ² of 0.984. This is a significant improvement over the circular disk model, which has a reduced χ² of 3.3 (Section 3.2.1). The best-fit profile of the multicomponent disk is displayed in the bottom left panel of Figure 4 together with the contributions from the inner and outer disks.

We confirm that the multicomponent model can replicate the blueshifted line core and the redshifted wings seen in CI Tau's profile, which is not possible for a circular disk. The "bump" at V = 80 km s⁻¹ is an exaggerated by-product of our model because of the lack of flux from the region centered on the discontinuity. A more detailed analysis can be performed in the future where this is filled with some realistic emission that should aid in smoothing out the "bump" and provide an even better fit.

The best-fit parameters for the multicomponent disk are presented in Table 4. The CO extends from r_in = 0.052 ± 0.002 au to r_out = 1.53 ± 0.16 au. Compared to the circular disk (Section 3.2.1), the inner radius marginally changed while the outer radius increased by 70%. The inner radius is also close to the truncation radius estimated by GRAVITY Collaboration et al. (2023; 0.034 ± 0.014 au) and may extend to within the star's corotation radius

$$\begin{eqnarray}&&{R}_{{\rm{rot}}}={\left(\displaystyle \frac{{{GM}}_{\ast }{P}_{{\rm{rot}}}^{2}}{4{\pi }^{2}}\right)}^{1/3}.\end{eqnarray} \tag{ 19 }$$

Assuming Keplerian rotation, we calculate the corotation radius to be 0.069 au. Similar results have been found for other T Tauri stars (e.g., Carr 2007). The same calculation for CI Tau's corotation radius was done by GRAVITY Collaboration et al. (2023) but there they assumed a different rotational period and stellar mass.

The best-fit values for the eccentricities of the inner and outer components are e_in = 0.056 ± 0.015 and e_out = 0.048 ± 0.008. These values can be seen as luminosity-weighted averages that are more representative of the inner edges of the components, which are more luminous. Overall, these eccentricities are much smaller than that of CI Tau b's orbit (e = 0.25; Flagg et al. 2019). This may be an indication of the planet's orbit being less eccentric, since one might expect the planet's and disk's eccentricities to be similar. This is discussed further in Section 3.3.

Our best fit gives Δ = 0.009 ± 0.003 au, which is the half-width of the disk break. As a reminder, it also represents a region in the disk (i.e., a "gap") where we assigned no emission. This is about 6% of r_b = 0.14 au, the location of the break, which implies that the transition occurs sharply. Taking into account the eccentricities of the disks, the full width of the gap represented by Δ varies between 0.0035 au at the apoapsis (periapsis) of the inner (outer) disk and 0.032 au at the periapsis (apoapsis) of the inner (outer) disk. Considering the uncertainties in our fit, this gap is too narrow for us to definitively confirm its existence. It mainly helps to separate the two disks just enough to avoid overlapping (see the right panel of Figure 4).

We find the argument of periapsis of the inner disk (ω_in) to have an upper limit at 40°. Again this describes the angle between the semimajor axis of the inner component and the line of sight projected onto the plane of the disk. This would correspond to the outer disk having a periapsis of ω_out = ω_in + π = 220°.

The outer component's intensity profile follows an ι_out = 3.11 dependence. This value is positive because a negative sign is incorporated in Equation (16). The large value of ι_out indicates a steep drop in flux—a sign that a large portion of the disk's emission originates from the inner component.

We would like to remind the reader that our results are sensitive to the mass of the star and the inclination of the disk. Here we assume an inclination of i = 71° (GRAVITY Collaboration et al. 2023) and a stellar mass of M_* = 1.02 M_☉ (Law et al. 2022) but these parameters have changed multiple times now. For example, the mass of the star was originally 0.8 M_☉ (Guilloteau et al. 2014) but then it was updated once to 0.9 M_☉ (Simon et al. 2019) and then again to the current value we used. It is possible that these values can be updated further and, if they are, then our best-fit radii would change by the following factor:

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{{\rm{new}}}}{{r}_{{71}^{\circ }}}={\left(\displaystyle \frac{\sin ({i}_{{\rm{new}}})}{\sin ({71}^{\circ })}\right)}^{2}\left(\displaystyle \frac{{M}_{{\rm{new}}}}{1.02\,{M}_{\odot }}\right),\end{eqnarray} \tag{ 20 }$$

where r_71° corresponds to the values assigned to the r_in, r_out, and r_b parameters. Changing these quantities by a constant factor will have no qualitative effect on our overall disk model.

3.3. Hydrodynamic Simulation

In this section, we present a proof-of-concept simulation that qualitatively demonstrates how a planet might be able to generate the disk features inferred by our model. The simulation is performed using PEnGUIn (Fung 2015) with a setup similar to the one used by Muley et al. (2019), but here the planet's mass (M_P) and orbit (a and e) are fixed. We pick a representative planet mass M_P = 1 M_J, and assign it an orbital semimajor axis that equals r_b in our multicomponent model (Table 4), and an orbital eccentricity of 0.05, similar to the eccentricities of both the inner and outer disks. We choose not to use a planet mass as large as the proposed mass of CI Tau b (11.6 M_J) because Kley & Dirksen (2006) had previously demonstrated that planets more massive than about 5 M_J would excite disk eccentricities much larger than those inferred by our model.

The simulation is locally isothermal and follows Equation (10) and Table 4 for the temperature profile. The surface density has an initial power-law profile the same as Equation (11) with β = −2, but the normalization is set to 1 in code units. Because we do not consider the self-gravity of the disk, the normalization to the surface density profile plays no role in the dynamics simulated. In Figure 5, we assign a physical value to the normalization based on the fact that we did not detect a planetary gap (Section 3.2.2). Also we assume a constant CO/H₂ ratio. The Sunyaev–Shakura viscosity parameter is assumed to be 0.01. Our choice describes a turbulent disk, which is plausible at short distances from the star where the ionization fraction is expected to be high and the magnetorotational instability (MRI; Balbus & Hawley 1991) is expected to be active; the choice of 0.01 is roughly similar to that for the most turbulent MRI-active disks (e.g., Guan & Gammie 2009; Simon & Hawley 2009).

Zoom In Zoom Out Reset image size

The final snapshot of the hydrodynamic simulation after 1000 planetary orbits is presented in Figure 5. The top left and top right panels are in Cartesian and polar coordinate systems, respectively. The inner and outer disks have developed a subtle eccentricity that generally appears similar in magnitude to the planet's eccentricity, but only near the planet's orbit. A slightly higher planet mass and/or a higher planetary eccentricity might better match our observation. However, the disk would develop an eccentricity that is too high if the planet were given the same parameters as CI Tau b.

Although the arguments of periapsis of the inner and outer disks do not appear to be exactly antiparallel as they were in our model (Figure 4), they are far from aligned. The precise geometry of the inner and outer disks is complicated. While parts of the disk further away from the planet's orbit may be approximated by eccentric orbits (Kley & Dirksen 2006), gas streamlines that co-orbit with the planet should follow horseshoe orbits that are deformed by the planet's eccentricity (Pan & Sari 2004). Even though our slab model does not capture all the nuances, we are confident that it does capture the fact that the inner and outer disks have distinct geometries.

In this simulation, the planet carves out a gap that is depleted in gas by about two orders order of magnitude and about 0.05 au in width. If the disk is sufficiently dense, the CO emitting layer can remain undepleted, producing the optically thick emission we observed, even inside the gap; but what defines "sufficiently" dense is a complicated matter that involves not only the H/CO ratio, but also the expected abundances of CO at low column densities (e.g., Bruderer 2013; Doppmann et al. 2017). To estimate how much CO is needed to produce optically thick emissions, we first set τ_J = 1 in Equation (13). From there, we obtain the rovibrational population N_J , which is then converted to the CO surface density, N, using Equation (12). In the bottom row of Figure 5 we plot the τ_J = 1 region ("optically thick boundaries") for the high-J P(20)–P(40) lines. The azimuthally averaged column density profile of the simulation (orange in Figure 5) is scaled so that it lies just above the τ_J = 1 regime and the detected CO surface density profile (green). This represents a possible version of the disk where a planet has carved out a gap that our observation is not sensitive to, though it is possible for the disk to have an even higher density.

Taking this further, we can see what constraint this assumption implies for the disk mass. In the bottom of Figure 5, we see that the unperturbed profile (orange dashed) of the simulated model needs to be about two orders of magnitude higher than Equation (11), the column density of the emitting layer. Or, more precisely, the full column density may be expressed as ${\rm{\Sigma }}=9.0\times {10}^{21}{(r/{r}_{\mathrm{in}})}^{-2}\,{\mathrm{cm}}^{-2}$. The total number of CO molecules can then be calculated by integrating this over the extent of the disk. Or,

$$\begin{eqnarray}&&{N}_{\mathrm{CO}}=2\pi {\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}9.0\times {10}^{21}\,{\mathrm{cm}}^{-2}{\left(\displaystyle \frac{r}{{r}_{\mathrm{in}}}\right)}^{-2}r\,{dr}.\end{eqnarray} \tag{ 21 }$$

Using the parameters from Table 4, we get N_CO ∼ 10⁴⁷ CO molecules. Assuming a CO/H₂ ratio of 1.6 × 10⁻⁴ (France et al. 2014), this translates to a mass of about 10⁻⁶ M_☉ or 10⁻³ M_J. This low disk mass may help explain why the planet can remain on an eccentric orbit over many orbital periods—eccentricity damping due to planet–disk interaction scales linearly with disk gas density (Artymowicz 1992). Future detailed modeling can better evaluate the strength of eccentricity damping and explore the implications for the planet's orbital evolution. Because this estimate of disk mass uses a density profile that lies just above our detected CO column, it can be seen as a lower limit.

The hot Jupiter companion CI Tau b was discovered by Johns-Krull et al. (2016) via IR radial velocity monitoring, and Flagg et al. (2019) detected CO directly from the planet's atmosphere. These studies cross-validated the planet's orbital period to be nine days. The host star's rotational period was found by Biddle et al. (2018) to be 6.6 days through analysis of the K2 light curves. However, the origin of these signatures was questioned by Donati et al. (2020), who claimed that both signatures can be attributed to the star itself. A follow-up analysis of the K2 photometry by Biddle et al. (2021) affirms that signature strengths cannot be replicated by stellar activity alone. Whether CI Tau b as a planet exists remains an on-going investigation.

One of our results that may validate CI Tau b is the seemingly nine-day period captured in the stellar accretion rate, derived from the hydrogen equivalent widths (Section 3.1). We also find some indirect evidence for a planet's existence (but not necessarily CI Tau b) through planet–disk interaction, where we find an inner and an outer disk that have separate eccentricities and arguments of periapsis (Section 3.2.2). Further investigations are needed to explain some other features we observed. Below we list a few that we find the most puzzling.

As illustrated in Figure 3, we observed the hydrogen Pfβ equivalent widths to vary by 35% during the 2022 epochs (#4–12; Section 3.1). This variation also encompasses the values observed in epoch #3 but not #1 and #2. Since it is ambiguous how representative our data are of the overall behavior, it remains unclear whether epochs #1 and #2 are atypical or not. For instance, they are around 45% from the 2022 mean and have relatively large uncertainties; statistically, it is possible that they are consistent. However, if they are meaningfully different then that may indicate changes in the star's luminosity and accretion rates. Classical TTSs are known to be quite variable over a wide array of time periods so CI Tau could have simply changed its luminosity or accretion rate between observation runs. Another possible explanation is that if the accretion flow is not isotropic, its direction might be related to the disk's eccentricity. Since eccentric disks driven by planet–disk interaction are known to precess over hundreds of planetary orbits (e.g., Kley & Dirksen 2006), the direction of the accretion flow may change with it and we would observe a variable accretion luminosity as the accretion flow goes in and out of alignment with our line of sight.

Unlike hydrogen, the ¹²CO equivalent widths (Figure 3) appear less variable (Section 3.1). In 2022 they varied by 12%, and that, along with having much smaller uncertainties, leads us to believe that the earlier epochs (#1–3), which are >20% from the 2022 mean, are different. Similar to the discussion in the previous paragraph, this may be due to the disk precessing and altering the amount of emission directed along our line of sight. Changes in the disk's eccentricity were observed over many planetary orbits (Section 3). However, there remains the possibility, again, that the star changed its luminosity and, as such, the temperature of the disk. If this were the case, then it means that from 2008 to 2022, the disk's luminosity decreased while the star's increased.

CI Tau's average emission line profile of 2022 was fitted with a disk that is constructed from two components (Figure 4) that have their own eccentricities (Section 3.2.2). Our derived eccentricities, which are about 0.05 (Table 4), are lower than that of CI Tau b, which is e = 0.25 ± 0.16. We also find that the disk has a very narrow gap if it exists; in other words, we find no evidence of gas at any radius traveling with an eccentricity as large as CI Tau b's. It is not a stable configuration for a planet to have a highly eccentric orbit and be simultaneously embedded in a near-circular disk—either the disk will damp the planet's eccentricity (e.g., Duffell & Chiang 2015) or the planet will force the disk to become more eccentric (e.g., Bitsch et al. 2013). Given the large uncertainty in CI Tau b's eccentricity, it may be possible that its true value is only about 0.05. According to the analysis by Donati et al. (2020), it is also possible that CI Tau b does not exist, and here we are instead observing the influence of another planet with a much lower eccentricity, such as the one simulated in Section 3.3.

CI Tau was observed for nine consecutive nights with the NASA IRTF in 2022 January (Section 2). These data were reduced (Figure 1) and paired with older data taken in 2008 using Keck and in 2018/2019 also with IRTF, thus giving us a total of 12 epochs (Table 2). For each epoch, we constructed the emission line profiles of the ¹²CO rovibrational transitions and the hydrogen Pfβ line (Figure 2). We fitted CI Tau's average high-J line profile from 2022 with flat, two-dimensional disk models, where we considered both a circular model (Section 3.2.1) and an eccentric model with two components (Section 3.2.2). From these fits, we constrained the disk properties (Table 4), such as the profiles of temperature (Equation (10)) and surface density (Equation (11)), and the eccentricities (e_in and e_out) of the different components. Our main conclusions are highlighted below.

• 1.  
  The core of the average stacked ¹²CO line profile is blueshifted while the wings are redshifted (Section 3.2.2). This is well modeled by introducing an inner and an outer disk component that are both eccentric but oppositely oriented (Figure 4). In our best-fit model, the components are separated near r_b = 0.14 au (Table 4). Both components also have their own eccentricities of about 0.05. The disk's structure may be explained by an embedded giant companion around r_b with a similar eccentricity (Section 3.3).
• 2.  
  A nine-day variability was observed in the hydrogen Pfβ line (Figure 3), which indicates variations in the stellar accretion rate (Section 3.1). A nine-day periodicity has been reported before by multiple groups, who attributed it to the presence of a giant companion (Johns-Krull et al. 2016; Biddle et al. 2021) or stellar activity (Donati et al.2020).
• 3.  
  The ¹²CO lines displayed varying asymmetries over hundreds of orbits (Figure 2). This might indicate that the inner disk around CI Tau is eccentric and precessing (Section 3).
• 4.  
  We did not detect ¹²CO emission from the circumplanetary disk around the proposed planet CI Tau b (Johns-Krull et al. 2016; Flagg et al. 2019; Biddle et al. 2021), and showed that we were likely not sensitive to it (Section 3). We also did not detect hydrogen Pfβ emission from the planet, which allowed us to place an upper limit on the planet's accretion rate of 1 × 10⁻⁸ Jupiter masses per year (Section 3.1).
• 5.  
  The inner radius of the protoplanetary disk was fitted to be 0.052 ± 0.002 au, which is consistent with the truncation radius of 0.034 ± 0.014 au reported by GRAVITY Collaboration et al. (2023) and may extend within the star's corotation radius of 0.069 au.
• 6.  
  The average accretion rate of CI Tau in 2022 is 3.1 × 10⁻⁸ M_☉ yr⁻¹. This aligns with what was previously reported by Donati et al. (2020; 2.5 × 10⁻⁸).

Being potentially the youngest host of a hot Jupiter companion, the CI Tau system can serve as the testing ground for theories pertaining to planet–disk interactions. Below, we chart a few directions in which to further investigate this system.

The CO and hydrogen equivalent widths (Table 3 and Figure 3) were observed to be markedly different in prior years. Future observations of the CI Tau system separated by a similar time span may be able to confirm and characterize this long-term behavior. Explaining it, such as whether it is caused by disk precession, would require detailed modeling.

Our analysis of the ¹²CO emissions (Section 3.2) utilized a simple two-dimensional slab model that does not adequately account for three-dimensional effects such as disk flaring and vertical temperature variations. Future analysis may incorporate a more rigorous disk geometry that may capture more complexities.

The hydrodynamic simulation in this study (Section 3.3) was a proof-of-concept for how an embedded giant planet may create the disk features observed. A more rigorous numerical study may better constrain the parameters of the planet–disk system.

We thank Jean-François Donati for a helpful discussion and providing input on the CI Tau system. This work includes data gathered at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.