High-resolution Millimeter Imaging of the CI Tau Protoplanetary Disk: A Massive Ensemble of Protoplanets from 0.1 to 100 au

Since the 1995 discovery of the first hot Jupiter (Mayor & Queloz 1995), it is now established that such gas giant planets orbiting at radii <0.1 au from their parent stars are found in around 1% of main-sequence solar-type stars (Wright et al. 2012). There is considerable debate as to whether these objects formed in situ or have instead migrated from larger radii, either from interaction with their natal protoplanetary disk (Kley & Nelson 2012) or planet–planet scattering after the disk has dispersed (Rasio & Ford 1996). With typical ages of up to several Gyr, most hot Jupiter hosts have long since lost their protoplanetary disks (typical lifetime of a few Myr; Haisch et al. 2001); arguments about the origin of hot Jupiters are thus usually based on theoretical models linking hypothetical initial conditions to present-day orbital parameters.

The recent discovery (Johns-Krull et al. 2016; Biddle et al. 2018), using the radial velocity technique, of a hot Jupiter in the young disk-bearing solar-type star CI Tau, has demonstrated that in at least this case the hot Jupiter is already in a very close orbit when the star is only ∼2 Myr old (Guilloteau et al. 2014). CI Tau is a well-studied system, with mass 0.92 M_⊙ (Simon et al. 2017), luminosity 0.93 L_⊙ (Guilloteau et al. 2014), and is already known to host a massive dust and gas disk extending many hundreds of au from the star (Andrews & Williams 2007; Guilloteau et al. 2011); additionally it displays a high accretion rate of gas onto the star (McClure et al. 2013). The hot Jupiter's mass is ∼11.3 M_Jupiter if its orbit is aligned with the outer disk (Guilloteau et al. 2014), consistent with the orbital alignment between hot Jupiters and outer planets found in mature exoplanetary systems (Becker et al. 2017). This mass places it in the top 5% of the main-sequence hot Jupiter population. Around half of mature hot Jupiter systems also contain companions (Knutson et al. 2014; Ngo et al. 2015) at less than 20 au which, if present at early times, would create structure in the protoplanetary disk. Although previous submillimeter observations have hinted at a possible gap in the disk around CI Tau at 100 au, they lacked the resolution to characterize this in detail or probe the inner disk where companions may be expected (Konishi et al. 2018).

Here we present high-resolution 1.3 mm Atacama Large Millimeter/submillimeter Array (ALMA) imaging of the disk surrounding CI Tau and report three pronounced annular gaps in emission between 10 and 100 au. We present visibility modeling and explore the origin of these structures via hydrodynamical modeling, arguing for the presence of three gas giants as outer companions to the hot Jupiter.

CI Tau was observed with ALMA on 2017 September 23 and 24 (Project ID: 2016.1.01370.S, PI: C. J. Clarke) with 40 antennas (baselines between 21 and 12145.2 m) and an on-source integration time of 32.35 minutes in both cases. The correlator used four spectral windows centered at 224, 226, 240, and 242 GHz in time division mode to measure the continuum in Band 6. Each spectral window used 128 channels and a bandwidth of 1.875 GHz, together providing an effective total continuum bandwidth of 7.5 GHz. To calibrate the visibilities, a set of standard calibrators were also observed. Calibration of the complex interferometric visibilities used the Common Astronomy Software Applications (CASA) v5.1.1 and the ALMA Pipeline. In panel (a) of Figure 1 we imaged the calibrated visibilities using the multi-scale CLEAN algorithm with scale parameters of 0'', 0048, 008, 024, and Briggs weighting with a robust parameter of 0.5 to obtain the optimal signal-to-noise ratio (S/N) and spatial resolution. The resulting synthesized beam is 005 × 003 with a position angle of 16°, while the achieved rms noise level is 13 μJy beam⁻¹.

Zoom In Zoom Out Reset image size

We characterize the CI Tau brightness by fitting the continuum visibilities with an axisymmetric parametric model consisting of an envelope and three gaps. For the envelope we use an exponentially tapered power law

