© P. Curone et al. 2022

1 Introduction

Planets have been detected around very low-mass (VLM) stars, defined as stars with masses ≲0.3 M_⊙ but still above the hydrogen-burning limit of ~0.08 M_⊙, below which the brown dwarf regime begins (Liebert & Probst 1987). Particularly fascinating are the cases of TRAPPIST-1, a ~0.085 M_⊙ star in the solar neighborhood with a system of seven rocky planets (Gillon et al. 2017), and Proxima Centauri, the closest star to the Sun, which has a mass of ~0.12M_⊙ and hosts a terrestrial planet, a super-Earth candidate, and a newly discovered sub-Earth candidate (Anglada-Escudé et al. 2016; Damasso et al. 2020; Faria et al. 2022). In general, data from the Kepler spacecraft show that the occurrence rate of small planets (with radii of 1.0 2.8 R_⊙) is 3.5 times higher for M dwarfs than FGK stars (Mulders et al. 2015). Even more intriguingly, giant planets have been confirmed orbiting around VLM stars and brown dwarfs. Chauvin et al. (2005) directly imaged a ~5 M_Jup planet around a ~25 M_Jup young brown dwarf, and Morales et al. (2019) discovered a planet with a minimum mass of 0.46 M_Jup orbiting a ~0.12 M_⊙ star. However, giant planet formation in this low-stellar-mas s regime remains a conundrum.

Planets originate in protoplanetary disks and, generally, the most supported theory to explain their formation is the core accretion model (Pollack et al. 1996). Assuming this scenario, the main problems for planet formation around VLM stars are the fast dust radial drift (Pinilla et al. 2013) and the apparent lack of material necessary for disks in the low-stellar-mass regime to generate planets (Testi et al. 2016; Sanchis et al. 2020). Several theoretical studies (Payne & Lodato 2007; Liu et al. 2020; Miguel et al. 2020) performed numerical simulations to assess planet formation around VLM stars and brown dwarfs in a core accretion scenario. They have shown that rocky planet formation in this condition is possible, but the emergence of gas giants is always excluded. A viable explanation for the formation of giant planets around VLM stars seems to be the fragmentation of a disk in its early stages due to gravitational instability (Mercer & Stamatellos 2020).

Observational evidence of protoplanetary disks around VLM stars and brown dwarfs has been collected for decades via multiple observatories (e.g., Comeron et al. 1998; Natta & Testi 2001; Natta et al. 2002; Mohanty et al. 2005). Nowadays, high-resolution and high-sensitivity observations from the Atacama Large Millimeter/submillimeter Array (ALMA) can provide essential information to shed light on planet formation in this low-stellar-mass regime (Ricci et al. 2014; Testi et al. 2016). Nonetheless, only a minimal sample of protoplanetary disks in the low-mass regime has been resolved at high resolution. The reason lies in their smaller size and fainter emission compared to disks around more massive stars (Hendler et al. 2017; Rilinger et al. 2019; Sanchis et al. 2020). Hence, they require more demanding observations.

Recently, Kurtovic et al. (2021) presented the first small survey of disks around VLM stars at high resolution observed by ALMA (~0.1″ at 0.87 mm, Band 7). The authors selected a sample of six low-mass disks located in the Taurus star-forming region, focusing on the brightest objects to maximize the likelihood of observing substructures. Three sources showed structures in their dust emission, which could be explained by the interaction with massive planets. Among these disks is CIDA 1 (2MASS J04141760+2806096), the subject of our work.

This source was first observed by Briceno et al. (1993), and the circumstellar disk was detected by Schaefer et al. (2009) in millimeter wavelengths. Then, with observations from ALMA Cycle 0 at an angular resolution of about 0.4″ (corresponding to ~55au at 137.5 pc, the estimated distance to CIDA 1¹), Ricci et al. (2014) showed a dust disk and detected the CO (J = 3−2) line emission.

An inner cavity with a radius of ~20 au in the dust emission was first detected by Pinilla et al. (2018), using observations in Band 7 from ALMA Cycle 3 with a resolution of 0.21″ × 0.12″ (~29 × 17 au). To explain such structure, the authors exploited the Crida et al. (2006) criterion for gap opening in the gas, assuming a Shakura & Sunyaev (1973) viscosity value of 10⁻⁴ and a stellar mass M_* = 0.1 M_⊙. They found that a planet with a minimum mass of ~0.3 M_Jup orbiting at a distance of 15 au from the central star could carve the observed cavity.

The most recent high-resolution observation of CIDA 1 was presented by Pinilla et al. (2021). They showed ALMA Cycle 6 observations (ALMA project 2018.1.00536.S, PI: A. Natta) for dust continuum in Band 7 (0.9 mm) and Band 4 (2.1 mm) at a resolution of ~0.050″ × 0.034″ (~7 × 5 au). These images clearly reveal a dust ring, whose emission peaks at ~20 au from the star, surrounding a gap with an unresolved inner disk (see Fig. 1). Furthermore, these observations also detected the gas line emission from ¹²CO (J = 3−2) and ¹³CO (J = 3−2). The authors used 1D dust evolution simulations and excluded that a planet more massive than M_* = 0.5 M_Jup could cause the observed dust gap.

In this work we aim to reproduce the observed dust and gas emission of the disk to evaluate whether the presence of the observed substructures can be explained by the interactions of the disk with an embedded planet. While Pinilla et al. (2021) assumed an analytical gap shape, we obtain it from hydrodynamical simulations. We employed full 3D modeling to take the dynamical and radiative aspects of the disk into account, along with the effects introduced by interferometric observations with ALMA. Furthermore, in our work we compared the results of our models not only with observed dust emission but also gas emission. We used 3D hydrodynamical simulations to study how planets with different masses influence the dynamics of gas and dust in a disk around a VLM star. After that, we used 3D Monte Carlo radiative transfer simulations to compute the dust continuum emission and the gas line emission. Finally, we treated the output results as real interferometric observations and compared our models to the collected data. The key parameters that we derive are the mass and position of the supposed planet around CIDA 1. The presence of a gas giant planet around a young VLM star would call into question our current theory of planet formation.

This paper is organized as follows: in Sect. 2 we describe the numerical simulations and procedures used to obtain the final synthetic observations. In Sect. 3, we present the results on the modeling of the dust and gas emission, comparing them with the ALMA observations of CIDA 1. Section 4 is devoted to the discussion of our findings: we analyze the simulated spectral index and optical depth, evaluate whether planet-induced perturbations could be detected in the gas line emission, and constrain the minimum and maximum mass of the putative planet. We also compare our modeling with previous mass estimates of the planet and investigate its possible origin. Finally, in Sect. 5 we summarize our results.

Fig. 1 ALMA observation in Band 7 (0.9 mm) of the system CIDA 1 (Pinilla et al. 2021), with a full width at half maximum (FWHM) beam size of 0.050″ × 0.034″, indicated by the white ellipse in the bottom-left corner.

2 Methods

2.1 Gas and Dust Hydrodynamical Simulation

We performed a set of five hydrodynamical simulations with the smoothed particle hydrodynamics (SPH) code PHANTOM (Price et al. 2018b). This code has been extensively used for the modeling of protoplanetary disks, with the intent of reproducing the observed structures in both the dust and gas emission (e.g., Dipierro et al. 2015; Ragusa et al. 2017; Pinte et al. 2018; Price et al. 2018a; Cuello et al. 2019; Ubeira Gabellini et al. 2019; Toci et al. 2020a,b; Veronesi et al. 2020). SPH formulation, and the PHANTOM implementation in particular, includes exact conservation of mass, momentum, and angular momentum, along with a self-consistent computation of the stellar and planet accretion rates and migrations (Price et al. 2018b). We used the one-fluid (Laibe & Price 2014; Ballabio et al. 2018) multigrain (Hutchison et al. 2018) algorithm, which discretizes the disk with a single set of SPH particles, each containing information of gas and dust with different grain sizes. This method is suitable for treating dust grains not fully decoupled from the gas. Such a condition is fulfilled if the midplane Stokes number St ∝ ρ₀a_d/Σ_g < 1, where ρ₀ is the dust grain intrinsic density, a_d the grain size, and Σ_g the gas surface density.

2.1.1 Disk Model

We selected the initial parameters for our simulations, listed in Table 1, based on the observations and estimates of Pinilla et al. (2018, 2021). The system consists of a central star surrounded by a disk composed of gas and dust, sampled by N_SPH = 10⁶ SPH particles. The disk has a gas mass of M_gas = 1.5 × 10⁻³ M_⊙ and initially extends from R_in = 1 au to R_out = 100 au. Following a common practice (e.g., Toci et al. 2020a), we assumed a power law with an inner and an outer taper for the initial gas surface density profile,

$${\displaystyle {\sum }_g\left(R\right)}={{\displaystyle {\sum }_{\text{c}}\left(\frac{R}{R_{\text{in}}}\right)}}^{-p}\exp \left[-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-p}\right]\left(1-\sqrt{\frac{R_{\text{in}}}{R}}\right),$$(1)

where Σ_c is a normalization constant depending on the total disk mass, R_c = 80 au is the characteristic radius of the outer exponential taper, and we set p = 1.5 (as assumed in Pinilla et al. 2018). Our model adopts a locally isothermal equation of state $P={\rho }_{\text{g}}{c}_{\text{s}}^2$, with the following sound speed radial profile:

$$c_{\text{s}}\left(R\right)=c_{\text{s,in}}{\left(\frac{R}{R_{\text{in}}}\right)}^{-q},$$(2)

where ρ_g is the gas volume density, c_s,in is the sound speed at the inner disk radius and q = 0.3 (derived from the spectral energy distribution; see Appendix A). The gas in the disk is in vertical hydrostatic equilibrium, which leads to the following aspect ratio:

$$\frac{H\left(R\right)}{R}=\frac{c_{\text{s}}(R)}{v_{\text{K}}\left(R\right)}=\frac{H\left(R_{\text{in}}\right)}{R_{\text{in}}}{\left(\frac{R}{R_{\text{in}}}\right)}^{1/2-q},$$(3)

where v_K is the Keplerian velocity and the aspect ratio at the inner radius is H(R_in)/R_in = 0.08. With this choice of parameters, the mean value of the ratio between the azimuthally averaged SPH smoothing length 〈h〉 and the scale height H is 〈h〉/H ≈ 0.14, implying that the disk is vertically resolved. We modeled the disk viscosity, which regulates the angular momentum transport, using the SPH artificial viscosity α_AV introduced by Lodato & Price (2010). We set the Shakura & Sunyaev (1973) viscosity α_SS = 5 × 10⁻³, corresponding to α_AV = 0.36. This value directly affects the timescale regulating the disk viscous evolution (higher viscosity generally leading to faster evolution) and on the interaction between the disk and a planet, especially regarding the gas gap carved by a planet (see Sect. 4.4).

