Ro-vibrational Spectroscopy of CI Tau -- Evidence of a Multi-Component Eccentric Disk Induced by a Planet

CI Tau is currently the only T Tauri star with an inner protoplanetary disk that hosts a planet, CI Tau b, that has been detected by a radial velocity survey. This provides the unique opportunity to study disk features that were imprinted by that planet. We present multi-epoch spectroscopic data, taken with NASA IRTF in 2022, of the ${}^{12}$CO and hydrogen Pf$\beta$ line emissions spanning 9 consecutive nights, which is the proposed orbital period of CI Tau b. We find that the star's accretion rate varied according to that 9~d period, indicative of companion driven accretion. Analysis of the ${}^{12}$CO emission lines reveals that the disk can be described with an inner and outer component spanning orbital radii 0.05-0.13~au and 0.15-1.5~au, respectively. Both components have eccentricities of about 0.05 and arguments of periapses that are oppositely aligned. We present a proof-of-concept hydrodynamic simulation that shows a massive companion on a similarly eccentric orbit can recreate a similar disk structure. Our results allude to such a companion being located around an orbital distance of 0.14~au. However, this planet's orbital parameters may be inconsistent with those of CI Tau b whose high eccentricity is likely not compatible with the low disk eccentricities inferred by our model.


INTRODUCTION
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). Amongst 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 gasrich, 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 & c are two known protoplanets around a TTS 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. Amongst 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 & 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 planetdisk 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;Rosotti et al. 2017;Duffell & Chiang 2015;Ragusa et al. 2018;Muley et al. 2019a).
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. 2022) 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 d (Biddle et al. 2021) translates to a semi-major axis of a = 0.085 au. The system parameters utilized in our equations are presented in Table 1.
It is expected for massive companions to carve out deep and wide gaps which are identifiable as a lack of 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. (2018a) 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 ro-vibrational NIR emission of the 12 CO; which CTTS systems like CI Tau commonly emits (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 elude to the presence of disk sub-structures 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 (Sec. 2) provides a brief overview of the data collection routine and reduction process. Sec. 3 describes the asymmetries in the emission line profiles and the variabilities are discussed in Sec. 3.1. CI Tau's average 2022 profile is analyzed and compared to simulated models in Sec. 3.2; first fitted to a circular disk model (Sec. 3.2.1) and then a twocomponent disk model (Sec. 3.2.2). Results from the profile fitting are compared to a hydrodynamic simulation in Sec. 3.3. Sec. 4 addresses lingering questions and how our results contribute to the debate around CI Tau's hot Jupiter companion. Lastly, Sec. 5 highlights our main results.

OBSERVATIONS AND REDUCTIONS
The spectra were collected with the iSHELL crossdispersion 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 03 Oct 2018 and 03 Jan 2019, were first included in Banzatti et al. (2022) and the rest, obtained during the period of 21-29 Jan 2022, are presented in this work for the first time. In addition, archival data was downloaded for one epoch, 10 Oct 2008, 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 9 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 second exposures with 12 coadds while the telluric standard was observed with 5 second exposures with 5 coadds. For the iSHELL observations, the slit width of 0.375 ′′ provided a spectral resolution of ∼92,000 (Banzatti et al. 2022). The average during the 2022 observation run was 0.84 ′′ . The NIRSPEC observations were acquired with a 0.432 ′′ slit width that provided a spectral resolution of ∼25,000.
The iSHELL data were reduced using the SpeXtool5 pipeline (Cushing et al. 2004;Vacca et al. 2003) with which the calibration frames (sky frame and master flat) are prepared and the spectra 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 NIR-SPEC data with a custom routine described in Brittain et al. (2003).

RESULTS
The average M -band spectrum from the 2022 observation run is presented in Fig. 1. It covers the 12 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ν J ≥ 2060 cm −1 (whereν 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.
The average emission line profiles of the 12 CO low-J (J≤20) lines, high-J (J>20) lines and the hydrogen Pfβ line are presented in Fig. 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 then 0.5 km s −1 ). This is done because the wavelength calibrations may not be the same between nights. The profiles are then shifted again to correct for their respective dates' Doppler motion to center the profiles at V = 0 km s −1 . Finally, the average profile is calculated while excluding the profiles that have significant corruption.
The average stacked line profiles for the 12 CO low-J lines are presented in the left column of Fig. 2. They possess an absorption component and a telluric contribution that we remove during the reduction pipelinethis 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: Fig. 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; Fig. 2) displayed no clear asymmetric feature that varied over the 9 days. The lack of a variable feature places a limit on the emission from a CPD associ-ated 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 limit CPD ) a CPD must have to be detected in our data. A best-fit temperature profile of the inner disk (Sec. 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 M-band centered wavenumber (2127.7 cm −1 ), the CPD's flux density is The intensity of the CPD, I CPD , and B J (ν J ,T ) are related by I CPD = ΩB J (ν J , T ), where Ω is the solid angle Ω = π(R 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 −15 erg s −1 cm −2 cm is the standard deviation of the spectrum noise. Finally, R limit CPD can be estimated as (see also Pineda et al. (2019)). We find R limit CPD = 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 (2) Since this value falls below our detection limit R limit CPD , 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 Sec. 4.
The hydrogen Pfβ profiles are presented in the right column of Fig. 2 with the 2022 average overplotted. Since this is an individual line at 2148.7 cm −1 , no averaging was performed. Because hydrogen probes regions much closer to the star their maximum velocities are much higher. 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 12 CO R(0) and R(1) lines. Because of this we limit the hydrogen profiles to velocities from V = -250 km s −1 to +200 km s −1 .
Similar to the CO profiles, the earlier epochs of the hydrogen profiles differ from the 2022 ones. Unlike CO,

Variability
The equivalent widths of the average line profiles of the CO transitions are presented in Table 3 and plotted as red points in Fig. 3. The 2022 average value is 21.8 km s −1 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 −1 ; as a result, we have not detected any significant variability over the 9d period. We do note that epoch #7 has a notably lower value of 19.46 km s −1 while epoch #11 has a higher value of 24.30 km s −1 . 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.
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 9 d period of the 2022 run. This is discussed further in Sec. 4. Also tabulated in Table 3 are the equivalent widths of the hydrogen Pfβ line alongside derived accretion rates (described later in this Sec.). The equivalent widths are plotted in Fig. 3 as blue points with the same units as CO. The hydrogen equivalent widths have a 2022 average of 73.7 km s −1 and it is plotted as a horizontal dashed blue line. The equivalent widths do display a no-ticeable amount of variation over the 9 d period. Starting with the higher value at epoch #4 with 99.2 km s −1 , the equivalent width then gradually dips to 64.3 km s −1 by epochs #8/#9 and then rises again to 86.6 km s −1 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 Sec. 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 hrs of coverage on 24 Feb 2010 and 24 hrs on 4 Sept 2010. 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 ISM-like grain properties, 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): 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 (Ṁ ) by equating L acc to the release of gravitational energy from R cor , the co-rotation radius in the disk, to R * , the host star's radius: The (1 − R * /R 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Ṁ = 3.1±2.2 ×10 −8 M ⊙ yr −1 which agrees with the previously reported value ofṀ = 2.5±1.8 ×10 −8 M ⊙ yr −1 (Donati et al. 2020). The accretion rate varied by 35% over the 9 d 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 (Ṁ limit p ) by applying Eqn. 4 where the (1−R * /R cor ) −1 factor is now 1 by assuming, for simplicity, the gas falls onto the planet from a distance much larger than the planet's radius. We can writė where we choose R p to be 2 Jupiter radii. According to Spiegel & Burrows (2012, 2013, a 2 Myr 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 limit acc in Eqn. 5 is related to our 5σ F detection limit (see the sixth paragraph in Sec. 3), or where ∆V is the width of the planet's hydrogen line. Assuming the planet's magnetospheric accretion flow speed  Table 3) with the 2022 averages plotted as the horizontal dashed lines ( 12 COavg = 21.8 km s −1 & Havg =73.7 km s −1 ). In the 2022 epochs (#4-12) we see the hydrogen vary with a possible 9 d period whereas the 12 CO stayed consistent. The 2022 epochs are noticeably different from the earlier ones.
is of order the free fall speed GM p /R p ≈ 100 km s −1 and the line width is roughly two times that, we get ∆V ≈ 200 km s −1 or about 1.3 cm −1 in wavenumber units. Finally, plugging these values into Eqn. 5 and 6, the upper limit for CI Tau b's accretion rate is about 1 × 10 −11 M ⊙ yr −1 , or 1 × 10 −8 M J yr −1 .

Circular Disk Model
We begin our analysis of the emission lines by assuming a circular disk. The spectrum (Fig. 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 We choose to make the step sizes in the radial grid correspond with a step size of 1 km s −1 in Keplerian velocities. 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 −1 . The projected velocities are where i is the system's inclination. The projected velocities are then made the center of normalized Gaussian profiles which are represented as Here the Gaussian line width, b, is given the value 8.8 km s −1 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 and Assuming thermodynamic equilibrium, the populations of the ro-vibrational states are 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 where A J andν J are the Einstein A coefficient and central wavenumber of the transition. g ′ J is the degeneracy of the adjacent state.
The flux densities are where c is the speed of light in a vacuum, and B J (ν J , T ) is the Planck function expressed in units of wavenumbers (cm −1 ). The Planck function is evaluated at the central wavenumber of each transition,ν J , and at the local temperature, T . The factorν J /c converts the units to Doppler-shifted velocities (km s −1 ). Finally, the line profiles, 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 where f norm normalizes the 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 Fig. 4). We perform a χ 2 optimization over six parameters: T 0 and α in the temperature profile (Eqn. 10); N 0 and β in the surface density profile (Eqn. 11); r in ; and r out . Our best-fit has a reduced χ 2 value of 3.3.
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 19 cm −2 , respectively. The radial extent of the disk goes from r in = 0.048 ± 0.001 au to r out = 0.898 ± 0.001 au.
In the top left panel of Fig. 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 −1 ), we see velocities where the model overshoots the data (V ∼ +50 km s −1 ) but then undershoots it (V ∼ +80 km s −1 ). This is peculiar because it cannot be matched even if we assign an eccentricity to the disk. Applying an eccentricity shifts the profile towards 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.