$$\begin{eqnarray}&&{I}_{\mathrm{env}}(R)={I}_{0}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{{\gamma }_{1}}\exp \left[-{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{{\gamma }_{2}}\right]\quad \mathrm{for}\quad R\leqslant {R}_{\mathrm{out}}\end{eqnarray} \tag{ 1 }$$

with I₀ a brightness normalization constant (Jy sr⁻¹). Each gap is parametrized as the difference of two logistic functions:

$$\begin{eqnarray}&&{I}_{\mathrm{gap}}(R)={\delta }_{\mathrm{gap}}\left[\displaystyle \frac{1}{1+{e}^{-{k}_{2}(R-{R}_{{\rm{g}}}-{w}_{2})}}-\displaystyle \frac{1}{1+{e}^{-{k}_{1}(R-{R}_{{\rm{g}}}+{w}_{1})}}\right],\end{eqnarray} \tag{ 2 }$$

where δ_gap describes the gap depth (δ_gap = 0 corresponding to no gap), R_g is the gap radial location, w₁ and w₂ are the left- and right-hand gap widths at half depth and k₁, k₂ express the steepness of the left and right gap profile. The brightness profile is given by

$$\begin{eqnarray}&&I(R)={I}_{\mathrm{env}}(1+{I}_{\mathrm{gap},1})(1+{I}_{\mathrm{gap},2})(1+{I}_{\mathrm{gap},3}),\end{eqnarray} \tag{ 3 }$$

involving five free parameters for I_env and six free parameters for each I_gap. We simultaneously fit the disk inclination i and position angle (PA; defined east of north) and the offset (ΔR.A., Δdecl.) from the phase center. We thus have a parameter set θ = (I₀, R_c, γ₁, γ₂, ...) described by 23 parameters for the brightness profile plus four for the system geometry. The computation of the visibilities V_mod for each model $\theta$ is performed using GALARIO⁷ (Tazzari et al. 2018), which first computes the 2D image of the disk for a given I(R) and then Fourier transforms and samples it in the observed (u, v) points. The likelihood of the observations V_obs given the model visibilities ${V}_{\mathrm{mod}}(\theta )$ is assumed Gaussian:

$$\begin{eqnarray}&&\mathrm{log}p({V}_{\mathrm{obs}}\,| \,\theta )=-\displaystyle \frac{1}{2}\,{\chi }^{2}=-\displaystyle \frac{1}{2}\displaystyle \sum _{j=1}^{N}| {V}_{\mathrm{mod}}-{V}_{\mathrm{obs}}{| }^{2}{w}_{j},\end{eqnarray} \tag{ 4 }$$

where N is the total number of visibility points and w_j is the weight⁸ of the jth visibility. The parameter space is explored with a Bayesian approach using the emcee Markov Chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013), providing an estimate of the posterior probability distribution of the model parameters given the observations

$$\begin{eqnarray}&&\mathrm{log}p(\theta \,| \,{V}_{\mathrm{obs}})=\mathrm{log}p({V}_{\mathrm{obs}}\,| \,\theta )+\mathrm{log}p(\theta )+C,\end{eqnarray} \tag{ 5 }$$

where p(θ) is the prior on the parameters, and C is a normalization constant. Because the parameters are independent, the priors can be written as $p(\theta )={\prod }_{i}p({\theta }_{i})$. We choose uniform priors on all parameters, except for the inclination for which $p(i)=\sin (i)$ for 0 ≤ i ≤ π/2.