The dust component is initially distributed as the gas, following the same surface density profile (Eq. (1)). The initial dust-to-gas mass ratio in the disk is 10⁻² (the standard value assumed for the interstellar medium, Bohlin et al. 1978), and it is constant throughout the disk. We simulated 11 different grain sizes, logarithmically spaced between 0.16 urn and 1 mm. The grain size distribution follows the power law $\text{d}N\left(a_{\text{d}}\right)/\text{d}{a}_{\text{d}}\propto a_{\text{d}}^{-\xi }$, truncated at the minimum and maximum grain size, where a_d indicates the grain size and N(a_d) is the number of dust grains per grain size. In our simulation, we assumed the typical value ξ = 3.5 for the distribution power law (Mathis et al. 1977). Simulating various grain sizes allows different levels of coupling between dust and gas to be tested. The intrinsic density of dust grains is ρ₀ = 2 g cm⁻³.

Fig. 2 Face-on gas and dust surface density, in units of g cm−2, from the Phantom simulation with an initial planet mass of 1.5 MJup and a final planet mass of 2.0 MJup after ≈4 × 104 yr (≈560 planet orbits). Gas density is displayed in the top-left corner, and the remaining panels present the surface density for the 11 simulated dust grain sizes. The bigger central green dot represents the star, and the smaller green dot indicates the planet position. Images were produced using the tool Splash (Price 2007).

2.1.2 Properties of the Central Star and the Embedded Planet

The central star and the planet are modeled as sink particles, which are free to migrate and accrete gas and dust (see Appendix B). The stellar mass of CIDA 1 is known to be in the range ~0.1−0.2 M_⊙ (Pinilla et al. 2021). Kurtovic et al. (2021) estimate a mass of −0.19 M_⊙ using the Pascucci et al. (2016) method with the distance of 137.5 pc (Pinilla et al. 2021). In our simulations, we adopted a stellar mass M_* = 0.2 M_⊙. The star accretion radius is R_*acc = 1 au and defines a region where incoming SPH particles are considered accreted onto the star. The planet accretion radius is chosen to be R_Pacc = 1/4 R_H, where R_H is the planet Hill radius, defined as

$$R_{\text{H}}={\left(\frac{1}{3}\frac{M_{\text{P}}}{M_{\star }}\right)}^{1/3}{d}_{\text{p}},$$(4)

where d_P is the distance of the planet from the central star.

We performed a set of SPH simulations, each with an embedded planet starting at a distance d_P,0 = 10.0 au from the star. We varied the initial planet mass M_P,0: the chosen values are 0.1, 0.5, 1.0, 1.5, and 2.0 M_Jup. We let the simulations evolve until they reach a stable configuration. We analyzed the simulation results at 4 × 10⁴ yr, corresponding to about 560 planet orbits. To check that no major transients are neglected with this choice, we also computed the dust emission synthetic observations for a longer evolution time of 8 × 10⁴ yr, equivalent to about 1120 planet orbits.

For illustrative purposes, we show in Fig. 2 the face-on surface densities of the gas and the 11 dust populations from the simulation with initial planet mass M_P,0 = 1.5 M_Jup, evaluated after 4 × 10⁴ yr. Here, the dust radial drift and the dust trap induced by the planet (Pinilla et al. 2012b) can be appreciated. This is because the massive planet creates a pressure bump in the gas outside of its orbit, which attracts nearby dust and intercepts dust migrating from the outer regions toward the central star. The faster motion of larger grains (from ~0.1 to 1 mm) results in ring-shaped dust distributions, whereas smaller grains present wider radial profiles of their density, because they have a slower radial drift being more coupled to the gas.

After this first suite of simulations, we aim at better constraining the minimum and maximum planet mass able to reproduce the observed dust emission by running additional simulations. In these cases, the embedded planet is placed at different initial radial distances from the star to find the best match with the peak intensity of the observed ring. The minimum mass of the planet is discussed in Sect. 4.4, where we present the results of two simulations: the first with M_P,0 = 0.5 M_Jup and d_P,0 = 12.0 au, the second with M_P,0 = 1.0 M_Jup and d_P,0 = 11.0 au. The maximum mass of the planet is assessed in Sect. 4.5, considering a simulation with M_P,0 = 4.0 M_Jup and d_P,0 = 9.0 au.

2.1.3 Limitations of the Modeling

High-resolution ALMA images of the dust continuum emission from CIDA 1 (Fig. 1) show a dust ring peaking at −20 au from the star, a large gap, and an unresolved inner disk at the center, indicating the presence of dust within a few astronomical units of the star. We focus on modeling the external dust ring, whereas reproducing the inner disk is beyond the aim of this work, given the setup we use for the SPH simulations.

In our models, we adopted R_*acc = 1 au; that is, we are not simulating the dynamics of the material within a radius of 1 au from the central star. Adopting a smaller Racc rapidly causes a substantial increase in computational time due to a higher dynamic range. SPH particles closer to the central star move faster, so smaller time steps and a greater computational cost are needed to simulate their motion.

Furthermore, an intrinsic consequence of the artificial viscosity implementation in Phantom is that α_ss ∝ 〈h〉 ∝ ρ^−1/3 (in 3D), meaning that the disk viscosity increases wherever the density decreases (Lodato & Price 2010). A local increase in viscosity leads to a shorter viscous timescale and, therefore, a faster transport of angular momentum in the disk. This description applies to our simulations: if the embedded planet is massive enough to carve a gap in the disk, in that region, the viscosity rises, and this eventually brings to a faster depletion of the inner material accreting onto the central star.

However, these numerical limitations only apply to the inner area, within the first few astronomical units. They do not affect the simulation in outer regions, where the observed dust ring is located and most of the dust flux and mass lie.

Fig. 3 Refinement of the hydrodynamical model for the radiative transfer simulation via the extrapolation and interpolation of the dust densities for a single SPH particle. In the left panel, starting from dust density values obtained with SPH simulations (red dots), we perform a linear fit for log ρd values associated with the two smallest simulated grain sizes to obtain densities for grains down to 5 nm (blue line). Then, we use a cubic spline interpolation on the remaining log ρd values (orange line). In the right panel, with the linear fit, we obtain ten dust densities related to grain sizes between 5 nm and 0.37 µm (to the left of the vertical dashed gray line), while with the cubic spline interpolation we obtain 50 densities for grain sizes between 0.37 µm and 1 mm (to the right of the vertical dashed gray line).

2.2 Radiative Transfer

2.2.1 Linking Phantom to Polaris

For the radiative transfer simulations, we employ the 3D Monte Carlo code Polaris (Reissl et al. 2016), coupled for the first time with Phantom SPH simulations using the sf3dmodels code as interface (Izquierdo et al. 2018). Polaris is a flexible code that can handle multiple dust compositions, each of them defined by dust material, intrinsic grain density, and grain size distribution between a minimum and maximum grain size. Each dust population is associated with a spatial density distribution. From each SPH simulation, we acquired density distributions for 11 grain sizes, logarithmically spaced between 0.16 µm and 1 mm. To have a more precise model, we obtained more dust density distributions from this initial data. Considering a single SPH particle, at first, we extrapolated densities for smaller grains. Under the assumption that small dust grains well-coupled to the gas maintain their initial size distribution (i.e., $\text{d}N\left(a_{\text{d}}\right)/\text{d}{a}_{\text{d}}\propto a_{\text{d}}^{-3.5}$, which can be rearranged as ${\rho }_{\text{d}}\propto a_{\text{d}}^{0.5}$, with ρ_d indicating the dust mass density), we fitted the values of log₁₀ ρ_d associated with the two smallest simulated grain sizes (0.16 µm and 0.37 µm) using a straight line with slope ~0.5. With this relation, we sampled 10 density values logarithmically spaced between 5 nm and 0.37 µm. For larger dust grains, we performed a cubic spline interpolation for log₁₀ ρ_d values related to the remaining simulated grain sizes. Thanks to this interpolation, we sampled 50 density values, logarithmically spaced between 0.37 um and 1mm. This procedure is illustrated in Fig. 3. We repeated this operation for every SPH particle and obtained in total 60 dust density spatial distributions, 10 for grain sizes between 5 nm and 0.37 um and 50 related to grain sizes between 0.37 um and 1 mm.

While in SPH the information is contained in point particles, Polaris needs data to be sampled on a 3D grid for the radiative transfer simulations. To maintain the spatial distribution of physical properties as close as possible to the simulation output we used a Voronoi grid. Given a distribution of N points, a Voronoi grid is a partition of space into N convex polygonal cells, each containing exactly one generating point such that any portion of space inside a specific cell is closer to its generating point than to any other (Okabe et al. 2000). To avoid numerical artifacts where fewer SPH particles are present, as in the case of gaps and cavities in disks, we introduced a set of “dummy” grid points (as in Izquierdo et al. 2021a), which are artificially added particles to which we assigned a negligible density value (10~⁴² g cm”³). We randomly distributed 1/5 N_Spu dummy points; then, we accepted only those located in low-density areas, less sampled by SPH particles. After that, we used the SPH particles and the accepted dummy points as generating points for the Voronoi grid.

2.2.2 Polaris Simulations

With the 60 dust distributions sampled on a Voronoi grid obtained from the initial SPH models, we then performed the radiative transfer simulations with polaris. We considered the same dust composition adopted in Ricci et al. (2010). Here, dust grains are spheres composed of astronomical silicates (10% in volume, optical properties from Draine 2003), carbonaceous materials (20%, Zubko et al. 1996), water ice (30%, Warren & Brandt 2008), and a porosity of 40%.

First, with Polaris, we computed the dust temperature using a Monte Carlo approach. From the estimates of Pinilla et al. (2021), CIDA 1 has a luminosity L_* = 0.15 ± 0.03 L_⊙. In our simulations, we modeled the radiation as 10⁸ photon packages emitted by the central star, which is assumed as a black body with a surface temperature T_eff = 3050 K and a radius R_* = 1.3 R_⊙, whose total luminosity equals the observed value. Second, the dust continuum emission and the gas line emission are computed via ray-tracing. We fixed the distance of CIDA 1 at 137.5 pc. From the analysis of Pinilla et al. (2021), we used an inclination angle i = 37.5° and a position angle PA = 11° for the isk.

