CI Tau: A controlled experiment in disk-planet interaction

CI Tau is a young (~2 Myr) T Tauri system with a substantial near-infrared excess in its SED, indicating that the protoplanetary disk extends very close to its star. This is seemingly at odds with the radial-velocity discovery of CI Tau b, a ~12 MJ planet at ~0.1 au, which would be expected to carve a wide, deep cavity in the innermost disk. To investigate this apparent contradiction, we run 2D hydrodynamics simulations to study the effect of the planet on the disk, then post-process the results with radiative transfer to obtain an SED. We find that at ~0.1 au, even such a massive companion has little impact on the near-infrared excess, a result that holds regardless of planetary eccentricity and dust size distribution. Conversely, the observed full-disk signature in CI Tau's SED is consistent with the existence of the hot super-Jupiter CI Tau b. As our simulations uncover, clear transition-disk signatures in SEDs are more likely to be signposts of nascent"warm"Jupiters, located at around 1 AU in the future habitable zones of their host stars.


INTRODUCTION
The spectral energy distribution (SED) of a planetforming system consists of three basic components: the stellar photosphere, accretion luminosity, and an excess at infrared and longer wavelengths emitted by dust in the protoplanetary disk. In some cases, the nearinfrared (NIR) portion of the disk excess is reduced, implying that dust has been depleted in the innermost, hottest part of the disk. Such SEDs are thought to reflect a "transitional" phase in disk evolution (e.g., Calvet et al. 2002Calvet et al. , 2005 where they clear from the inside out, a picture increasingly borne out in recent years by high-resolution observations (e.g., Andrews et al. 2011;Espaillat et al. 2014;van der Marel et al. 2018) Given their ubiquity, planets are expected to play a key role in carving transition cavities, and conversely, transition cavities are thought to be signposts of planet formation. These prospects have received substantial theoretical (Dodson-Robinson & Salyk 2011;Zhu et al. 2011Zhu et al. , 2012Dong et al. 2012a;Muley et al. 2019;Bae et al. 2019) and observational (Sallum et al. 2015;Keppler et al. 2018) interest. However, planets are difficult to directly image in transition disks, and often have large, model-dependent uncertainties in mass that make their role in carving a given cavity unclear (see Müller et al. 2018;Asensio-Torres et al. 2021, and references therein).
In this context, the ∼2 Myr old CI Tauri system is of particular interest. It is one of the rare protoplanetarydisk systems where a super-Jupiter planet candidate has been detected independently of imaging-in this case, by radial velocity (Johns-Krull et al. 2016). The detection of CI Tau b is not uncontroversial at the moment (Donati et al. 2020). The companion is predicted to have a mass of ∼12M J , a semimajor axis of 0.08 au, and an eccentricity of 0.25 ± 0.15 (Flagg et al. 2019;Donati et al. 2020). However, the NIR excess of the system is characteristic of a full, unperturbed disk, with no obvious signs of disk-planet interactions (e.g., McClure et al. 2013a,b). Determining that the cavity excavated by CI Tau b would not substantially impact the SED would resolve this apparent contradiction. This would provide a starting point to test which disk and planet parameters do impact the SED, and more broadly, under what circumstances transitional disk-like SEDs can be seen as signposts of inner massive planets. The existence of hot Jupiters in ∼Myr old systems would pose a direct test of formation models such as in-situ growth or disk migration (Dawson & Johnson 2018).
To this end, we carry out hydrodynamical simulations of disk-planet interactions, then post-process the results with a radiative-transfer code, in order to obtain disk structures and SEDs. We vary parameters such as planetary semimajor axis, eccentricity, and dust size distribution. By focusing on cavities 5 au in width extending all the way to the star, and by self-consistently importing a dust distribution from the hydrodynamics, our study complements previous work by Zhu et al. (2012) and Dong et al. (2012b), who used semi-analytic modeling to study the SED of disks with ∼20 au and ∼65 au cavities respectively. Section 2 contains our methods, Section 3 our results, and Section 4 our conclusions.

Full-disk case
Before investigating the effect of a planet, we first use the HOCHUNK3D Monte Carlo ray-tracing code (Whitney et al. 2013) to find an unperturbed, analytical disk model that best reproduces the observed photometric SED. We largely follow the CI Tau model in Ballering & Eisner (2019) in this step. We use a photospheric temperature of T * = 4000 K and a total luminosity of L tot = 1.18L (Andrews et al. 2013); 1 we assume that 75% of this luminosity (L * = 0.89L ) comes from direct photospheric emission, while 25% (L acc = 0.30L ) comes from accretion onto the star, in order to match the modeling work of McClure et al. (2013a) at λ < 1µm. The photospheric luminosity corresponds to a stellar radius of R * = 1.96R , typical of T Tauri stars. Using this radius, and using a stellar mass of M * = 0.9M (Guilloteau et al. 2014;Flagg et al. 2019) we estimate a mass accretion rate ofṀ acc = 2.09 × 10 −8 M yr −1 , within order-unity of that in McClure et al. (2013a). As they did, we model the accretion luminosity as an 8000 K blackbody emanating evenly from the stellar surface, and find likewise that we reproduce the observed SED in the optical. We use 5×10 8 photon packets to perform ray-tracing.
Our models use Kim et al. (1994) small interstellar medium grains (hereafter ISM), which provide a good fit to the SED photometry for wavelengths λ 100µm, as well as the 10 µm feature observed by Spitzer/IRS. Such grains, with Stokes numbers St 1, are well-coupled to the gas, and so do not require separate modeling in our hydrodynamics simulations (Section 2.2). We use a dust sublimation temperature T sub = 1700 K, with a sublimation radius at r sub = 0.08 au and a curved wall with r curve = 0.003 au to account for pressure-mediated vertical dependence in the sublimation radius (Isella & Natta 2005). We assume a dust-to-gas ratio of 0.01, and fiducially assign 10% of this mass to small grains (Dong et al. 2012b); we plot the relevant opacity in Figure 1. The rest of the dust mass is assumed to be in bigger grains. Because our focus is on the NIR excess, we need not model this population of larger grains, which supply only a negligibly small fraction of total opacity at λ 10µm.
Our disk is flared, with an aspect ratio of where c s is the sound speed, v k the Keplerian velocity, r the distance from the star, H 0 = 0.065 the aspect ratio at 1 au, and β = 0.125 the flaring index. These parameters are constrained by fitting the disk SED (Ballering & Eisner 2019). The disk profile thus does not reflect that obtained by simply taking the radiative transfer output temperature and assuming hydrostatic equilibrium. Such parametrizations are often adopted in disk modeling (e.g., Andrews et al. 2011;Dong et al. 2012b;Zhang et al. 2014), as certain physical processes that affect the thermal state of the disk, including planetary shocks and viscous heating in a deep cavity, are not readily and self-consistently treated in radiative transfer codes. We use a total (gas and dust) mass of M disk = 0.02M inside its outer boundary r out = 100 au. The surface density scales as where Σ 0 = 283 g/cm 2 and r trunc = 0.08 au-where the rotation period of the star (and therefore, that of its frozen-in magnetic field) coincides with the Keplerian period, leading to a magnetospheric truncation of the disk. This happens to coincide with the dust sublimation radius.
2.2. Planet-carved cavity case Figure 1. Opacity for our dust prescription, at a typical dust-to-gas ratio of 0.01 with 10% small grains by mass. The lack of large grains means that the optical thickness and emission of our disk at wavelengths longer than ∼ 10µm, which we do not focus on, is understated.
To simulate the effect of CI Tau b on its disk, we run 2D hydrodynamics simulations with PEnGUIn (Fung 2015), a GPU-accelerated, Lagrangian-remap code that employs the piecewise parabolic method (PPM) of fluid reconstruction. As in Muley et al. (2021), we solve the viscous, compressible Navier-Stokes equations, with gravitational forces from both the star and the planet. We relax the local temperature in each grid cell to the background temperature on a timescale t c = 10 −4 local dynamical times, making our simulations effectively locally isothermal (Miranda & Rafikov 2020). This facilitates comparison with previous simulations of super-Jupiters (e.g, Dunhill et al. 2013;Ragusa et al. 2018;Muley et al. 2019;Teyssandier & Lai 2020), and simplifies the complex thermal physics (Rafikov 2016) and gas-grain coupling issues (Bae et al. 2021) that would emerge for the strong shocks they would generate.
Per Flagg et al. (2019), we set CI Tau b's mass M p = 11.6M J and the stellar mass M * = 0.9M . We use a semimajor axis a p = 0.08 au, and test eccentricities of e p = {0, 0.25, 0.4}, bracketing the e p,obs = 0.25 ± 0.15 range expected from observations. These parameters are fixed throughout the simulation, and do not change due to the gravitational back-reaction of disk material onto the planet; however, similar values are naturally achievable for super-Jupiter planets when the planet is free to migrate (Papaloizou et al. 2001;Bitsch et al. 2013;Dunhill et al. 2013;Ragusa et al. 2018;Muley et al. 2019). White dots indicate planet location; solid lines its orbit; and dotted lines the region between perihelion and aphelion, which is subject to accretional clearing. The faded region is inside the sublimation radius r sub = 0.08 au; we include it for accuracy in our hydrodynamics, but in itself it makes no contribution to the SED.
We allow the planet to accrete disk material according to the prescription of Muley et al. (2019), at a ratė where k = 10, the typical free-fall timescale t ff = 2r 3 acc /GM p , and the typical accretion radius r acc ≡ r H /3 = (r p /3)(M p /M * ) 1/3 (r, φ are cylindrical coordinates of a grid cell within the disk, while r p is the instantaneous radial location and r H is the Hill radius of the planet). This rate is limited, spatially, by a Gaussian function centered on the planet: (4) where we pickṀ max = 25M J /yr, allowing essentially unlimited accretion within the planet's Hill radius. The accreted mass is removed from the domain entirely and is not added to that of the planet.
We use a Shakura-Sunyaev viscosity prescription ν = αc s h with α = 10 −3 , a value commonly used in hydrodynamical simulations of protoplanetary disks (e.g., Fung et al. 2014;Fung & Chiang 2016;Bae et al. 2019) and emerging naturally from small-scale turbulence in nearly-isothermal disks such as the one we use (Manger et al. 2021). Following Ballering & Eisner (2019) we set the background temperature equal to that implied by Equation 1, scaling as T ∝ r −0.75 (note that the disk is being heated by both stellar irradiation and accretion). The initial surface density is given by Equation 2, augmented with an additional exponential truncation factor exp(−(r/r out ) 2 ) for r out = 2.5 au, to minimize spurious addition and removal of mass and angular momentum at the (circularly symmetric) outer boundary. As shown later, we mainly focus on the inner disk structure within 1 AU and photometry at wavelengths shorter than 10µm.
We run hydrodynamics simulations to 25 years, or 1143 orbits of CI Tau b, long enough to establish a quasi-steady state. Our domain spans a radial range of {0.03, 5} au, and the full 2π in azimuth; we use a zero-gradient outflow condition at the inner boundary, and a fixed condition at the outer one. Numerical tests (not shown here) reveal that shifting the inner boundary inward by a factor of 2 does not impact cavity depth. We cover this range with 1000 (r) × 1228 (φ) cells, logarithmically spaced in radius and evenly spaced in azimuth, yielding an effective resolution of ∼9 cells per scale height at the location of the planet. This resolution is similar to that of the CI Tau simulations in Teyssandier & Lai (2020), which in turn is based on detailed numerical tests by Teyssandier & Ogilvie (2017). . Surface density in the innermost 1 au of our hydrodynamical simulations, assuming an ap = 0.08. Thick lines are azimuthally averaged, while thin lines are radial cuts of the closest and farthest extent of the cavity wall. Dotted lines are inside the dust sublimation radius, and therefore make no contribution to the SED. Transluscent bars indicate the perhielion-aphelion range of the planet for each eccentricity. As expected, an eccentric planet creates a wider, and more eccentric, gap.

Hydrodynamics
Planets open gaps in disks by the action of ordinary Lindblad torques on the surrounding disk material. When the planet-star mass ratio q ≡ M p /M * 0.003, a secondary effect-the eccentric Lindblad torquesbecome prominent and excite eccentricity in the disk (Kley & Dirksen 2006). Gas parcels with large semimajor axes can thus intersect the planet's orbit and get consumed by accretion, leaving a wider gap than ordi-nary Lindblad torques alone would indicate. These effects are enhanced substantially when the planet is itself eccentric and exposed to a wider swathe of the disk, as visible in Figure 2. This is similar to the result in Muley et al. (2019), but the locking of disk and planet pericenters means the disk remains eccentric in our present work, rather than circularizing as in theirs. Our simulations also show "horseshoes", azimuthal asymmetries in surface density formed by companions with mass ratios q 0.01 (Ragusa et al. 2020).
In Figure 3, we show more quantitative disk surface density profiles for all the eccentricities we test, with radial cuts of the closest and farthest approaches of the cavity edge. Our plot assumes the fiducial a p = 0.08 au, corresponding to CI Tau's orbital parameters; we note that maintaining a disk mass of M disk = 0.02M means that plotted densities scale down, and distances scale up, by factors of 10 0.5 (a p = 0.26) and 10 (a p = 0.81). For large eccentricities, the cavity width is up to 5× the planetary semimajor axis at its widest point, and can be depleted by a factor of ∼10 2 − 10 3 ; in the zeroeccentricity case, the cavity is only about 2.5× as wide as the planetary, but deeper than in the eccentric case. Our cavity depths are broadly consistent with those of Teyssandier & Lai (2020), who focused on pulsed accretion, also observed here. However, we concentrate on gap depths, so we use a zero-gradient inner boundary that allows material to flow into and out of the domain. It is not as well-suited to studying pulsed accretion as their "diode" boundary, which only allows material to exit; this enforces anṀ disk ≤ 0 at all times, at the cost of some excess cavity depletion.

Radiative transfer
In Figure 4, we plot model SEDs for the three planetary eccentricities we test, along with a full-disk SED, and compare them with observed photometry points (Andrews et al. 2013). For the fiducial cavity (a p = 0.08 au, leading to r cav ≈ 0.2 − 0.4 au, increasing with eccentricity), shown in the upper panel, eccentricity has limited observational impact, with the full disk and e p = {0, 0.25, 0.4} cases all agreeing with each other, as well as the observed CI Tau SED at λ 10µm. We conclude that the existence of a companion in CI Tau is compatible with the system's SED. Zhu et al. (2012) used hydrodynamics and radiative transfer to study transition cavities between ∼2 − 20 au. We connect to their parameter space by rescaling our simulation domain such that a p = 0.26, 0.81 au. The order-of-magnitude range comfortably accommodates not only changes in planetary position, but also any gap-width variations that might be expected due In the upper panel ap = 0.08 AU, as in the case of CI Tau. In the other two panels we explore cavities opened by planets with larger ap. The planet has little effect for ap = 0.08 au, but clearly reduces NIR excess for ap = 0.81 au. In the latter case, the reduction is more pronounced for ep = 0.25, 0.4 than for ep = 0.
to variations in scale height or viscosity (Dong & Fung 2017). When we do so, the differences between the fulldisk and planetary cases become more pronounced than for a p = 0.08. Given the opacities in Figure 1, and the surface densities in Figure 3 (as rescaled for the planet location), planetary depletion is sufficient to make the deepest parts of the cavity optically thin; furthermore, it lowers the altitude of the radial τ = 1 surface, meaning cavity material is less irradiated (e.g., Dong 2015). Both of these effects contribute to a weaker NIR excess at λ 10µm, while the latter opens up the cavity wall to direct stellar irradiation, creating an excess at 10µm λ 100µm; the precise cutoffs shift outward with increasing cavity radius. Effects are somewhat stronger for e p = 0.25, 0.4, on account of the wider cavity in those cases, than for e p = 0. We stress that all of these results were obtained with an accreting, massive ∼12M J companion, and thus represent a best-case scenario for carving cavities with transitional disk-like SEDs.
Besides planetary parameters, dust properties also play a role in transition-disk cavity visibility. Coagulation of dust into larger-sized grains reduces the effective surface area per unit mass, and thus opacity, while fragmentation does the reverse (Birnstiel et al. 2012). In our radiative-transfer models, we parametrize these effects by scaling the small dust abundance to 10×, 1× and 0.1× its fiducial ratio; this is equivalent to assigning 100%, 10%, and 1% of dust mass to small grains, respectively. Alternatively, these rescalings provide a computationally inexpensive way to bracket the effects that a change in gap depth-and thus, viscosity and scale height Fung et al. (2014) -might have on the SED. As elsewhere in this work, we do not include emission from large grains.
We plot our results in Figure 5 for the full disk and all of a p = {0.08, 0.26, 0.81} au, choosing the e p = 0.4 models as they represent the best-case scenario for transition-cavity formation. Given the same fractional cavity depletion, the altitude of the radial τ = 1 surface is more strongly reduced, and more of the cavity becomes optically thin, when small grains are fewer to begin with. As expected, for a p = 0.08, the observational impact of this is minimal, but for a p = 0.81, a small-dust fraction 0.1× our fiducial value can lead to a transition disk-like dip that almost fully exposes the bare stellar photosphere. The same process increases the fraction of the cavity wall that is directly irradiated, creating a larger proportional increase in the excess at 10 − 100µm.

CONCLUSION
Inspired by the CI Tau system, we use hydroynamical simulations and radiative-transfer post-processing to study the properties of cavities carved by closely orbiting super-Jupiters in protoplanetary disks. All of our Figure 5. Test of different small-dust-to-gas fraction; different colors correspond to 10×, 1×, and 0.1× our fiducial small-dust-to-gas ratio of 0.001. All planet models have an ep = 0.4. In the upper panel ap = 0.08 AU, as in the case of CI Tau. In the other two panels we explore cavities opened by planets with larger ap. In each case, we give a comparison to the expected full-disk model in the same color, but with a thin line. The lower the fraction of small grains, the more pronounced the transition-disk cavity's signature.
simulations use the same, freely accreting 11.6M J companion orbiting a 0.9 M star. We test three different orbital eccentricities, e p = {0, 0.25, 0.4}. While the fiducial model mimics the CI Tau system with a planet at 0.08 au, we experiment with different semimajor axis and dust-to-gas ratio by testing a p = {1, 10 0.5 , 10}×0.08 au and f SD = {0.1, 1, 10} × 10 −3 , with respect to the fiducial values. Our SED modeling builds on previous work (e.g., Zhu et al. 2012;Dong et al. 2012b;Muley et al. 2019) both by studying cavities 4 au in radius, and by self-consistently taking dust surface density from the hydrodynamics rather than parametrizing it with piecewise analytic functions.
By design, our setup represents an extremely optimistic scenario for transition-cavity formation. Our cavity is depleted by a factor of ∼10 −4 − 10 −3 and excavated to a width of ∼2.5 − 5× the planet's orbit, the larger figures contingent on the existence of planetary eccentricity. For a planetary semimajor axis ∼1 au, this makes the deepest part of the cavity optically thin, while lowering the radial τ = 1 surface to expose the cavity wall. This, in turn, creates the characteristic transition disk-like dip from 1 − 10µm, and an additional excess at 10µm, which become more pronounced for high planetary eccentricities and low initial fractions of small dust. While the same mechanisms are at play for a closer semimajor axis of ∼0.1 au, the 100×-smaller surface area of the excavated cavity means that the NIR excess is indistinguishable from that of an unperturbed disk-a result that holds even in tests with a cavity artificially depleted by a factor of 10 −6 (not shown here). We thus infer that the SED photometry of the CI Tau system is consistent with the existence of the close-in, multi-Jupiter-mass companion, CI Tau b.
While transitional and pre-transitional SED features are not formed by hot Jupiters per se, they are good signposts for nascent "warm" super-Jupiters in the eventual habitable zones of their host stars. Future theoretical studies of such systems may sample more finely in time, in order to better understand observed variability in near-infrared excess (e.g., Espaillat et al. 2011;Flaherty et al. 2012). Additionally, a coupled dust-gas code would accurately model large-grain dynamics in such systems, and allow radiative transfer to reproduce the long-wavelength part of the SED. Observationally, such systems may be targeted for radial-velocity surveys if appropriately inclined, or for interferometric observations that reveal the telltale eccentric cavity (e.g., Kraus et al. 2013, Fig. 2 in this work). Beyond individual systems, comparing the populations of disks with transitional SED features to those of warm Jupiters-as van der Marel & Mulders (2021) did for resolved transition disks and massive planets more generally-would help elucidate formation mechanisms for these systems.