We run the MCMC sampler with 120 walkers for 50 × 10³ steps after a burn in phase of 30 × 10³ steps. We assessed convergence through visual inspection of the chains' trace plots and also by estimating the autocorrelation time (Foreman-Mackey et al. 2013), resulting in ∼150 steps on average for all parameters. From the 6 × 10⁶ samples in the MCMC chain, we select as best-fit model the maximum likelihood model, i.e., that with lowest normalized χ² ≃ 1, as given by the following parameters: I₀ = 10.72 Jy sr⁻¹, R_c = 046, γ₁ = −0.39, γ₂ = 1.50, R_g,1 = 012, w_1g,1 = 005, w_2g,1 = 001, k_1g,1 = 222.04 arcsec⁻¹, k_2g,1 = 52.82 arcsec⁻¹, δ_g,1 = 0.98, R_g,2 = 025, w_1g,2 = 029, w_2g,2 = 011, k_1g,2 = 95.08 arcsec⁻¹, k_2g,2 = 59.77 arcsec⁻¹, δ_g,2 = 0.70, R_g,3 = 088, w_1g,3 = 064, w_2g,3 = 006, k_1g,3 = 11.68 arcsec⁻¹, k_2g,3 = 21.87 arcsec⁻¹, δ_g,3 = 0.84, R_out = 277, i = 4924, PA = 1128, ΔR.A. = 033, Δdecl. = −009. This maximum likelihood model falls in the central 68% interval of the posterior distribution of all parameters and, indeed, its brightness profile is representative of the density of models generated by the posterior (see panel (d) of Figure 1).

In panel (b) of Figure 1 we compare the observed visibilities and those of the best-fit model as a function of deprojected baseline. Panel (d) shows a family of 5 × 10³ models drawn from the inferred posterior (red lines) and the best-fit model (black dotted–dashed line): assuming a distance of 140 pc (Guilloteau et al. 2014), the brightness profile is tightly constrained between 20 and 100 au (i.e., the spatial scales probed by most of the interferometric baselines in the data set) and more uncertain for R < 20 au. Thus we cannot firmly constrain the detailed shape of the innermost gap, whose width is comparable to the beam (∼7 au). We explored in greater detail this degeneracy with a dedicated model suite and found an upper limit on the ratio of the flux inside to outside the gap of ∼0.28: the gray shaded area in panel (d) highlights the range of brightness values that is compatible with the data.

The synthesized image of the residuals obtained for the best-fit model is shown in panel (c) of Figure 1: there are virtually no residuals (<3σ) in most of the disk at radii R > 25 au, confirming the axisymmetry of the brightness profile. The residuals are most significant (up to a 12σ level) at the disk center and in the north–west of the innermost ring. The central residuals reflect the fact that the functional form that we have adopted is insufficiently flexible to correctly capture the emissivity profile in the innermost disk. The latter non-axisymmetric residuals might be caused by a combination of optical depth effects owing to the viewing angle of the observations (i ∼ 49°) and a genuine difference in the local dust temperature.

We note that any perturbation of the disk caused by the hot Jupiter at 0.1 au would occur on a scale of a few times its orbital radius and would thus be indistinguishable within the synthesized beam of 7 × 4 au.

Structure in protoplanetary disks can derive from many causes. The non-planetary mechanisms that have been proposed to date, however, are not well matched to CI Tau (e.g., photoevaporation produces holes rather than gaps; Ercolano et al. 2017), while simulations of the vertical shear instability (Flock et al. 2017) and of non-ideal magnetohydrodynamic (MHD) effects (Flock et al. 2015) do not produce the well-spaced multiple rings seen in CI Tau. While gaps may also arise from opacity effects associated with ice sublimation fronts (Zhang et al. 2016), the outermost two rings are well outside of the sublimation fronts of even the least volatile species, N₂ and CO (Kwon et al. 2015). Moreover, in the innermost gap, our modeling implies a depletion of the optical depth by a factor ∼50 compared with adjacent regions, considerably more than can be attributed to opacity variations. We therefore focus on the planetary hypothesis.

While gap width can be used to infer the required planet mass (Rosotti et al. 2016), this conversion depends on the turbulence in the disk, as parameterized by the Shakura–Sunyaev α parameter (Shakura & Sunyaev 1973), which controls the transport properties of both dust and gas. The level of disk turbulence is difficult to constrain observationally: some estimates based on turbulent line broadening have suggested very low α values that then struggle to reproduce observed accretion rates (Flaherty et al. 2015), although the universality of this result has been questioned (e.g., Teague et al. 2016).