The synthetic images are computed at the same wavelengths observed by ALMA. We simulated the dust continuum emission in Band 7 (0.9 mm) and Band 4 (2.1 mm). The gas line emission is calculated for ¹²CO and ¹³CO in the J = 3−2 transition. For the radiative transfer simulation of the gas emission, we assumed that the system is in local thermodynamic equilibrium and that the gas temperature is the same as that of the dust. Works by Bruderer et al. (2012), Bruderer (2013), and Facchini et al. (2018) show that gas and dust in protoplanetary disks in certain conditions may not be fully thermally coupled. However, given the uncertainty in this context, mainly due to a strong dependence on disk turbulence and dust grain size distribution, we made the assumption that gas and dust temperature are coupled and do not introduce further parameters in our simulations. The systemic velocity in our models is fixed at 6.25 km s⁻¹³, because this value led to the best match with the spatial distribution of the observed channel maps. We adopted a gas turbulent velocity of 100 m s⁻¹ . We included freeze-out and photo-dissociation for CO following the parameterizations in Williams & Best (2014). The condition T_gas < 20 K defines the freeze-out region (Bergin et al. 2007). CO is assumed photo-dissociated by UV radiation from the central star and the interstellar radiation field where the gas column density (calculated in the vertical direction) is lower than a critical value N_dissoc = 1.3 × 10²¹ cm⁻² (Visser et al. 2009). Therefore, in the warm molecular layer between the frozen-out midplane and the photo-dissociated regions at higher altitudes, we set the ¹²CO abundance relative to H₂ [¹²CO/H₂] = 5 × 10⁻⁵, and the isotopologue abundance ratio [¹²CO/¹³CO] = 70 (Wilson & Rood 1994); elsewhere, the CO abundance is zero.

To match the velocity resolution of the observed gas channel maps, namely 0.5 km s⁻¹ for ¹²CO and 1.0 km s⁻¹ for ¹³CO, we first computed the simulated channel maps at a higher velocity resolution of 0.045 km s⁻¹. Then, to obtain the final spectrally convoluted synthetic channels, we averaged 11 contiguous simulated ¹²CO channels, centered at the same velocities of the observed channels. We repeated the same procedure for the computed ¹³CO channels, averaging 22 channels to match the velocity resolution of the observations.

2.3 Synthetic Observations

To obtain the final beam-convoluted synthetic images from our simulations, we treated the full-resolution dust continuum images and gas channel maps computed by POLARIS as the actual ALMA observations. Initially, we employed the task samplelmage from the package GALARIO (Tazzari et al. 2018) to compute the synthetic visibilities of the dust continuum images from our models at the same (u, v) points of the ALMA observations, in Band 7 and Band 4, respectively. After that, we calculated the residuals between the observations and the simulations. Finally, we performed the imaging and acquired the synthetic observations using the tclean algorithm from the software CASA, version 5.7 (McMullin et al. 2007). We applied the same parameters used to obtain the dust continuum Band 7 and Band 4 images in Pinilla et al. (2021): the Briggs weighting scheme and a robust parameter of 0.5.

At this point, we can compare the flux in our synthetic images to the one from the observations. We consider the total dust mass of our models as a free parameter that can be rescaled by a constant factor in the radiative transfer simulations to match the flux emitted in Band 7 (0.9 mm) from the observed external dust ring (see Appendix C). It is possible to apply this procedure as long as the dust-to-gas mass ratio remains ≪1, so that the dust back-reaction onto the gas is negligible. In our models, this condition is always fulfilled.

We performed the same method to obtain the ¹²CO and ¹³CO synthetic channel maps, working channel by channel. In the cleaning process, we used a Briggs robust parameter of 1 and applied a uv-tapering with a Gaussian of 0.035″ to recover the same beam size of the observations.

3 Results

The resulting parameters of our models are listed in Table 2, for each simulation with different initial planet mass and for the two considered times t_Sim = 4 × 10⁴ yr and t_Sim = 8 × 10⁴ yr. We note little difference in planet mass (MP) and distance from the central star (dP) in each simulation at the two times, meaning that 4 × 10⁴ yr are sufficient for the planet to reach a slowly evolving configuration. Such a finding is confirmed by the plots in Appendix B showing the time evolution for the planet-star distance and the planet mass in the case of the simulation with an initial planet mass of 1.5 M_Jup. Here, we note a fast transient in the first 10⁴ yr, when the planet is carving the gap, followed by a more stable evolution. At this later stage, the system has adapted from the arbitrary initial conditions and reached a quasi-equilibrium state. If the planet is massive enough, such a situation corresponds to the case in which the planet has already carved an internal cavity both in the gas and in the dust, as in the case illustrated in Fig. 2.

3.1 Dust Continuum Emission

The first aim of our work is reproducing the dust thermal emission of the dust ring in CIDA 1 observed by ALMA in Band 7 (0.9 mm) and Band 4 (2.1 mm), and reported in Pinilla et al. (2021). We present in Fig. 4 the comparison at both wavelengths between the dust continuum image from the ALMA observation and the synthetic image from our simulation with a final planet mass of 2.0 M_Jup after 4 × 10⁴ yr, together with their residuals. In the case of this particular simulation, we note that our model is able to generally reproduce the morphology and the intensity of the observed emission from the ring. As expected, our synthetic images do not present the inner disk, which persists in the residuals. Moreover, we note that the residuals in both bands show some areas, located in the region of the ring and the gap, above the 5σ value or below the −5σ value. Slight shifts in the simulated image centroids could not eliminate these differences. Such inequalities between our model and the observations might be caused by an imperfect estimate of the disk inclination and position angle that we assumed in the simulations or by some actual non-axisymmetric features in the observations.

A thorough comparison between the dust continuum emission from the observations and all our set of simulations after 4 × 10⁴ yr is presented in Fig. 5. We show the azimuthally averaged radial profiles of the intensity in the two bands from the observations, all our simulations, and their residuals.

The smallest simulated planet, with a final mass of 0.2 M_Jup, does not produce any major effect on the disk morphology. The 0.9 M_Jup planet is massive enough to carve a central cavity, but the intensity peak is significantly shifted inward compared to the peak in the observation profiles. The third model, with a planet mass of 1.4 M_Jup, shows a good match with the observation, except for a slight offset in the location of the ring. The two remaining models, hosting the most massive planets with a mass of 2.0 and 2.5 M_Jup, respectively, best reproduce the observed emission in both Band 7 and Band 4. Looking at Table 2 and considering the results after an evolution time of 4 × 10⁴ yr, we note that these last three cases also share consistent estimates of the total dust mass needed to match the observed flux, ~7 × 10⁻⁶ M_⊙, with a dust-to-gas mass ratio of ~9 × 10⁻³.

Results for the dust emission from the simulations after a longer evolution time of 8 × 10⁴ yr are presented in Appendix D. The synthetic profiles for simulations with a minimum planet mass of 2.0 M_Jup still show a good agreement with the observations in Band 7 (which we take as reference; Appendix C). However, the simulated profiles in Band 4 have a slightly lower peak intensity compared to the observed ones because the bigger millimeter-sized dust grains, which emit the most at this wavelength and have the highest radial drift, accreted more onto the central star. Nonetheless, the evident similarities between dust emission after 4 × 10⁴ and 8 × 10⁴ yr prove that, at these times, the disks in our simulations have already overcome the initial fast transient phase and are in a dynamically slowly evolving condition, in the absence of dust evolution.

Fig. 4 Band 7 (0.9 mm, top panels) and Band 4 (2.1 mm, bottom panels) comparison between dust continuum images from the ALMA observation and our synthetic observations obtained from the simulation with a final planet mass of 2.0 MJup after 4 × 104 yr, along with their residuals. The synthesized beam (0.050″ × 0.034″ for Band 7, 0.048″ × 0.032″ for Band 4, FWHMs) is shown in the bottom-left corner of each image. The contours in the residuals represent the 5σ and −5σ flux levels, respectively, for Band 7 and Band 4.

Fig. 5 Comparison between radial profiles of azimuthalfy averaged intensities from the observations, simulations, and corresponding residuals, after the simulations have evolved for 4 × 104 yr. The left-hand panels show results in Band 7, while the right-hand panels refer to Band 4. Each row contains the intensity profiles from the same simulation with a specific final planet mass, ranging from 0.2 MJup (top row) to 2.5 MJup (bottom row). The distance of the planet from the central star is indicated by the vertical dash-dotted brown line. The mean FWHMs of the synthesized beam in Band 7 and Band 4 are shown as horizontal bars in the top row panels. The radial intensity profiles are obtained by deprojecting the disk, using the disk inclination angle i = 37.5°, and dividing it into overlapping annuli, whose widths correspond to the averaged synthesized beam diameter. For each bin, uncertainties are calculated by taking the intensity standard deviation for each pixel in the annulus and dividing it by the square root of the number of synthesized beams contained in the area of the annulus.

Fig. 6 Comparison between the ALMA observation (top panels) and the synthetic images (bottom panels) of CIDA 1 12CO (J = 3−2) channel maps. Color scales are the same in the top and bottom panels. Synthetic images are obtained from the simulation with a final planet mass of 2.0 MJup after an evolution time of 4 × 104 yr. We fixed the systemic velocity at 6.25 km s−1. The contour level traces the 15σ emission from the Band 7 continuum image. The velocity resolution is 0.5 km s−1, and the central velocities of each channel are indicated on the top of each panel. The synthesized beam (0.101″ × 0.075″, FWHM) is shown in the lower-left corner of each panel.

3.2 Gas Line Emission

The dust continuum modeling shown in Fig. 5 proves that, while obtaining a good match in the case of a 1.4 M_Jup planet, the minimum planet mass able to best replicate the observed dust emission is 2.0 M_Jup. Therefore, we decided to use the model with such a planet to compute the ¹²CO (J = 3−2) and ¹³CO (J = 3−2) channel maps.

In Fig. 6, we compare the ¹²CO channel maps observed by ALMA to our synthetic gas emission images. In our model, we set a systemic velocity of 6.25 km s⁻¹. Strong cloud absorption is evident in the central channels of the observation, from ~5.5 to 7.5 km s⁻¹. We analyze the cloud contamination effects in Appendix E, using previous observations of the gas emission from the Taurus molecular cloud by Davis et al. (2010).