Multi-component 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 periapses. This is motivated by the fact that the line core appears blue-shifted while the wings are red-shifted (Fig. 4). If CI Tau's disk is divided into inner and outer components, the outer component can account for the blue-shifted core while the inner component accounts for the red-shifted 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 periapse ω 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 semi-major 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 (Fig.4) so, to reduce the number of free parameters, we chose to fix them to be anti-parallel (ω out = ω in + π). The validity of this assumption is discussed in Sec. 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 Sec. 3.2.1 where we fitted the overall spectrum with a temperature and surface density profile (Eqns. 10 and 11) as a means to compute the disk's intensity (Eqn. 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: where I 0 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 δ = log(1 + ∆/r 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 semi-major axis space, where the distance to the star is now dependent on the semi-major axis, eccentricity and azimuthal angle: where the variables e and ω are e in /e out or ω in /ω out for the inner/outer components, respectively. Like 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 periapse of the annulus. The new expression for the projected velocities is V P (a, θ) = GM * a(1 − e 2 ) (sin θ + e sin ω) sin i .
An 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 Eqn. 9. This is done over the extent of the disk, summed together, and normalized to CI Tau's profile by adjusting I 0 in Eqn. 16.
Our χ 2 optimization finds a best-fit that has a reduced χ 2 of 0.984. This is a significant improvement over the circular disk model that has a reduced χ 2 of 3.3 (Sec. 3.2.1). The best-fit profile of the multi-component disk is displayed in the bottom left panel of Fig. 4 together with the contributions from the inner and outer disks.
We confirm that the multi-component model can replicate the blue-shifted line core and the red-shifted wings seen in CI Tau's profile; which is not possible for a circular disk. The "bump" at V = 80 km s −1 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 multi-component 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 (Sec. 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 co-rotation radius 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 co-rotation radius was done by GRAVITY Collaboration et al. (2023) but there they assumed a different rotational period and stellar mass.
The best-fit value 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 luminosityweighted 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 Sec. 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, to 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 Fig. 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 semi-major 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 a ι out = 3.11 dependence. These values are positive because a negative sign is incorporated in Eqn. 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 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 do, then our best-fit radii would change by the following factor 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.  Table 4 for the best-fit parameters). The bottom-left panel is the multi-component disk fit (Sec. 3.2.2) which allowed for the inner+outer components of the disk to contribute independently and provide a much better fit (see Table 4 for the best-fit parameters). The right figure is a visual of the best-fit multi-component disk where the inner+outer components are eccentric and anti-parallel.

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. (2019b), 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 semi-major axis that equals r b in our multi-component model (Table  4), and an orbital eccentricity of 0.05, similar to the eccentricities of both the inner and outer disk. We choose not to use a planet mass as large as the proposed CI Tau b's mass (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 Eqn. 10 and Table 4 for the temperature profile. The surface density has an initial power-law profile the same as Eqn. 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 Fig. 5, we assign a physical value to the normalization based on the fact that we did not detect a planetary gap (Sec. 3.2.2). Also we assume a constant CO/H 2 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 the most turbulent MRI-active disks (e.g., Simon & Hawley 2009;Guan & Gammie 2009).
The final snapshot of the hydrodynamic simulation after 1000 planetary orbits is presented in Fig. 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 bet- Figure 5. The top row is the final snapshot of the hydrodynamic simulation (Sec. 3.3) with a hypothetical 1MJ planet that may reproduce the disk eccentricities we observed (0.05; Fig. 4). The color scale is logarithmic and in units of CO molecules cm −2 . The left and right panels are the same picture in Cartesian and polar coordinates, respectively. The white dashed line traces the orbit of the simulated planet. In this case, the outer disk develops an eccentricity similar to the planet's, while the eccentricity of the inner disk is more subtle. The bottom row plots the azimuthally averaged CO column density profile. The orange solid curve is taken from the same simulation and time-averaged over the last orbit; the orange dashed curve represents the unperturbed disk; and the green solid curve is the emission we detected (Table 4; Sec. 3.2.1). A τ =1 region was calculated for the P(20)-P(40) transitions and plotted as the red region (Sec. 3.3). The simulated profile (orange) is normalized so that, even at the bottom of the planetary gap, it lies above our detected CO column, illustrating that even though the gap is depleted by a factor of ∼100, its 12 CO emission can remain optically thick (Sec. 3.3).
ter match our observation. However, the disk would develop an eccentricity that is too high if the planet was given the same parameters as CI Tau b.
Although the arguments of periapses of the inner and outer disks do not appear to be exactly anti-parallel as they were in our model (Fig. 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 Eqn. 13. From there, we obtain the ro-vibrational populations N J , which is then converted to the CO surface density, N , using Eqn. 12. In the bottom row of Fig. 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: Fig. 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 Fig. 5, we see that the unperturbed profile (orange dashed) of the simulated model needs to be about two orders of magnitude higher than Eqn. 11, the column density of the emitting layer. Or, more precisely, the full column density may be expressed as Σ = 9.0 × 10 21 (r/r in ) −2 cm −2 . The total number of CO molecules can then be calculated by integrating this over the extent of the disk. Or, N CO = 2π rout rin 9.0 × 10 21 cm −2 r r in −2 r dr . (21) Using the parameters from Table 4, we get N CO ∼ 10 47 CO molecules. Assuming a CO/H 2 ratio of 1.6×10 −4 (France et al. 2014), this translates to a mass of about 10 −6 M ⊙ or 10 −3 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 on the planet's orbital evolution. Because this disk mass estimate uses a density profile that lies just above our detected CO column, it can be seen as a lower limit.

DISCUSSION
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 9 d. The host star's rotational period was found by Biddle et al. (2018b) to be 6.6 d through analysis of the K2 lightcurves. 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 piece of our results that may validate CI Tau b is the seemingly 9 d period captured in the stellar accretion rate, derived from the hydrogen equivalent widths (Sec. 3.1). We also find some indirect evidence for a planet's existence (but not necessarily CI Tau b) through planetdisk interaction, where we find an inner and outer disk that have separate eccentricities and arguments of periapses (Sec. 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 Fig. 3, we observed the hydrogen Pfβ equivalent widths to vary by 35% during the 2022 epochs (#4-12; Sec. 3.1). This variation also encompasses the values observed in epoch #3 but not #1-2. Since it is ambiguous how representative our data is of the overall behavior, it remains unclear if epochs #1-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. CTTSs 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 light-of-sight.
Unlike hydrogen, the 12 CO equivalent widths (Fig. 3) appear less variable (Sec. 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 dif-ferent. 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-ofsight. Changes in the disk's eccentricity was observed over many planetary orbits (Sec. 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 (Fig. 4) that have their own eccentricities (Sec 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 Sec. 3.3.

SUMMARY AND CONCLUSIONS
CI Tau was observed for 9 consecutive nights with the NASA IRTF in January of 2022 (Sec. 2). This data was reduced (Fig. 1) and paired with older data taken in 2008 using Keck and 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 12 CO rovibrational transitions and the hydrogen Pfβ line (Fig.  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 (Sec. 3.2.1) and an eccentric model with two components (Sec. 3.2.2). From these fits, we constrained the disk properties (Table 4), such as the temperature (Eqn. 10) and surface density (Eqn. 11) profiles, and the eccentricities (e in and e out ) of the different components. Our main conclusions are highlighted below: • The core of the average stacked 12 CO line profile is blue-shifted while the wings are red-shifted (Sec. 3.2.2). This is well modeled by introducing an inner and outer disk component that are both eccentric but oppositely oriented (Fig. 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 (Sec. 3.3).
• A 9 d variability was observed in the hydrogen Pfβ line (Fig. 3), which indicates variations in the stellar accretion rate (Sec. 3.1). A 9 d 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).
• The 12 CO lines displayed varying asymmetries over hundreds of orbits (Fig. 2). This might indicate that the inner disk around CI Tau is eccentric and precessing (Sec. 3).
• We did not detect 12 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 (Sec. 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 −8 Jupiter-mass per year (Sec. 3.1).
• 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.
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 to further investigate this system.
The CO and hydrogen equivalent widths (Table 3; Fig.  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 12 CO emissions (Sec. 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 (Sec. 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.

ACKNOWLEDGEMENTS
We thank Jean-Francois Donati for a helpful discussion and providing input on the CI Tau system. This work includes data gathered at the Infrared Tele-scope 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.