We explore hydrodynamical models in which α and the gas surface density are constrained by the observed high accretion rate onto the star ($\dot{M}=3\times {10}^{-8}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ McClure et al. 2013), assuming that this accretion is driven by some form of turbulent viscosity; moreover, the highly axisymmetric image implies that the disk is not self-gravitating. Together these constraints favor a rather high α value (>10⁻²). We also assume that the maximum grain size, a_max, is locally determined by the minimum of two limits imposed by radial drift and fragmentation, assuming a fragmentation velocity of 10 m s⁻¹ (Birnstiel et al. 2012). We compute the corresponding dust opacity (Tazzari et al. 2016), assuming a population of compact silicate grains with size distribution dn/da ∝ a⁻³ for a < a_max, and in our emissivity modeling adopt the temperature profile derived by Kwon et al. (2015). We assume an initial dust to gas ratio of 0.01.

Below we describe simulations of our fiducial model with parameters α = 0.014 and total disk mass within 200 au of 20 M_Jupiter. These parameters imply strong accretion in the disk (∼10⁻⁸ M_⊙ yr⁻¹) and yield a profile of a_max that is compatible with measurements of the disk-averaged spectral index in CI Tau.⁹

We use these parameters in hydrodynamical simulations where we insert three planets in the disk, using the 2D version of the FARGO3D code (Benéz-Llambay & Masset 2016) with our implementation of drag-coupled dust (Rosotti et al. 2016). We employ 350 logarithmically spaced cells in the radial direction (from 5.6 to 378 au) and 512 cells azimuthally, producing approximately square cells at each location. We adopt an initial gas surface density profile Σ_gas ∝ 1/r normalized at 8 g cm⁻² at 25 au, steepening to a 1/r² profile beyond 60 au. The local value of a_max is computed as above.

We then compute a synthetic emissivity profile (using the temperature profile and calculation of opacity as a function of a_max detailed above) for direct comparison with our GALARIO-derived profile. The purple dashed line in panel (d) of Figure 1 presents the brightness profile of our fiducial model, where planets of mass 0.75, 0.15 and 0.4 M_Jupiter are located at orbital radii of 14, 43, and 108 au. We also produce a synthesized image (Figure 2), generating model visibilities via the ft task in CASA with exactly the same uv-plane coverage and observational setup as the actual observations, and then CLEANing the image using the same imaging parameters as the observed image. Note that for a given disk model, the planet mass within the innermost gap is only determined to within around a factor of two because this gap is poorly resolved, while the values are constrained to within around 30% in the outer two gaps.

Zoom In Zoom Out Reset image size

5.1. Observational Tests of the Fiducial Model

5.2. Evolutionary Scenarios for the Fiducial Model: Formation and Migration

5.3. Comparison with Other Gapped Disks

High-resolution ALMA data of the disk in the young star CI Tau has revealed three prominent annular emission gaps which we have interpreted as an ensemble of massive planets spanning a factor 1000 in orbital radius. The wealth of supplementary data available on CI Tau has allowed us to construct models that are consistent with all of the data on this system available to date. The inferred planetary architecture suggests that the observed association between hot Jupiter and other companions may be in place at very early times. We note that the outer two planets (sub-Jovian planets at radii of 43 and 108 au) present a challenge to current planet formation models.

This work is supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. M.K. acknowledges funding from the European Unions Horizon 2020 program, Marie Sklodowska-Curie grant agreement No. 753799, and F.M. from a Leverhulme Trust Early Career Fellowship, the Isaac Newton Trust and a Royal Society Dorothy Hodgkin Fellowship. This Letter uses ALMA data ADS/JAO.ALMA#2016.1.01370.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Software: GALARIO v1.2 (Tazzari et al. 2018), emcee v2.2.0 (Foreman-Mackey et al. 2013), CASA v5.1.1 (McMullin et al. 2007).