Our simulation reproduces appropriately the spatial extent of the channels centered at the velocities v = 4.0, 4.5, 5.0, 8.0, 8.5, 9.0 km s⁻¹. This is more evident in Fig. 7, which shows for these channels the comparison between the observed and simulated ¹²CO emission in terms of intensity contour levels at 5, 10, 15mJy beam⁻¹ and intensity profiles along the PA = 11°. The channels v = 4.5, 5.0, 8.0 km s⁻¹ are characterized by a significant difference in intensity values between the observation and the model. This effect can be explained by a still-present influence of the cloud on these channels, which are just outside the most absorbed velocities. Intensities are more similar in the channels v = 4.0, 8.5, 9.0 km s⁻¹, where the cloud contamination in the observation is weaker, and the emission comes from regions of the disk well sampled by our simulation.

The observed channel maps also show gas emission from the innermost part of the disk, near the inner disk detected in the dust continuum. This material is also the cause of the central emission in high-velocity channels, such as those at v = 3.5 km s⁻¹ and v = 9.5−10.0 km s⁻¹. As we explain in Sect. 2.1.3, our model does not aim at describing this inner region.

We follow the same procedure for the analysis of the ¹³CO emission. Figure 8 presents the comparison between the ¹³CO channel maps observed by ALMA and the synthetic ones computed from our simulation. Here, we have a lower velocity resolution of 1.0 km s⁻¹. In the observations, we note that the 6.5 km s⁻¹ channel is completely affected by cloud absorption, which also impacts the channel at v = 5.5 km s⁻¹. In Fig. 9, we compare the contour levels at 3, 6, 9 mly beam⁻¹ and the intensity profiles at PA = 11° from the observed and simulated ¹³CO emission. Even in this case, our models are able to reproduce the spatial distribution of the observed emission while presenting higher intensity values, except for the 8.5 km s⁻¹ channel. In particular, we note a strong cloud absorption for the observed intensity profile of the 5.5 km s⁻¹ channel, but the corresponding contour plot shows a good agreement in the extent of the emission between observation and simulation. As for the ¹²CO, there is some observed inner emission for the high-velocity channels at 3.5 and 9.5 km s⁻¹, which cannot be reproduced by our model due to its empty inner cavity.

From all of these results, we realize that our single model with a final planet mass of 2.0 M_Jup is able to globally reproduce both the dust and gas emission observed in CIDA 1. In particular, this case perfectly recovers the external dust ring as observed both in Band 7 and Band 4, along with replicating the spatial distribution of the gas emission as traced by ¹²CO (J = 3−2) and ¹³CO (J = 3−2).

Fig. 7 Contour levels (left panel in each subplot) at 5, 10, 15mJy beam−1, and intensity profile comparisons (right panel in each subplot) between the 12CO observation and simulation with a 2.0 MJup planet mass after 4 × 104 yr, in the six channels not completely obscured by the cloud (v = 4.0, 4.5, 5.0, 8.0, 8.5, 9.0 km s−1). Intensity profiles are taken along the PA = 11° of the disk, and the direction is indicated by the dotted gray line in the contour levels plots. The uncertainty in the observation intensity profile is the 1σ noise level (1.6 mJy beam−1 ).

Fig. 8 Comparison between the ALMA observation (top panels) and the synthetic images (bottom panels) of CIDA 1 13CO (J = 3-2) channel maps. Color scales are the same in the top and bottom panels. Synthetic images are obtained from the simulation with a final planet mass of 2.0 MJup after an evolution time of 4 × 104 yr. We fixed the systemic velocity at 6.25 km s−1. The contour level traces the 15σ emission from the Band 7 continuum image. The velocity resolution is 1 km s−1, and the central velocities of each channel are indicated on the top of each panel. The synthesized beam (0.104″ × 0.076″, FWHM) is shown in the lower-left corner of each panel.

4 Discussion

4.1 Stellar Mass and Systemic Velocity

Our models allow fundamental information regarding the CIDA 1 system to be recovered. In the observations, the ¹²CO and ¹³CO channel maps are strongly affected by cloud absorption in the channels near the systemic velocity (Sect. 3.2 and Appendix E). Hence, we compare our synthetic images to the observations only in those channels probing higher velocities in the disk (Figs. 7 and 9). Models are consistent with observations, being able to reproduce the spatial distribution of the emission. The stellar mass = 0.2 M_⊙ and systemic velocity v_sys = 6.25 km s⁻¹ assumed in our simulations are therefore reasonable estimates of these properties in CIDA 1. This stellar mass is compatible with the known range ~0.1−0.2 M_⊙ (Pinilla et al. 2021) and is in very good agreement with the estimate of 0.19 M_⊙ by Kurtovic et al. (2021). We also simulated the case with = 0.1 M_⊙, but the morphology of the resulting gas channel maps differed significantly from the observed ones; therefore, we decided not to include these models in the simulation sample.

4.2 Spectral Index and Optical Depth

Having dust emission data at two different wavelengths allows us to compare the observed and simulated spectral index. Given the fluxes and frequencies corresponding to Band 7 and Band 4, respectively F_B7, v_B7, and F_B4, v_B4, the spectral index α_B7,_B4 is calculated as

$${\alpha }_{\text{B7,B4}}=\frac{{\log }_{10}\left(F_{\text{B7}}/F_{\text{B4}}\right)}{{\log }_{10}\left(v_{\text{B7}}/v_{\text{B4}}\right)}.$$(5)

Assuming that opacity follows a power-law k_ν ∝ ν^β, the emission is optically thin and within the Rayleigh-Jeans regime, then the spectral index can be approximated as α_B7,B4 ≈ 2 + β (Testi et al. 2014). The value of β depends on the maximum dust grain size, while not being influenced by the minimum grain size. A low β (i.e., β < 1, meaning that α_B7,B4 < 3) could be indicative of the presence of grains with size ≳0.1 mm, suggesting that grain growth in the disk may have occurred (Draine 2006). Protoplanetary disks around VLM stars or brown dwarfs with low spectral indices have been observed (e.g., Ricci et al. 2014; Pinilla et al. 2017).

The left and central panels of Fig. 10 present the azimuthally averaged radial profiles of the intensity from our simulations with final planet mass of 2.0 M_]up in Band 7 and Band 4, along with the corresponding profiles of the optical depth along the line of sight. All profiles were produced directly from the full-resolution images computed by Polaris, without any convolution with a synthesized beam. Here, we note that the emission at both wavelengths is partially optically thick, that is, the optical depth τ is above 1. We calculate that ≈51% of the emission in Band 7 is optically thick, whereas this fraction decreases to ≈38% in Band 4. These results are consistent with the findings in Tazzari et al. (2021), showing similar optically thick fractions in a sample of disks located in the Lupus star-forming region. Such a mixture of optically thin and thick emission in our model leads us to interpret with caution the spectral index in CIDA 1. Moreover, the dust temperature predicted in our simulation (~20 K in the disk midplane, where most of the settled dust is located) lies at the limit of the Rayleigh-Jeans domain. In this spectral region, and for a temperature of 20 K, dropping the Rayleigh-Jeans approximation but retaining the optically thin one, we calculate that the spectral index is α_{B7, B4} ≈ 1.7 + β.

Pinilla et al. (2021) find that the spatially integrated spectral index of CIDA 1 is 2.0 ± 0.2 (the uncertainty includes the observed root mean square of the noise and a 10% flux error from calibration). Table 3 compares the flux in Band 7, Band 4, and the spatially integrated spectral index between the observation (considering only flux coming from the external dust ring) and our simulation with a final planet mass of 2.0 M_Jup after of our model is consistent with the observed value within the uncertainties.

A more detailed comparison is presented in the right panel of Fig. 10, which shows the azimuthally averaged radial profiles of the spectral index from the observation and, again, from the simulation with a final planet mass of 2.0 M_Jup after 4 × 10⁴ yr. Our model best reproduces the observed spectral index profile between ~10 au and ~30 au, where the dust ring is detected. Since we focus on reproducing the external dust ring, we do not display the spectral index profiles within a distance of 5 au from the central star. Then, we note that the spectral indices start diverging beyond ~30 au. This effect is due to a slight intensity excess in the model in Band 7 that does not perfectly match the observed Band 7 profile at these radii (see Fig. 5). Between ~14au and ~20au the dust emission is optically thick in both Band 7 and Band 4, meaning that the low spectral index in this region cannot be linked to the presence of millimeter-sized grains. However, this is not the case for the region between ~20 au and ~30 au, still showing a spectral index around ~22.5 associated with optically thin emission. Such behavior is expected in our simulation, where we know that millimeter-sized grains are present at the dust ring location, but this strong match between our model and the observation can indicate that grain growth and migration have occurred in CIDA 1.

Fig. 9 Contour levels (left panel in each subplot) at 3, 6, 9 mJy beam−1, and intensity profile comparisons (right panel in each subplot) between the 13CO observation and simulation with a 2.0MJup planet mass after 4 × 104 yr, in the four channels not completely obscured by the cloud (v = 4.5, 5.5, 7.5, 8.5 km s−1). Intensity profiles are taken along the PA = 11° of the disk, and the direction is indicated by the dotted gray line in the contour levels plots. The uncertainty in the observation intensity profile is the 1σ noise level (1.2 mly beam−1).

Fig. 10 Optical depth in our models and comparison between observed and simulated spectral indices. Left and central panels: Azimuthally averaged radial profiles of the optical depth along the line of sight and the intensities from the full-resolution dust emission image of the simulation (i.e., not convoluted with the synthesized beam) with a final planet mass of 2.0 MJup after 4 × 104 yr, in Band 7 and Band 4. The vertical dash-dotted brown line indicates the location of the planet, and the horizontal dotted gray line marks the optical depth τ = 1 level. The region where the emission is optically thick (optical depth τ ≥ 1) is shown in beige. Right panel: azimuthally averaged radial profiles of the spectral index from the ALMA observation and our simulation with a final planet mass of 2.0 MJup after 4 × 104 yr. Values are obtained using Eq. (5) and the synthetic image intensity profiles in Fig. 5. Uncertainties include the original error in the intensity profiles and a 10% flux error due to calibration. As a reference, the distance of the planet from the central star is indicated by the vertical dash-dotted brown line, and the position of the peak intensity is displayed by the vertical dotted blue line. The beige-shaded area marks the radii where the full-resolution simulated emission is optically thick (τ ≥ 1) in both Band 7 and Band 4.

4.3 Planet-induced Perturbations in Synthetic Gas Channel Maps

The presence of substructures, such as rings, gaps, and cavities, in the dust continuum images of protoplanetary disks can be interpreted as the result of tidal interactions between the disk and a planet. A more direct way to reveal the effects of planet-disk interactions is to look for localized velocity perturbations in the gas channel maps, namely “kinks” (Pinte et al. 2018, 2019, 2020; Izquierdo et al. 2021b; Bollati et al. 2021).

In the channel maps of CIDA 1 observed by ALMA in ¹²CO (J = 3−2) and ¹³CO (J = 3−2) (top panels in Fig. 7 and Fig. 9, respectively), there is no clear sign of a deviation from the Keplerian rotation. Here, we investigate whether, with sufficient sensitivity and both spectral and spatial resolution, a perturbation could be detected in one of the channels that are not affected by cloud absorption.

We consider our simulation with a final planet mass of 2.0 M_Jup after 4 × 10⁴ yr. We explore different velocity resolutions and beam sizes. The synthetic observations are computed by convolving the full-resolution images obtained from the radiative transfer simulations with a circular Gaussian beam. In the¹² CO emission, there is no evidence of kinks in any channels. On the other hand, a local perturbation nearby the planet is visible in the ¹³CO emission. Figure 11 compares the simulated ¹³CO emission at three different angular resolutions: 0.03″, 0.02″, and 0.01″. We show the channels at v = 5.0km s⁻¹ and v = 7.5 km s⁻¹. These velocities, equidistant from the systemic velocity, are just outside the window where the cloud totally absorbs ¹³CO emission. We only show the results with a channel width of 0.5 km s⁻¹ since a higher velocity resolution does not lead to any major change. The planet is located at the azimuthal angle ϕ = 145° with respect to the disk major axis (see Fig. 11). We note that a velocity perturbation near the planet is visible at a resolution of 0.02″ in the 5.0 km s⁻¹ channel, becoming more evident at 0.01″. Instead, no kink is detected in the 7.5 km s⁻¹ channel. Therefore, this analysis proves that a planet-induced kink would be visible only in ¹³CO emission and at a very high resolution of at least 0.02″, assuming that the planet is located in a favorable position that allows the perturbation to be detected in a channel not obscured by the cloud. The sensitivity needed to clearly detect ¹³CO emission at such a high angular resolution leads to an inaccessible observing time (about several tens of days) with current ALMA capabilities.

We identify two main reasons why the planet-induced perturbation is only detected in ¹³CO and not in ¹²CO in our simulations. First, the emission of ¹²CO comes from regions at high altitudes, whose dynamics could be less affected by the planet orbiting on the midplane (Rabago & Zhu 2021), whereas ¹³CO is emitted closer to the midplane. Second, the spatial resolution of our models naturally becomes coarser at higher altitudes, where ¹²CO is emitted, thus making it more challenging to detect localized perturbations induced by a planet. So far, published detections of planet-induced kinematical kinks in real ALMA observations of disks have involved only ¹²CO (Pinte et al. 2018, 2019, 2020). The explanation is that ¹²CO has a bright line emission, allowing for very high signal-to-noise ratios even at high spectral and spatial resolution. Kinks in ¹³CO have a larger kinematical signal, but are more difficult to detect because of its fainter line brightness.

4.4 Minimum Planet Mass

The comparison in Fig. 5 shows that final planet masses of 0.9 M_Jup and 1.4 M_Jup are able to open a gap in the dust, though producing a ring that is shifted with respect to the observed one. While all planets considered in Fig. 5 have an equal initial distance d_P,₀ from the central star of 10.0 au, now, we aim to better constrain the minimum planet mass able to explain the observed dust ring by modifying the parameter d_P,0. Using the same initial planet masses that eventually resulted in 0.9 M_Jup and 1.4 M_Jup, that is, M_P,0 = 0.5 and 1.0 M_Jup, respectively, we explore different values of d_P,0. The best cases are presented in Fig. 12, after an evolution time of 4 × 10⁴ yr.

The planet with M_P = 0.9 M_Jup and d_P = 11.0 au is not able to reproduce the observed profiles. As explained in Sect. 2.1.3, our numerical resolution prevents us from trusting the dust morphology inside the first few astronomical units of the disk. In this case, therefore, we do not deem the apparent presence of an inner ring within 4 au significant (which might be reminiscent of the observed inner disk in CIDA 1), but we clearly see an excess emission between ~4 and ~ 10 au (a region properly mapped by SPH particles in our simulation), which is significantly higher than what is observed. In addition, the morphology of the external dust ring is slightly different, particularly in Band 4. Conversely, we note an excellent match between the observed dust ring and the one obtained from a planet with M_P = 1.4 M_Jup and d_P = 10.3 au in both Band 7 and Band 4. Therefore, we conclude that a minimum mass of ~1.4 M_Jup can explain the observed dust ring.

Such constraints on the planet mass have been obtained by varying the initial planet mass and location but always fixing the initial value of the disk viscosity at α_SS = 5 × 10⁻³. Through a set of hydrodynamical simulations of disks with a planet, Facchini et al. (2018) prove that the dust gap width (which they define as the distance between the peak of the dust ring and the center of the gas gap) remains substantially unchanged when varying the disk viscosity by an order of magnitude (from α_SS = 10⁻³ to 10⁻⁴), thus showing its effective independence from viscosity (see Fig. 11 in their paper). On the contrary, the CO gap width appears to be dependent on viscosity (Fig. 13 in Facchini et al. 2018)

Beyond viscosity, in our simulations we maintained the same initial values of all the other parameters characterizing the disk (see Table 1) and the same dust composition. Given the computational cost of each simulation, a full exploration of all these parameters is out of the scope of this work.

Fig. 11 Synthetic images of the 13CO (J = 3−2) emission from our simulation with MP = 2.0 MJup after an evolution time of 4 × 104 yr. We compute the images at two different channels: v = 5.0 km s−1 (top panels) and v = 7.5km s−1 (bottom panels), with a velocity resolution of 0.5km s−1. We use three different circular beam sizes: 0.03″, 0.02″, and 0.01″ (FWHMs, from left to right), indicated by the white circles in the lower-left cornel of each plot. The blue dot represents the planet position at the azimuthal angle ϕ = 145° (relative to the angle ϕ = 0°, which is located at the disk major axis indicated by the dashed white line in the top-left panel).

Fig. 12 As in Fig. 5, but for the cases of simulations with final planet masses 0.9 and 1.4 MJup with different initial distances from the star. For the simulation in the first row, MP,0 = 0.5 MJup and dP,0 = 12.0 au, and, after an evolution time of 4 × 104 yr, MP = 0.9MJup and dP = 11.0 au. Second row: MP,0 = 1.0 MJup and dP,0 = 11.0 au, and, after the same evolution time, MP = 1.4 MJup and dP = 10.3 au.

Fig. 13 As in Fig. 5, but for the case of a simulation with higher planet mass. The planet started with MP,0 = 4.0 MJup and dP,0 = 9.0 au, and, after an evolution time of 4 × 104 yr, it reached MP = 4.5 MJup and dP = 8.5 au.

4.5 Maximum Planet Mass

The simulations reported in Fig. 5 do not allow an estimation of the maximum planet mass able to recover the dust ring in the observations. The cases with the highest final planet masses, 2.0 and 2.5 M_Jup, starting at a distance of 10.0 au from the central star, both produce an accurate match with the observed dust emission morphology. We try to constrain the maximum planet mass by simulating a new case with a higher initial planet mass of M_P,0 = 4.0M_Jup. Choosing an initial location d_P,0 = 9.0 au, after 4 × 10⁴ yr the planet reached d_P = 8.5 au with M_P = 4.5 M_Jup, and also this model reproduces the observed dust ring in both Band 7 and Band 4 (see Fig. 13).

Thus, we followed another approach to give an estimate of the possible maximum planet mass. Lodato et al. (2019) report an empirical relation according to which the width of an observed dust gap scales with the planet Hill radius (see Eq. (4)) in case of low disk viscosity (α_SS ≲ 0.01; Dodson-Robinson & Salyk 2011; Pinilla et al. 2012a; Rosotti et al. 2016; Fung & Chiang 2016; Facchini et al. 2018). Averaging results from hydrodynamical simulations (Clarke et al. 2018; Liu et al. 2019), they derive

$$\Delta =5.5R_{\text{H}},$$(6)

where Δ is defined as the distance between the minimum intensity in the gap and the peak intensity in the ring. Hence, we take the value of 2Δ as the total dust gap width. This relation assumes that the planet position coincides with the location of the dust gap, and that the gap is opened by a single planet. With such an analytical approach, we do not suffer from the numerical limitations of the hydrodynamical simulations (see Sect. 2.1.3). Instead, we can place an upper limit to the planet mass by requiring the dust gap that it produces not to be so wide as to prevent the formation of the observed inner disk (which we did not reproduce in the hydrodynamical simulations).

From the intensity radial profile of CIDA 1 observed in both Band 7 and Band 4 (see Fig. 5), the total dust gap width is 19.2 au. We measured this distance from the peak of the inner disk emission (i.e., the center of the disk) to the peak of the dust ring. We assumed that the observed inner disk is not a fast transient feature, but a long-lived structure. We fixed the star mass at 0.2 M_⊙ and adopted a planet location between 9.5 au, consistent with the values in Table 2, and 7.5 au, assuming that a higher-mass planet needs a shorter distance from the star to correctly reproduce the observed dust ring, as occurred with the case of the 4.5 M_Jup planet (Fig. (13)). Employing Eq. (6) with these values, we obtained the total dust gap width 2Δ versus the planet mass, as depicted in Fig. 14. Comparing the computed gap widths with the observed one, we find that the maximum planet mass able to create the CIDA 1 gap should be ~4−8 M_Jup. A more massive planet would carve a wider gap, either generating the dust ring at a greater distance from the star or depleting the inner disk.

4.6 Comparison with Previous Estimates of the Planet Mass

Assuming a stellar mass of 0.1 M_⊙, a disk viscosity α_SS= 10⁻³ and that the gap is carved by a single planet located at the minimum of the gap (~8 au), Pinilla et al. (2021) employed the Crida et al. (2006) criterion and estimated that the minimum planet mass to open a gap in the gas is ~0.4 M_Jup. However, such a planet would generate a pressure maximum at ~13 au. So, either a more massive planet or a higher distance from the star is needed to retrieve the observed dust ring at ~ 20 au. Following this argument, the authors computed dust evolution simulations considering two cases: a 2. 4 MJup planet located at 8 au and a 0. 5 MJup planet at 12 au, assuming a gap shape obtained analytically. Their 1D models evolved for 10⁶ yr, starting from micron-sized dust particles and then considering dust growth, fragmentation and erosion. Both planets were able to form a dust ring at the observed location. The 2. 4 MJup planet case was excluded: the formed gap was too deep to allow replenishment of dust from the outer region to an inner disk, resulting in an empty internal cavity. The 0. 5 MJup planet case, instead, was able to reproduce the contrast between the inner disk and the ring after 10⁵ yr of evolution. At longer times, however, the inner disk became fainter.

The models presented in this work differ significantly from those in Pinilla et al. (2021). Our hydrodynamical simulations take into account all the complex dynamical interactions between the disk gas and dust components and the planet in 3D. The computational effort needed inevitably leads to evolution times shorter than 10⁶ yr. We did not include the effects of dust evolution such as grain growth or fragmentation, but we assumed that dust growth has already occurred, simulating the dynamical behavior of different dust populations with fixed grain sizes, ranging from submicron to millimeter scales. In our models, most of the dust mass is contained in the bigger grain sizes, allowing us to reproduce the observed low value of the spectral index. On the other hand, the simulations of Pinilla et al. (2021) struggled to produce large dust grains, leading to higher values of the spatially-integrated spectral index (~2.6). The numerical limitations described in Sect. 2.1.3 prevent us from reproducing the inner disk, but the final synthetic images in our models allow, nonetheless, a careful match with the observed dust ring and the gas emission morphology. Moreover, thanks to the comparison between the simulated and observed gas channel maps, we find that a better estimate of the stellar mass is 0.2 M_⊙, whereas Pinilla et al. (2021) assumed a 0.1 M_⊙ central star. This difference in the stellar mass, along with a disk viscosity α_SS = 10⁻³, lower than the initial value of 5 × 10⁻³ adopted in our simulations, might explain why Pinilla et al. (2021) predict a planet mass of ~ 0. 5 MJup while our models require a minimum planet mass of ~1.4 M_Jup to reproduce the observed dust ring.

The two modeling approaches are different, and neither is fully complete. Future studies need to take into account the phenomena of dust growth and fragmentation within comprehensive hydrodynamical and radiative transfer simulations. Including all of these effects, each with its typical timescale, may end up reducing the uncertainty on the mass of the planet that can generate the observed substructures in CIDA 1. Moreover, it is important to note that our work and the one by Pinilla et al. (2021) aim at reproducing the observed disk morphology with a single planet. It remains to be explored whether a multi-planet system could reproduce the disk emission, and future simulations including more than one embedded planet are needed to test such a scenario.

Fig. 14 Total gap width in the dust as a function of the planet mass, derived from Eq. (6). The red-shaded area indicates the maximum planet mass that could fit the observed total gap width, assuming that the planet is located between 7.5 and 9.5 au from the star.

4.7 Implications for Planet Formation Around VLM Stars

Our analysis of the CIDA 1 system is consistent with the presence of a planet more massive than Jupiter around a 0. 2 M_⊙ star. Even considering that CIDA 1 may be an outlier in terms of the initial disk mass or its properties, it can still provide relevant constraints on planet formation theories. It is thus important to compare what we find in this system with the prediction of the leading theories.

Core accretion model (Pollack et al. 1996) is based on the dynamics of dust in disks and how it evolves from submicron dust grains in the interstellar medium to kilometer-sized planetesimals through collisions. During this process, when dust grains reach millimeter size, they should rapidly drift toward the central star, preventing grains from growing into planetesimals and planetary cores (Weidenschilling 1977). This barrier for planet formation is even harder to overcome for disks around VLM stars since dust radial drift is more efficient in these environments (Pinilla et al. 2013; Zhu et al. 2018). Furthermore, the disk mass content is crucial to determine whether planet formation may occur. Disk population studies in nearby star-forming regions (e.g., Barenfeld et al. 2016; Ansdell et al. 2016; Testi et al. 2016; Pascucci et al. 2016; Cieza et al. 2019; Williams et al. 2019; Sanchis et al. 2020; Villenave et al. 2021) prove that the scaling relation between disk mass and stellar mass is steeper than linear. Therefore, it appears that disks in the low-stellar-mass regime do not possess enough material to form the known planetary systems (Manara et al. 2018), unless the planet-formation efficiency is 100% (Mulders et al. 2021). However, these studies consider disks with a mean age >10⁶yr, whereas planet formation via gravitational instabilities should occur at earlier times (Kratter & Lodato 2016). The properties of disks in the early Class 0/I protostars is an active area of research and the properties of such disks are highly uncertain and debated at the moment, both on observational and modeling grounds (Maury et al. 2019; Tychoniec et al. 2020; Lebreuilly et al. 2021). Nevertheless, it seems likely that at least some young disks will be prone to undergo gravitational instabilities and possibly form massive planets early on through this path.

Various theoretical studies employing numerical simulations have been performed to investigate planet formation through core accretion or gravitational instability in the low-stellar-mass regime. Payne & Lodato (2007) assessed planet formation around brown dwarfs adapting models for higher stellar masses based on core accretion. Through Monte Carlo simulations, they found that Earth-like planets can form in this condition and the planet mass depends strongly on the disk mass. However, none of their simulations showed a planetary rocky core accreting a gaseous envelope to form a giant planet. A way to overcome the radial drift barrier is a rapid rocky core growth. Pebble accretion is a mechanism able to speed up significantly the giant planet formation process (i.e., Lambrechts & Johansen 2012; Bitsch et al. 2015). Liu et al. (2020) carried out a theoretical study on planet formation driven by pebble accretion in the (sub)stellar mass range between 0.01 and 0.1 M_⊙. First, they calculated the initial masses of protoplanets by extrapolating previous numerical simulations conducted in previous literature. Next, they performed a population synthesis study to track the growth and migration of a large sample of protoplanets under the influence of pebble accretion. Their results show that, around a 0.01 M_⊙ brown dwarf, planets can grow up to 0.1−0.2 M_⊕, while, around 0.1 M_⊙ stars, planets can reach a maximum mass of 2−3 M_⊕. Findings from this study show that even pebble accretion does not seem to be sufficient to form gas giants around VLM stars and brown dwarfs. Miguel et al. (2020) used a population synthesis approach based on planetesimal accretion to explore planet formation in the stellar mass range between 0.05 and 0.25 M_⊙. They let the synthetic population of planetary systems evolve for 10⁸ yr. The authors find that to form planets with masses higher than 0.1 M_⊕ they need stars of at least 0.07 M_⊙, implying that planet formation around brown dwarf may not be a usual outcome. Then, stars with masses higher than 0.15 M_⊙ are necessary to form planets more massive than the Earth. Therefore, from all of these studies, we conclude that core accretion model currently cannot explain the presence of gas giants around VLM stars or brown dwarfs. Either the core accretion theory is incomplete, or another mechanism for planet formation is needed. This is confirmed by Lodato et al. (2005), who discussed the origin of the 5 M_Jup planet detected around the 25 M_Jup brown dwarf 2MASSW J1207334-393254 (Chauvin et al. 2005). They found that the core accretion mechanism is far too slow to generate such a planet in less than 10⁷ yr, the estimated age of the system. Therefore, the authors proposed gravitational instabilities arising during the early phases of the disk lifetime as a viable possibility for the formation of the planet.

Haworth et al. (2020) assessed the susceptibility for disk fragmentation due to self-gravity depending on the properties of the disk and the host star. They performed a set of SPH simulations and accounted for the stellar irradiation. Their results show that disks around VLM stars are less prone to fragment with respect to disks around more massive stars. Fragmentation can occur around VLM stars only if the disk is optically thick to stellar irradiation and there is a high disk-to-star mass ratio (M_disk/M_* ≳ 0.3). Mercer & Stamatellos (2020) focused on the role of disk instability for planet formation around low-mass stars. They conducted a set of SPH simulations to study the fragmentation conditions of initially gravitationally stable disks, whose mass was slowly but steadily incremented. They considered stars with masses between 0.2 and 0.4 M_⊙. The authors find that, via gravitational instabilities, protoplanets form very fast, within a few thousand years, provided that the disk is massive enough (M_disk/M_* ≳ 0.3). The formed planets have masses between 2 and 6 M_Jup and they are initially distant from the central star, between 15 and 105 au, but may migrate afterward. Hence, a plausible origin of a giant planet around CIDA 1 might be the fast fragmentation of the disk due to gravitational instabilities, assuming the disk was massive enough in its early stages.

5 Conclusions

In this work, we inspect the origin of the substructures observed by ALMA in the protoplanetary disk around CIDA 1, one of the very few young VLM stars with high-resolution data. Assuming the presence of an embedded planet, we self-consistently reproduced the observed continuum emission from the dust ring and the gas line emission morphology with a single model. We performed 3D hydrodynamical simulations of gas and multi-grain dust, considering initial planet masses from 0.1 to 4.0 M_Jup and starting at a radial distance of between 9 and 11 au from the central star. Then, we computed the dust and gas emission using radiative transfer modeling. Finally, we treated the obtained images as actual ALMA observations to acquire the final synthetic images. Here we summarize our main findings:

• We compared the observed and simulated dust emission profiles in Band 7 (0.9 mm) and Band 4 (2.1 mm). We find that the observed dust ring can be explained by the presence of a giant planet with a minimum mass of ~1.4 M_Jup located at a distance of 10.3 au from the central star. The formation of such a massive planet around a VLM star such as CIDA 1 challenges our current theoretical models based on core accretion theory. A valid explanation for its origin may involve gravitational instabilities in the early stages of disk evolution.

• our numerical models cannot directly constrain the maximum planet mass consistent with observations. Thus, we used an empirical relation between the observed dust gap width and the Hill radius of the planet and estimate a maximum planet mass of ~4–8 M_Jup.

• Assuming that dust is composed of porous grains with a mixture of astronomical silicates, carbonaceous materials, and water ice, we match the observed flux with a total dust mass in our models of ~7 × 10⁻⁶ M_⊙, corresponding to a dust-to-gas mass ratio of ~10⁻².

• The observed ¹²Co and ¹³Co channel maps are strongly affected by cloud absorption. Nonetheless, our model with a final planet mass of 2.0 M_Jup can recreate the spatial extent of the gas emission in the channels where it was detected, except for the channels that probe the highest velocities in the innermost regions, which are not considered in our model. Our simulation include a dust-to-gas mass ratio of ~10⁻², a ¹²CO abundance of [¹²CO/H₂] = 5 × 10⁻⁵, and an isotopologue abundance ratio [¹²CO/¹³CO] = 70. Since we can recreate the spatial extent and rotation pattern of the disk, we conclude that the assumed stellar mass of 0.2 M_⊙ and systemic velocity of 6.25 km s⁻¹ are solid estimates for the actual properties of CIDA 1.

• In the case of our simulation with a final planet mass of 2.0 M_Jup, we are able to reproduce the low spectral index (~2) observed at the location of the dust ring (10−30 au). The simulated dust emission in both Band 7 and Band 4 is optically thick in the range 14−20 au, with a total fraction of optically thick emission corresponding to ~40−50%. However, the emission from our model is optically thin between 20 and 30 au. The low spectral index in this region is expected in our simulations that contain millimeter-sized dust particles, but the match with the observation hints at the occurrence of grain growth and migration in CIDA 1.

• The 2.0 M_Jup planet creates a velocity perturbation that is predicted to be extremely arduous to observe. It would require ¹³CO channel maps at a very high angular resolution of at least 0.02″, with a velocity resolution of 0.5 km s⁻¹.

Our results suggest that the observed structures in CIDA 1 can indeed be explained by the presence of a massive embedded planet. It remains to be understood whether CIDA 1 constitutes a unique case or if other disks around VLM stars may also have undergone giant planet formation.

Acknowledgements

We thank the anonymous referee for the constructive comments that helped improve this paper. This work was partly supported by the Italian Ministero dell Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52I13000140001), by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) – Ref no. 325594231 FOR 2634/1 TE 1024/1-1, and by the DFG cluster of excellence Origins (www.origins-cluster.de). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00536.S, ADS/JAO.ALMA#2015.1.00934.S, and ADS/JAO.ALMA#2016.1.01511.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. M.T. has been supported by the UK Science and Technology research Council (STFC) via the consolidated grant ST/S000623/1, and by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project).

Appendix A Sound Speed Power Law Exponent from the SED

The profile of the sound speed in the SPH simulations follows a power law (Eq. 2). It is possible to derive the value of its exponent from the disk spectral energy distribution (SED), assuming that the disk emits as a multicolor blackbody (Beck-with et al. 1990). In the frequency range k_B T_out ≪ hv ≪ k_B T_in (where k_B is the Boltzmann constant, T_in and T_out are the disk temperatures at the inner and outer radii of the disk, respectively), usually corresponding to a wavelength range roughly between 2 µm and 100 µm, the slope of the SED is connected to the exponent q of the disk sound speed profile by this relation:

$$\frac{{\text{dlog}}_{10}\left(\lambda F_{\lambda }\right)}{{\text{dlog}}_{10}\left(\lambda \right)}=-4+\frac{1}{q}.$$(A.1)

The measured SED of CIDA 1 is reported in Hendler et al. (2017). We performed a linear fit of the value of log₁₀(λ F_λ) with respect to log₁₀(λ) in the wavelength range 2.5 − 100 µm (Fig. A.1), obtaining a slope of ≈ −0.62. Using the (A.1), we derive that q ≈ 0.3.

Fig. A.1 SED of CIDA 1. Data are from Hendler et al. (2017), and a typical uncertainty of 20% is assumed for each flux measurement. The orange line represents the linear fit performed on the data in the wavelength range 2.5 − 100 µm, between the two vertical dashed gray lines.

Appendix B Planet Accretion and Migration

In SPH simulations with Phantom, the planet accretion of gas and dust contained in the disk is self-consistently computed, and the planet migration is directly calculated from the disk gravitational torque, without the need for any prescriptions (Price et al. 2018b). In Fig. B.1, we report how planet-star distance and planet mass change with time, in the case of the simulation with an initial planet mass of 1.5 M_Jup. In both plots, a fast transient is evident in the first 10⁴ yrs, corresponding to the phase of gap opening. After that, the planet stabilizes in a slow type II migration with a smaller increase in mass. Therefore, we can notice that at 4 × 10⁴ yrs, the evolution time at which most of our analyses are performed, the simulation has overcome the initial transient and is in a slowly evolving phase.

From the upper plot in Fig. B.1, a spread in the planet-star distance is visible, indicating an elliptical orbit of the planet in our simulation. Knowing the periastron r_p and apoastron r_a, we can calculate the planet orbital eccentricity:

$$e=\frac{r_{\text{a}}-r_{\text{p}}}{r_{\text{a}}+r_{\text{p}}},$$(B.1)

obtaining a value e ≈ 0.01 after 8 × 10⁴ yrs.

Fig. B.1 Time evolution of planet-star distance (top panel) and planet mass (bottom panel) from the simulation with an initial planet mass of 1.5 MJup.

Appendix C Dust Mass Rescaling

For obtaining the final synthetic dust continuum images from our simulations, we rescale the total dust masses in our models by a constant factor to match the observed flux. In our simulations, the dust-to-gas mass ratio always remains ≪ 1. This condition allows this method to be applied without altering the disk dynamics since the dust back-reaction onto the gas remains insignificant. As explained in Sect. 2.1.3, in our modeling, we focus on reproducing the observed external dust ring. For this reason, we take as reference flux the one emitted only by the external ring in the ALMA Band 7 observation. In particular, we consider the area contained in the two green dashed ellipses depicted in the top panel of Fig. C.1. The outer ellipse is positioned in the image center, has a semimajor axis of 0.27″, a semiminor axis of 0.22″, and an inclination of 11°; for the inner ellipse, we only change the semiaxis length to 0.13″ and 0.10″. The flux contained in this elliptical annulus is 29.0mJy. In our synthetic images, without an inner disk, we consider the flux contained within the red dashed ellipse in the bottom panel of Fig. C.1, which has the same parameters as the outer ellipse used for the observation image. We accept the total dust masses that lead to fluxes from our synthetic images from 28.5 to 29.5 mJy.

Fig. C.1 Visual representation of the areas considered for the flux calculation to rescale the total dust masses in our simulations. We take as a reference the flux from the ALMA observation in Band 7, which is contained in the region between the two dashed green ellipses in the top panel. Then, we rescale the total dust mass in our simulations so that, from the synthetic images in Band 7, the flux within the dashed red ellipse shown in the bottom panel matches the reference flux from the observation. In the bottom panel, we display, as an example, the image obtained from the simulation where the planet reaches a mass of 2.0 MJup after 4 × 104 yrs.

Appendix D Dust Emission Profiles at Longer Evolution Times

Fig. D.1 shows the dust emission radial profiles in Band 7 and Band 4 of each simulation after an evolution time of 8 × 10⁴ yrs. Compared to the results after 4 × 10⁴ yrs (Fig. 5), dust grains had more time to undergo radial drift and possibly accrete onto the star.

We still note that a final planet mass of at least 2.0 M_Jup best reproduces the observations. The slight intensity excess just outside the dust ring has been removed by radial drift. The small mismatch in the peak flux for the simulation with M_P = 2.0 M_Jup and M_P = 2.6 M_Jup in Band 4 is due to the fact that, at this time, the bigger millimeter-sized grains with faster radial drift accreted more onto the central star, thus reducing the continuum emission at this longer wavelength.

Fig. D.1 As in Fig. 5, but for the case of simulations evolved for 8 × 104 yrs.

Appendix E Cloud Absorption

We aim here to qualitatively estimate the cloud absorption affecting the observed gas channel maps. Davis et al. (2010) mapped the ¹²CO emission at the transition J = 3 − 2 in the region of the Taurus molecular cloud, using the Heterodyne Array Receiver Programme (HARP) on the James Clerk Maxwell Telescope. Figure E.1 shows the ¹²CO map of Taurus that they obtained. The authors divided the region into five different subregions, and in a couple of them, namely L1595 E and L1495 SE, they also collected the ¹³CO J = 3 − 2 emission. CIDA 1 is located in the subregion L1495 W. Since in this area the ¹³CO data are missing, we consider the nearby subregion L1595 E for our estimates.

Figure E.2 reports the ¹²CO and ¹³CO spectra of the antenna temperature $\left(T_A^{*}\right)$. As stated in Davis et al. (2010), assuming the abundance ratio [¹²CO/¹³CO] = 70, the ¹²CO optical depth τ₁₂ can be calculated as:

$$\frac{{\displaystyle \int T_{\text{A}}^{*}\left({}^{12}\text{CO}\right)\text{d}v}}{{\displaystyle \int T_{\text{A}}^{*}\left({}^{13}\text{CO}\right)\text{d}v}}=\frac{1-{\text{e}}^{-{\tau }_{12}}}{1-{\text{e}}^{-{\tau }_{13}}},$$(E.1)

where τ₁₃ = τ₁₂/70 is the ¹³CO optical depth. We integrate the antenna temperatures along velocity intervals with the same centers and width of the observed ¹²CO channel maps (top panels in Fig. 6). The resulting ¹²CO optical depth is presented in Fig. E.3. The morphology of the optical depth spectrum follows approximately the cloud absorption behavior on the¹² CO observation, being strongest at ~ 6.0 − 7.0 km s⁻¹ while fading at ~ 4.0 km s⁻¹ and ~ 8.5 km s⁻¹, where there is a better match between the intensity profiles from our model and the ones from the observation (Fig. 7). Nonetheless, it must be stressed that our estimates on the ¹²CO spectrum are based on data from a region nearby the one hosting CIDA 1 and, furthermore, these data refer to the total emission of the cloud along the line of sight, while the position of CIDA 1 in this direction, and so the gas column density between us and the source, is unknown.

Fig. E.1 Map of the Taurus molecular cloud obtained by the James Clerk Maxwell Telescope using the 12CO (J = 3−2) emission (integrated between 2.0 and 11.0 km s−1). Areas also mapped with 13CO (J = 3−2) emission (integrated between 4.9 and 9.1 km s−1) are indicated by rectangular boxes. We calculated the optical depth property in the subregion L1495 E, whereas CIDA 1 is located in the region L1495 W. The image is adapted from right-hand panel of Fig. 1 in Davis et al. (2010).

Fig. E.2 12CO (J = 3−2) and 13CO (J = 3−2) spectra with respect to the local standard of rest (LSR) in the Taurus subregion L1495 E. The spectra are averaged over the entire extent of the subregion. Data are from Davis et al. (2010).

Fig. E.3 12CO (J = 3−2) optical depth in the Taurus subregion L1495 E, calculated from the spectra in Fig E.2 using Eq. E.1. We focus only on the velocity interval where line emission is detected. Bins replicate the central velocity and resolution of the observed 12CO channel maps.

References

All Tables

All Figures

	Fig. 1 ALMA observation in Band 7 (0.9 mm) of the system CIDA 1 (Pinilla et al. 2021), with a full width at half maximum (FWHM) beam size of 0.050″ × 0.034″, indicated by the white ellipse in the bottom-left corner.
In the text

Fig. 2 Face-on gas and dust surface density, in units of g cm−2, from the Phantom simulation with an initial planet mass of 1.5 MJup and a final planet mass of 2.0 MJup after ≈4 × 104 yr (≈560 planet orbits). Gas density is displayed in the top-left corner, and the remaining panels present the surface density for the 11 simulated dust grain sizes. The bigger central green dot represents the star, and the smaller green dot indicates the planet position. Images were produced using the tool Splash (Price 2007).

In the text

Fig. 3 Refinement of the hydrodynamical model for the radiative transfer simulation via the extrapolation and interpolation of the dust densities for a single SPH particle. In the left panel, starting from dust density values obtained with SPH simulations (red dots), we perform a linear fit for log ρd values associated with the two smallest simulated grain sizes to obtain densities for grains down to 5 nm (blue line). Then, we use a cubic spline interpolation on the remaining log ρd values (orange line). In the right panel, with the linear fit, we obtain ten dust densities related to grain sizes between 5 nm and 0.37 µm (to the left of the vertical dashed gray line), while with the cubic spline interpolation we obtain 50 densities for grain sizes between 0.37 µm and 1 mm (to the right of the vertical dashed gray line).

In the text

Fig. 4 Band 7 (0.9 mm, top panels) and Band 4 (2.1 mm, bottom panels) comparison between dust continuum images from the ALMA observation and our synthetic observations obtained from the simulation with a final planet mass of 2.0 MJup after 4 × 104 yr, along with their residuals. The synthesized beam (0.050″ × 0.034″ for Band 7, 0.048″ × 0.032″ for Band 4, FWHMs) is shown in the bottom-left corner of each image. The contours in the residuals represent the 5σ and −5σ flux levels, respectively, for Band 7 and Band 4.

In the text

Fig. 5 Comparison between radial profiles of azimuthalfy averaged intensities from the observations, simulations, and corresponding residuals, after the simulations have evolved for 4 × 104 yr. The left-hand panels show results in Band 7, while the right-hand panels refer to Band 4. Each row contains the intensity profiles from the same simulation with a specific final planet mass, ranging from 0.2 MJup (top row) to 2.5 MJup (bottom row). The distance of the planet from the central star is indicated by the vertical dash-dotted brown line. The mean FWHMs of the synthesized beam in Band 7 and Band 4 are shown as horizontal bars in the top row panels. The radial intensity profiles are obtained by deprojecting the disk, using the disk inclination angle i = 37.5°, and dividing it into overlapping annuli, whose widths correspond to the averaged synthesized beam diameter. For each bin, uncertainties are calculated by taking the intensity standard deviation for each pixel in the annulus and dividing it by the square root of the number of synthesized beams contained in the area of the annulus.

In the text

Fig. 6 Comparison between the ALMA observation (top panels) and the synthetic images (bottom panels) of CIDA 1 12CO (J = 3−2) channel maps. Color scales are the same in the top and bottom panels. Synthetic images are obtained from the simulation with a final planet mass of 2.0 MJup after an evolution time of 4 × 104 yr. We fixed the systemic velocity at 6.25 km s−1. The contour level traces the 15σ emission from the Band 7 continuum image. The velocity resolution is 0.5 km s−1, and the central velocities of each channel are indicated on the top of each panel. The synthesized beam (0.101″ × 0.075″, FWHM) is shown in the lower-left corner of each panel.

In the text

Fig. 7 Contour levels (left panel in each subplot) at 5, 10, 15mJy beam−1, and intensity profile comparisons (right panel in each subplot) between the 12CO observation and simulation with a 2.0 MJup planet mass after 4 × 104 yr, in the six channels not completely obscured by the cloud (v = 4.0, 4.5, 5.0, 8.0, 8.5, 9.0 km s−1). Intensity profiles are taken along the PA = 11° of the disk, and the direction is indicated by the dotted gray line in the contour levels plots. The uncertainty in the observation intensity profile is the 1σ noise level (1.6 mJy beam−1 ).

In the text

Fig. 8 Comparison between the ALMA observation (top panels) and the synthetic images (bottom panels) of CIDA 1 13CO (J = 3-2) channel maps. Color scales are the same in the top and bottom panels. Synthetic images are obtained from the simulation with a final planet mass of 2.0 MJup after an evolution time of 4 × 104 yr. We fixed the systemic velocity at 6.25 km s−1. The contour level traces the 15σ emission from the Band 7 continuum image. The velocity resolution is 1 km s−1, and the central velocities of each channel are indicated on the top of each panel. The synthesized beam (0.104″ × 0.076″, FWHM) is shown in the lower-left corner of each panel.

In the text

Fig. 9 Contour levels (left panel in each subplot) at 3, 6, 9 mJy beam−1, and intensity profile comparisons (right panel in each subplot) between the 13CO observation and simulation with a 2.0MJup planet mass after 4 × 104 yr, in the four channels not completely obscured by the cloud (v = 4.5, 5.5, 7.5, 8.5 km s−1). Intensity profiles are taken along the PA = 11° of the disk, and the direction is indicated by the dotted gray line in the contour levels plots. The uncertainty in the observation intensity profile is the 1σ noise level (1.2 mly beam−1).

In the text

Fig. 10 Optical depth in our models and comparison between observed and simulated spectral indices. Left and central panels: Azimuthally averaged radial profiles of the optical depth along the line of sight and the intensities from the full-resolution dust emission image of the simulation (i.e., not convoluted with the synthesized beam) with a final planet mass of 2.0 MJup after 4 × 104 yr, in Band 7 and Band 4. The vertical dash-dotted brown line indicates the location of the planet, and the horizontal dotted gray line marks the optical depth τ = 1 level. The region where the emission is optically thick (optical depth τ ≥ 1) is shown in beige. Right panel: azimuthally averaged radial profiles of the spectral index from the ALMA observation and our simulation with a final planet mass of 2.0 MJup after 4 × 104 yr. Values are obtained using Eq. (5) and the synthetic image intensity profiles in Fig. 5. Uncertainties include the original error in the intensity profiles and a 10% flux error due to calibration. As a reference, the distance of the planet from the central star is indicated by the vertical dash-dotted brown line, and the position of the peak intensity is displayed by the vertical dotted blue line. The beige-shaded area marks the radii where the full-resolution simulated emission is optically thick (τ ≥ 1) in both Band 7 and Band 4.

In the text

Fig. 11 Synthetic images of the 13CO (J = 3−2) emission from our simulation with MP = 2.0 MJup after an evolution time of 4 × 104 yr. We compute the images at two different channels: v = 5.0 km s−1 (top panels) and v = 7.5km s−1 (bottom panels), with a velocity resolution of 0.5km s−1. We use three different circular beam sizes: 0.03″, 0.02″, and 0.01″ (FWHMs, from left to right), indicated by the white circles in the lower-left cornel of each plot. The blue dot represents the planet position at the azimuthal angle ϕ = 145° (relative to the angle ϕ = 0°, which is located at the disk major axis indicated by the dashed white line in the top-left panel).

In the text

Fig. 12 As in Fig. 5, but for the cases of simulations with final planet masses 0.9 and 1.4 MJup with different initial distances from the star. For the simulation in the first row, MP,0 = 0.5 MJup and dP,0 = 12.0 au, and, after an evolution time of 4 × 104 yr, MP = 0.9MJup and dP = 11.0 au. Second row: MP,0 = 1.0 MJup and dP,0 = 11.0 au, and, after the same evolution time, MP = 1.4 MJup and dP = 10.3 au.

In the text

	Fig. 13 As in Fig. 5, but for the case of a simulation with higher planet mass. The planet started with MP,0 = 4.0 MJup and dP,0 = 9.0 au, and, after an evolution time of 4 × 104 yr, it reached MP = 4.5 MJup and dP = 8.5 au.
In the text

	Fig. 14 Total gap width in the dust as a function of the planet mass, derived from Eq. (6). The red-shaded area indicates the maximum planet mass that could fit the observed total gap width, assuming that the planet is located between 7.5 and 9.5 au from the star.
In the text

	Fig. A.1 SED of CIDA 1. Data are from Hendler et al. (2017), and a typical uncertainty of 20% is assumed for each flux measurement. The orange line represents the linear fit performed on the data in the wavelength range 2.5 − 100 µm, between the two vertical dashed gray lines.
In the text

	Fig. B.1 Time evolution of planet-star distance (top panel) and planet mass (bottom panel) from the simulation with an initial planet mass of 1.5 MJup.
In the text

Fig. C.1 Visual representation of the areas considered for the flux calculation to rescale the total dust masses in our simulations. We take as a reference the flux from the ALMA observation in Band 7, which is contained in the region between the two dashed green ellipses in the top panel. Then, we rescale the total dust mass in our simulations so that, from the synthetic images in Band 7, the flux within the dashed red ellipse shown in the bottom panel matches the reference flux from the observation. In the bottom panel, we display, as an example, the image obtained from the simulation where the planet reaches a mass of 2.0 MJup after 4 × 104 yrs.

In the text

	Fig. D.1 As in Fig. 5, but for the case of simulations evolved for 8 × 104 yrs.
In the text

Fig. E.1 Map of the Taurus molecular cloud obtained by the James Clerk Maxwell Telescope using the 12CO (J = 3−2) emission (integrated between 2.0 and 11.0 km s−1). Areas also mapped with 13CO (J = 3−2) emission (integrated between 4.9 and 9.1 km s−1) are indicated by rectangular boxes. We calculated the optical depth property in the subregion L1495 E, whereas CIDA 1 is located in the region L1495 W. The image is adapted from right-hand panel of Fig. 1 in Davis et al. (2010).

In the text

	Fig. E.2 12CO (J = 3−2) and 13CO (J = 3−2) spectra with respect to the local standard of rest (LSR) in the Taurus subregion L1495 E. The spectra are averaged over the entire extent of the subregion. Data are from Davis et al. (2010).
In the text

	Fig. E.3 12CO (J = 3−2) optical depth in the Taurus subregion L1495 E, calculated from the spectra in Fig E.2 using Eq. E.1. We focus only on the velocity interval where line emission is detected. Bins replicate the central velocity and resolution of the observed 12CO channel maps.
In the text