# A giant planet shaping the disk around the very low-mass star CIDA 1

P. Curone, A. F. Izquierdo, L. Testi, G. Lodato, S. Facchini, A. Natta, P.

Pinilla, N. T. Kurtovic, C. Toci, M. Benisty, et al.

► To cite this version:

P. Curone, A. F. Izquierdo, L. Testi, G. Lodato, S. Facchini, et al. A giant planet shaping the disk around the very low-mass star CIDA 1. Astronomy and Astrophysics - A&A, 2022, 665, 10.1051/0004-6361/202142748. insu-03860311

HAL Id: insu-03860311

https://insu.hal.science/insu-03860311

Submitted on 19 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A giant planet shaping the disk around the very low-mass star CIDA 1

P. Curone<fnref n="1" />, A. F. Izquierdo<fnref n="2" />,<fnref n="3" />, L. Testi<fnref n="2" />,<fnref n="4" />, G. Lodato<fnref n="1" />, S. Facchini<fnref n="1" />,<fnref n="2" />, A. Natta<fnref n="5" />, P. Pinilla<fnref n="6" />, N. T. Kurtovic<fnref n="6" />, C. Toci<fnref n="1" />, M. Benisty<fnref n="7" />,<fnref n="8" />, M. Tazzari<fnref n="9" />, F. Borsa<fnref n="10" />, M. Lombardi<fnref n="1" />, C. F. Manara<fnref n="2" />, E. Sanchis<fnref n="11" />,<fnref n="2" />, and L. Ricci<fnref n="12" />

Received 25 November 2021 / Accepted 17 May 2022

## ABSTRACT

**Context.** Exoplanetary research has provided us with exciting discoveries of planets around very low-mass (VLM) stars ($0.08 M_{\odot} \le M_{\star} \le 0.3 M_{\odot}$; e.g., TRAPPIST-1 and Proxima Centauri). However, current theoretical models still strive to explain planet formation in these conditions and do not predict the development of giant planets. Recent high-resolution observations from the Atacama Large Millimeter/submillimeter Array (ALMA) of the disk around CIDA 1, a VLM star in Taurus, show substructures that hint at the presence of a massive planet.

**Aims.** We aim to reproduce the dust ring of CIDA 1, observed in the dust continuum emission in ALMA Band 7 (0.9 mm) and Band 4 (2.1 mm), along with its ¹²CO (J = 3–2) and ¹³CO (J = 3–2) channel maps, assuming the structures are shaped by the interaction of the disk with a massive planet. We seek to retrieve the mass and position of the putative planet, through a global simulation that assesses planet-disk interactions to quantitatively reproduce protoplanetary disk observations of both dust and gas emission in a self-consistent way.

**Methods.** Using a set of hydrodynamical simulations, we model a protoplanetary disk that hosts an embedded planet with a starting mass of between 0.1 and 4.0 $M_{\text{Jup}}$ and initially located at a distance of between 9 and 11 au from the central star. We compute the dust and gas emission using radiative transfer simulations, and, finally, we obtain the synthetic observations, treating the images as the actual ALMA observations.

**Results.** Our models indicate that a planet with a minimum mass of ~1.4 $M_{\text{Jup}}$ orbiting at a distance of ~9–10 au can explain the morphology and location of the observed dust ring in Band 7 and Band 4. We match the flux of the dust emission observation with a dust-to-gas mass ratio in the disk of ~$10^{-2}$. We are able to reproduce the low spectral index (~2) observed where the dust ring is detected, with a ~40–50% fraction of optically thick emission. Assuming a ¹²CO abundance of $5 \times 10^{-5}$ and a ¹³CO abundance 70 times lower, our synthetic images reproduce the morphology of the ¹²CO (J = 3–2) and ¹³CO (J = 3–2) observed channel maps where the cloud absorption allowed a detection. From our simulations, we estimate that a stellar mass $M_{\star} = 0.2 M_{\odot}$ and a systemic velocity $v_{\text{sys}} = 6.25 \text{ km s}^{-1}$ are needed to reproduce the gas rotation as retrieved from molecular line observations. Applying an empirical relation between planet mass and gap width in the dust, we predict a maximum planet mass of ~4–8 $M_{\text{Jup}}$.

**Conclusions.** Our results suggest the presence of a massive planet orbiting CIDA 1, thus challenging our understanding of planet formation around VLM stars.

**Key words.** protoplanetary disks – planet-disk interactions – stars: individual: CIDA 1 – planets and satellites: formation – hydrodynamics – radiative transfer

## 1. Introduction

Planets have been detected around very low-mass (VLM) stars, defined as stars with masses $\le 0.3 M_{\odot}$ but still above the hydrogen-burning limit of $\sim 0.08 M_{\odot}$, below which the brown

dwarf regime begins (Liebert & Probst 1987). Particularly fascinating are the cases of TRAPPIST-1, a ~0.085 $M_{\odot}$ 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.12 $M_{\odot}$ and hosts a terrestrial planet, a

<footnote n="1">¹ Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy
e-mail: pietro.curone@unimi.it</footnote>

<footnote n="2">² European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany</footnote>

<footnote n="3">³ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="5">⁵ School of Cosmic Physics, Dublin Institute for Advanced Studies, 31 Fitzwilliams Place, Dublin 2, Ireland</footnote>

<footnote n="6">⁶ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117, Heidelberg, Germany</footnote>

<footnote n="7">⁷ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS, UMI 3386. Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="8">⁸ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="9">⁹ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="10">¹⁰ INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate (LC), Italy</footnote>

<footnote n="11">¹¹ Freie Universität Berlin, Institute of Geological Sciences, Malteserstr. 74-100, 12249 Berlin, Germany</footnote>

<footnote n="12">¹² Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91130, USA</footnote>
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<sub>Jup</sub> planet around a ~25 M<sub>Jup</sub> young brown dwarf, and Morales et al. (2019) discovered a planet with a minimum mass of 0.46 M<sub>Jup</sub> orbiting a ~0.12 M<sub>☉</sub> star. However, giant planet formation in this low-stellar-mass 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 ~55 au at 137.5 pc, the estimated distance to CIDA 1¹),

<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

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<sub>⋆</sub> = 0.1 M<sub>☉</sub>. They found that a planet with a minimum mass of ~0.3 M<sub>Jup</sub> 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<sub>⋆</sub> = 0.5 M<sub>Jup</sub> 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

<footnote n="1">¹ Throughout the work, we assume this value as the distance to CIDA 1 for consistency with Pinilla et al. (2021). Recently, Bailer-Jones et al. (2021) used Gaia EDR3 (Gaia Collaboration 2021) and derived a distance of 133.4<sup>+0.7</sup><sub>-0.9</sub> pc. Such a distance differs by ≈3% from the value in Pinilla et al. (2021), implying a systematic difference of the same order in the spatial scales derived from our models and a discrepancy of ~6% in the estimated value of the disk mass.</footnote>
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.

## 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; Ueira 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 \propto \rho_0 a_d / \Sigma_g < 1$, where $\rho_0$ is the dust grain intrinsic density, $a_d$ the grain size, and $\Sigma_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^6$ SPH particles. The disk has a gas mass of $M_{gas} = 1.5 \times 10^{-3} M_{\odot}$ 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,

$$ \Sigma_g(R) = \Sigma_c \left(\frac{R}{R_{in}}\right)^{-p} \exp\left[-\left(\frac{R}{R_c}\right)^{2-p}\right] \left(1 - \sqrt{\frac{R_{in}}{R}}\right), \quad (1) $$

where $\Sigma_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_g c_s^2$, with the following sound speed radial profile:

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

<table><caption>**Table 1.** Initial parameters used for the simulations with PHANTOM.</caption><thead><tr><th>Parameters</th><th>Value</th></tr></thead><tbody><tr><td>$N_{SPH}$</td><td>10<sup>6</sup></td></tr><tr><td>$M_{\star}(M_{\odot})$</td><td>0.2</td></tr><tr><td>$R_{\star \text{acc}}$ (au)</td><td>1</td></tr><tr><td>$R_{\text{in}}$ (au)</td><td>1</td></tr><tr><td>$R_{\text{out}}$ (au)</td><td>100</td></tr><tr><td>$R_c$ (au)</td><td>80</td></tr><tr><td>$p$</td><td>1.5</td></tr><tr><td>$q$</td><td>0.3</td></tr><tr><td>$H(R_{\text{in}})/R_{\text{in}}$</td><td>0.08</td></tr><tr><td>$\alpha_{SS}$</td><td>5 &times; 10<sup>-3</sup></td></tr><tr><td>$\alpha_{AV}$</td><td>0.36</td></tr><tr><td>$\langle h \rangle / H$</td><td>0.14</td></tr><tr><td>$M_{\text{gas}}(M_{\odot})$</td><td>1.5 &times; 10<sup>-3</sup></td></tr><tr><td>$M_{\text{dust}}/M_{\text{gas}}$</td><td>10<sup>-2</sup></td></tr><tr><td>$a_d$ (cm)</td><td>$[a_{\text{min}} = 1.6 \times 10^{-5}, a_{\text{max}} = 0.1, N = 11]$</td></tr><tr><td>$\xi$</td><td>3.5</td></tr><tr><td>$\rho_0$ (g cm<sup>-3</sup>)</td><td>2</td></tr><tr><td>$R_{P \text{acc}}$ (au)</td><td>1/4 of Hill radius</td></tr><tr><td>$d_{P,0}$ (au)</td><td>10.0</td></tr><tr><td>$M_{P,0}$ ($M_{\text{Jup}}$)</td><td>0.1, 0.5, 1.0, 1.5, 2.0</td></tr></tbody></table>

**Notes.** Each entry is commented on in the main text.

where $\rho_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(R)}{R} = \frac{c_s(R)}{v_K(R)} = \frac{H(R_{in})}{R_{in}} \left(\frac{R}{R_{in}}\right)^{1/2-q}, \quad (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 $\langle h \rangle$ and the scale height $H$ is $\langle h \rangle / H \approx 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 $\alpha_{AV}$ introduced by Lodato & Price (2010). We set the Shakura & Sunyaev (1973) viscosity $\alpha_{SS} = 5 \times 10^{-3}$, corresponding to $\alpha_{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^{-2}$ (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 µm and 1 mm. The grain size distribution follows the power law $dN(a_d)/da_d \propto a_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 $\xi = 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 $\rho_0 = 2$ g cm$^{-3}$.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Face-on gas and dust surface density, in units of g cm⁻², from the PHANTOM simulation with an initial planet mass of 1.5 M<sub>Jup</sub> and a final planet mass of 2.0 M<sub>Jup</sub> after ≈4 × 10⁴ 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).</figcaption></figure>

### 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<sub>☉</sub> (Pinilla et al. 2021). Kurtovic et al. (2021) estimate a mass of ~0.19 M<sub>☉</sub> 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<sub>⋆</sub> = 0.2 M<sub>☉</sub>. The star accretion radius is R<sub>⋆ acc</sub> = 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<sub>P acc</sub> = 1/4 R<sub>H</sub>, where R<sub>H</sub> is the planet Hill radius, defined as

$$R_H = \left( \frac{1}{3} \frac{M_P}{M_\star} \right)^{1/3} d_P, \qquad (4)$$

where d<sub>P</sub> 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<sub>P,0</sub> = 10.0 au from the star. We varied the initial planet mass M<sub>P,0</sub>: the chosen values are 0.1, 0.5, 1.0, 1.5, and 2.0 M<sub>Jup</sub>. We let the simulations evolve until they reach a stable configuration. We analyzed the simulation results at 4 × 10<sup>4</sup> 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<sup>4</sup> 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<sub>P,0</sub> = 1.5 M<sub>Jup</sub>, evaluated after 4 × 10<sup>4</sup> 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<sub>P,0</sub> = 0.5 M<sub>Jup</sub> and d<sub>P,0</sub> = 12.0 au, the second with M<sub>P,0</sub> = 1.0 M<sub>Jup</sub> and d<sub>P,0</sub> = 11.0 au. The maximum mass of the planet is assessed in Sect. 4.5, considering a simulation with M<sub>P,0</sub> = 4.0 M<sub>Jup</sub> and d<sub>P,0</sub> = 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<sub>⋆ acc</sub> = 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 R<sub>acc</sub> rapidly causes a substantial increase in computational time due to a higher
<figure><img src="image_4.png" /><figcaption>**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 \rho_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 \rho_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).</figcaption></figure>

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 $\alpha_{SS} \propto \langle h \rangle \propto \rho^{-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.

## 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., $dN(a_d)/da_d \propto a_d^{-3.5}$, which can be rearranged as $\rho_d \propto a_d^{0.5}$, with $\rho_d$ indicating the dust mass density), we fitted the values of $log_{10} \rho_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_{10} \rho_d$ values related to the remaining simulated grain sizes. Thanks to this interpolation, we sampled 50 density values, logarithmically spaced between 0.37 µm and 1 mm. 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 µm and 50 related to grain sizes between 0.37 µm 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^{-42} \text{ g cm}^{-3}$). We randomly distributed 1/5 $N_{SPH}$ 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 \pm 0.03 L_\odot$. In our simulations, we modeled the radiation as $10^8$ photon packages emitted by the central star, which is assumed as a black body with a surface temperature $T_{\text{eff}} = 3050$ K and a radius $R_* = 1.3 R_\odot$, 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^{\circ}$ and a position angle PA = 11$^{\circ}$ for the disk.

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 $^{12}$CO and $^{13}$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$^{-1}$, 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$^{-1}$. We included freeze-out and photo-dissociation for CO following the parameterizations in Williams & Best (2014). The condition $T_{\text{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_{\text{dissoc}} = 1.3 \times 10^{21}$ cm$^{-2}$ (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 $^{12}$CO abundance relative to H$_2$ [$^{12}$CO/H$_2$] = $5 \times 10^{-5}$, and the isotopologue abundance ratio [$^{12}$CO/$^{13}$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$^{-1}$ for $^{12}$CO and 1.0 km s$^{-1}$ for $^{13}$CO, we first computed the simulated channel maps at a higher velocity resolution of 0.045 km s$^{-1}$. Then, to obtain the final spectrally convoluted synthetic channels, we averaged 11 contiguous simulated $^{12}$CO channels, centered at the same velocities of the observed channels. We repeated the same procedure for the computed $^{13}$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 `sampleImage` 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

<table><caption>**Table 2.** Resulting parameters for our simulations with different initial planet masses.</caption><thead><tr><th><i>M</i><sub>P,0</sub><br>(<i>M</i><sub>Jup</sub>)</th><th><i>t</i><sub>Sim</sub><br>(yr)</th><th><i>M</i><sub>P</sub><br>(<i>M</i><sub>Jup</sub>)</th><th><i>d</i><sub>P</sub><br>(au)</th><th><i>M</i><sub>dust</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>dust</sub>/<i>M</i><sub>gas</sub></th></tr></thead><tbody><tr><td rowspan="2">0.1</td><td>4 × 10<sup>4</sup></td><td>0.2</td><td>9.4</td><td>6.0 × 10<sup>-6</sup></td><td>6.3 × 10<sup>-3</sup></td></tr><tr><td>8 × 10<sup>4</sup></td><td>0.2</td><td>9.2</td><td>4.3 × 10<sup>-6</sup></td><td>5.6 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">0.5</td><td>4 × 10<sup>4</sup></td><td>0.9</td><td>9.3</td><td>4.8 × 10<sup>-6</sup></td><td>5.8 × 10<sup>-3</sup></td></tr><tr><td>8 × 10<sup>4</sup></td><td>0.9</td><td>9.4</td><td>3.3 × 10<sup>-6</sup></td><td>5.1 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">1.0</td><td>4 × 10<sup>4</sup></td><td>1.4</td><td>9.4</td><td>7.0 × 10<sup>-6</sup></td><td>8.7 × 10<sup>-3</sup></td></tr><tr><td>8 × 10<sup>4</sup></td><td>1.5</td><td>9.4</td><td>4.0 × 10<sup>-6</sup></td><td>6.4 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">1.5</td><td>4 × 10<sup>4</sup></td><td>2.0</td><td>9.4</td><td>7.0 × 10<sup>-6</sup></td><td>8.9 × 10<sup>-3</sup></td></tr><tr><td>8 × 10<sup>4</sup></td><td>2.0</td><td>9.5</td><td>5.4 × 10<sup>-6</sup></td><td>8.8 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">2.0</td><td>4 × 10<sup>4</sup></td><td>2.5</td><td>9.4</td><td>7.3 × 10<sup>-6</sup></td><td>9.4 × 10<sup>-3</sup></td></tr><tr><td>8 × 10<sup>4</sup></td><td>2.6</td><td>9.6</td><td>6.4 × 10<sup>-6</sup></td><td>1.1 × 10<sup>-2</sup></td></tr></tbody></table>

Notes. Planet masses are evaluated at two different times during the simulation ($t_{\text{Sim}}$), namely $4 \times 10^4$ yr and $8 \times 10^4$ yr, corresponding to about 560 and 1120 orbits of the embedded planet. $M_{P,0}$ indicates the initial planet mass, while $M_P$ and $d_P$ are the planet mass and distance from the central star at the considered $t_{\text{Sim}}$. All planets started at $d_{P,0} = 10.0$ au.

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 $\ll 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 $^{12}$CO and $^{13}$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_{\text{Sim}} = 4 \times 10^4$ yr and $t_{\text{Sim}} = 8 \times 10^4$ yr. We note little difference in planet mass ($M_P$) and distance from the central star ($d_P$) in each simulation at the two times, meaning that $4 \times 10^4$ 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_{\text{Jup}}$. Here, we note a fast transient in the first $10^4$ 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
<figure><img src="image_5.png" /><figcaption>**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 $M_{\text{Jup}}$ after $4 \times 10^4$ 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\sigma$ and $-5\sigma$ flux levels, respectively, for Band 7 and Band 4.</figcaption></figure>

(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_{\text{Jup}}$ after $4 \times 10^4$ 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\sigma$ value or below the $-5\sigma$ 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 \times 10^4$ 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_{\text{Jup}}$, does not produce any major effect on the disk morphology. The $0.9 M_{\text{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_{\text{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_{\text{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 \times 10^4$ yr, we note that these last three cases also share consistent estimates of the total dust mass needed to match the observed flux, $\sim 7 \times 10^{-6} M_{\odot}$, with a dust-to-gas mass ratio of $\sim 9 \times 10^{-3}$.

Results for the dust emission from the simulations after a longer evolution time of $8 \times 10^4$ yr are presented in Appendix D. The synthetic profiles for simulations with a minimum planet mass of $2.0 M_{\text{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 \times 10^4$ and $8 \times 10^4$ 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.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison between radial profiles of azimuthally averaged intensities from the observations, simulations, and corresponding residuals, after the simulations have evolved for $4 \times 10^4$ 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 $M_{\text{Jup}}$ (top row) to 2.5 $M_{\text{Jup}}$ (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^\circ$, 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.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Contour levels (**left panel** in each subplot) at 5, 10, 15 mJy beam⁻¹, and intensity profile comparisons (**right panel** in each subplot) between the ¹²CO observation and simulation with a 2.0 *M*<sub>Jup</sub> planet mass after 4 × 10⁴ 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⁻¹). 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⁻¹).</figcaption></figure>

Our simulation reproduces appropriately the spa-
tial extent of the channels centered at the velocities
$v = 4.0, 4.5, 5.0, 8.0, 8.5, 9.0 \text{ km s}^{-1}$. This is more evident
in Fig. 7, which shows for these channels the compari-
son between the observed and simulated ¹²CO emission in
terms of intensity contour levels at 5, 10, 15 mJy beam⁻¹
and intensity profiles along the PA = 11°. The channels
$v = 4.5, 5.0, 8.0 \text{ km s}^{-1}$ 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 \text{ km s}^{-1}$, 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 com-
puted 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 mJy beam⁻¹ and the inten-
sity 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 present-
ing 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 emis-
sion 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<sub>Jup</sub> 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).

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
<figure><img src="image_8.png" /><figcaption>**Fig. 9.** Contour levels (left panel in each subplot) at 3, 6, 9 mJy beam⁻¹, and intensity profile comparisons (right panel in each subplot) between the ¹³CO observation and simulation with a 2.0 M<sub>Jup</sub> planet mass after 4 × 10⁴ yr, in the four channels not completely obscured by the cloud (v = 4.5, 5.5, 7.5, 8.5 km s⁻¹). 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 mJy beam⁻¹).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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 M<sub>Jup</sub> after 4 × 10⁴ 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 M<sub>Jup</sub> after 4 × 10⁴ 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.</figcaption></figure>

respectively F<sub>B7</sub>, v<sub>B7</sub>, and F<sub>B4</sub>, v<sub>B4</sub>, the spectral index α<sub>B7,B4</sub> is
calculated as

$$
\alpha_{B7, B4} = \frac{\log_{10} [F_{B7}/F_{B4}]}{\log_{10} [v_{B7}/v_{B4}]}. \qquad (5)
$$

Assuming that opacity follows a power-law $\kappa_v \propto v^\beta$, the emission is optically thin and within the Rayleigh-Jeans regime, then the spectral index can be approximated as $\alpha_{B7,B4} \approx 2 + \beta$ (Testi et al. 2014). The value of $\beta$ depends on the maximum dust grain size, while not being influenced by the minimum grain size. A low $\beta$ (i.e., $\beta < 1$, meaning that $\alpha_{B7,B4} < 3$) could be indicative of the presence of grains with size $\ge 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<sub>Jup</sub> 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 frac-
tion decreases to ≈38% in Band 4. These results are consistent
with the findings in Tazzari et al. (2021), showing similar opti-
cally 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 pre-
dicted in our simulation (~20 K in the disk midplane, where most
<table><caption>**Table 3.** Fluxes in Band 7 and Band 4 and the resulting spatially integrated spectral index for the ALMA observation (accounting only for flux from the external dust ring) and our model from the simulation with a final planet mass of 2.0 $M_{\text{Jup}}$ after $4 \times 10^4$ yr.</caption><thead><tr><th></th><th>Observation<br>(external ring)</th><th>Simulation<br>$M_P = 2.0 M_{\text{Jup}}$</th></tr></thead><tbody><tr><td>$F_{B7}$ (mJy)</td><td>29.0 ± 2.9</td><td>28.6 ± 2.9</td></tr><tr><td>$F_{B4}$ (mJy)</td><td>4.7 ± 0.5</td><td>4.0 ± 0.4</td></tr><tr><td>Spectral index $\alpha_{B7, B4}$</td><td>2.08 ± 0.16</td><td>2.26 ± 0.16</td></tr></tbody></table>

**Notes.** The areas considered when calculating the spectral index are the same used for the dust mass rescaling process (Appendix C). The uncertainties for the observations are obtained as the root sum square of the observed noise ($\sigma$) and a 10% flux error due to calibration. The uncertainties for our simulation only include the 10% flux error.

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 $\alpha_{B7, B4} \approx 1.7 + \beta$.

Pinilla et al. (2021) find that the spatially integrated spectral index of CIDA 1 is $2.0 \pm 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_{\text{Jup}}$ after $4 \times 10^4$ yr. We note that the spatially integrated spectral index 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_{\text{Jup}}$ after $4 \times 10^4$ 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 ~14 au and ~20 au 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 ~2-2.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.

### 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 $^{12}$CO ($J = 3-2$) and $^{13}$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_{\text{Jup}}$ after $4 \times 10^4$ 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 $^{12}$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 $^{13}$CO emission. Figure 11 compares the simulated $^{13}$CO emission at three different angular resolutions: 0.03", 0.02", and 0.01". We show the channels at $v = 5.0 \text{ km s}^{-1}$ and $v = 7.5 \text{ km s}^{-1}$. These velocities, equidistant from the systemic velocity, are just outside the window where the cloud totally absorbs $^{13}$CO emission. We only show the results with a channel width of $0.5 \text{ km s}^{-1}$ since a higher velocity resolution does not lead to any major change. The planet is located at the azimuthal angle $\phi = 145^{\circ}$ 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 \text{ km s}^{-1}$ channel, becoming more evident at 0.01". Instead, no kink is detected in the $7.5 \text{ km s}^{-1}$ channel. Therefore, this analysis proves that a planet-induced kink would be visible only in $^{13}$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 $^{13}$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 $^{13}$CO and not in $^{12}$CO in our simulations. First, the emission of $^{12}$CO comes from regions at high altitudes, whose dynamics could be less affected by the planet orbiting on the midplane (Rabago & Zhu 2021), whereas $^{13}$CO is emitted closer to the midplane. Second, the spatial resolution of our models naturally becomes coarser at higher altitudes, where $^{12}$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 $^{12}$CO (Pinte et al. 2018, 2019, 2020). The explanation is that $^{12}$CO has a bright line emission, allowing for very high signal-to-noise ratios even at high spectral and spatial resolution. Kinks in $^{13}$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_{\text{Jup}}$ and $1.4 M_{\text{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,0}$ 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_{\text{Jup}}$ and $1.4 M_{\text{Jup}}$,
¹³CO Simulation

<figure><img src="image_10.png" /><figcaption>Fig. 11. Synthetic images of the ¹³CO (J = 3–2) emission from our simulation with M<sub>P</sub> = 2.0 M<sub>Jup</sub> after an evolution time of 4 × 10<sup>4</sup> yr. We compute the images at two different channels: v = 5.0 km s<sup>-1</sup> (top panels) and v = 7.5 km s<sup>-1</sup> (bottom panels), with a velocity resolution of 0.5 km s<sup>-1</sup>. 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 corner 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).</figcaption></figure>

that is, M<sub>P,0</sub> = 0.5 and 1.0 M<sub>Jup</sub>, respectively, we explore different values of d<sub>P,0</sub>. The best cases are presented in Fig. 12, after an evolution time of 4 × 10<sup>4</sup> yr.

The planet with M<sub>P</sub> = 0.9 M<sub>Jup</sub> and d<sub>P</sub> = 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<sub>P</sub> = 1.4 M<sub>Jup</sub> and d<sub>P</sub> = 10.3 au in both Band 7 and Band 4. Therefore, we conclude that a minimum mass of ~1.4 M<sub>Jup</sub> 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 α<sub>SS</sub> = 5 × 10<sup>-3</sup>. 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 α<sub>SS</sub> = 10<sup>-3</sup> to 10<sup>-4</sup>), 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.

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<sub>Jup</sub>, 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<sub>P,0</sub> = 4.0 M<sub>Jup</sub>. Choosing an initial location d<sub>P,0</sub> = 9.0 au, after 4 × 10<sup>4</sup> yr the planet reached d<sub>P</sub> = 8.5 au with M<sub>P</sub> = 4.5 M<sub>Jup</sub>, 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 (α<sub>SS</sub> ≲ 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}}, \qquad (6) $$

where Δ is defined as the distance between the minimum intensity in the gap and the peak intensity in the ring. Hence, we
<figure><img src="image_11.png" /><figcaption>**Fig. 12.** As in Fig. 5, but for the cases of simulations with final planet masses 0.9 and 1.4 *M*<sub>Jup</sub> with different initial distances from the star. For the simulation in the *first row*, *M*<sub>P,0</sub> = 0.5 *M*<sub>Jup</sub> and *d*<sub>P,0</sub> = 12.0 au, and, after an evolution time of 4 × 10<sup>4</sup> yr, *M*<sub>P</sub> = 0.9 *M*<sub>Jup</sub> and *d*<sub>P</sub> = 11.0 au. *Second row*: *M*<sub>P,0</sub> = 1.0 *M*<sub>Jup</sub> and *d*<sub>P,0</sub> = 11.0 au, and, after the same evolution time, *M*<sub>P</sub> = 1.4 *M*<sub>Jup</sub> and *d*<sub>P</sub> = 10.3 au.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 13.** As in Fig. 5, but for the case of a simulation with higher planet mass. The planet started with *M*<sub>P,0</sub> = 4.0 *M*<sub>Jup</sub> and *d*<sub>P,0</sub> = 9.0 au, and, after an evolution time of 4 × 10<sup>4</sup> yr, it reached *M*<sub>P</sub> = 4.5 *M*<sub>Jup</sub> and *d*<sub>P</sub> = 8.5 au.</figcaption></figure>

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*<sub>Jup</sub> 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*<sub>Jup</sub>. 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 α<sub>SS</sub> = 10<sup>-3</sup>, 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
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

Crida et al. (2006) criterion and estimated that the minimum planet mass to open a gap in the gas is ~0.4 $M_{\text{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 $M_{\text{Jup}}$ planet located at 8 au and a 0.5 $M_{\text{Jup}}$ 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 $M_{\text{Jup}}$ 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 $M_{\text{Jup}}$ 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_{\odot}$, whereas Pinilla et al. (2021) assumed a 0.1 $M_{\odot}$ central star. This difference in the stellar mass, along with a disk viscosity

$\alpha_{\text{SS}} = 10^{-3}$, lower than the initial value of $5 \times 10^{-3}$ adopted in our simulations, might explain why Pinilla et al. (2021) predict a planet mass of ~0.5 $M_{\text{Jup}}$ while our models require a minimum planet mass of ~1.4 $M_{\text{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.

### 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_{\odot}$ 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 for-
mation 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_{\odot}$. First, they calculated the initial
masses of protoplanets by extrapolating previous numerical sim-
ulations 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 peb-
ble accretion. Their results show that, around a 0.01 $M_{\odot}$ brown
dwarf, planets can grow up to 0.1–0.2 $M_{\oplus}$, while, around 0.1 $M_{\odot}$
stars, planets can reach a maximum mass of 2–3 $M_{\oplus}$. Find-
ings 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 synthe-
sis approach based on planetesimal accretion to explore planet
formation in the stellar mass range between 0.05 and 0.25 $M_{\odot}$.
They let the synthetic population of planetary systems evolve
for $10^8$ yr. The authors find that to form planets with masses
higher than 0.1 $M_{\oplus}$ they need stars of at least 0.07 $M_{\odot}$, implying
that planet formation around brown dwarf may not be a usual
outcome. Then, stars with masses higher than 0.15 $M_{\odot}$ are nec-
essary 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 incom-
plete, 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^7$ 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 sim-
ulations 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. Fragmenta-
tion 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_{\text{disk}}/M_{\star} \ge 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 consid-
ered stars with masses between 0.2 and 0.4 $M_{\odot}$. The authors
find that, via gravitational instabilities, protoplanets form very
fast, within a few thousand years, provided that the disk is mas-
sive enough ($M_{\text{disk}}/M_{\star} \ge 0.3$). The formed planets have masses
between 2 and 6 $M_{\text{Jup}}$ and they are initially distant from the cen-
tral 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. Assum-
ing 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_{\text{Jup}}$
and starting at a radial distance of between 9 and 11 au from
the central star. Then, we computed the dust and gas emis-
sion 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_{\text{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 max-
  imum 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_{\text{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_{\odot}$, 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_{\text{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 iso-
  topologue 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_{\odot}$ 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_{\text{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_{\text{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 C52113000140001), 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 Skłodowska-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 Skłodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project).

## References

Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766

Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 751

Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, MNRAS, 504, 5444

Briceno, C., Calvet, N., Gomez, M., et al. 1993, PASP, 105, 686

Bruderer, S. 2013, A&A, 559, A46

Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91

Chauvin, G., Lagrange, A. M., Dumas, C., et al. 2005, A&A, 438, L25

Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Clarke, C. J., Tazzari, M., Juhasz, A., et al. 2018, ApJ, 866, L6

Comeron, F., Rieke, G. H., Claes, P., Torra, J., & Laureijs, R. J. 1998, A&A, 335, 522

Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587

Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

Damasso, M., Del Sordo, F., Anglada-Escudé, G., et al. 2020, Sci. Adv., 6, eaax7467

Davis, C. J., Chrysostomou, A., Hatchell, J., et al. 2010, MNRAS, 405, 759

Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73

Dodson-Robinson, S. E., & Salyk, C. 2011, ApJ, 738, 131

Draine, B. T. 2003, ApJ, 598, 1026

Draine, B. T. 2006, ApJ, 636, 1114

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104

Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115

Fung, J., & Chiang, E. 2016, ApJ, 832, 105

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456

Haworth, T. J., Cadman, J., Meru, F., et al. 2020, MNRAS, 494, 4130

Hendler, N. P., Mulders, G. D., Pascucci, I., et al. 2017, ApJ, 841, 116

Hutchison, M., Price, D. J., & Laibe, G. 2018, MNRAS, 476, 2186

Izquierdo, A. F., Galván-Madrid, R., Maud, L. T., et al. 2018, MNRAS, 478, 2505

Izquierdo, A. F., Smith, R. J., Glover, S. C. O., & van 2012a, MNRAS, 500, 5268

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021b, A&A, 650, A179

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Kurtovic, N. T., Pinilla, P., Long, F., et al. 2021, A&A, 645, A139

Laibe, G., & Price, D. J. 2014, MNRAS, 440, 2136

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Liebert, J., & Probst, R. G. 1987, ARA&A, 25, 473

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75

Liu, B., Lambrechts, M., Johansen, A., Pascucci, I., & Henning, T. 2020, A&A, 638, A88

Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212

Lodato, G., Delgado-Donate, E., & Clarke, C. J. 2005, MNRAS, 364, L91

Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, A3

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, 376, 127

Mercer, A., & Stamatellos, D. 2020, A&A, 633, A116

Miguel, Y., Cridland, A., Ormel, C. W., Fortney, J. J., & Ida, S. 2020, MNRAS, 491, 1998

Mohanty, S., Jayawardhana, R., & Basri, G. 2005, ApJ, 626, 498

Morales, J. C., Mustill, A. J., Ribas, I., et al. 2019, Science, 365, 1441

Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130

Mulders, G. D., Pascucci, I., Ciesla, F. J., & Fernandes, R. B. 2021, ApJ, 920, 66

Natta, A., & Testi, L. 2001, A&A, 376, L22

Natta, A., Testi, L., Comerón, F., et al. 2002, A&A, 393, 597

Okabe, A., Boots, B., Sugihara, K., & Chiu, S. N. 2000, Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, 2nd edn., Series in Probability and Statistics (John Wiley and Sons, Inc.)

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Payne, M. J., & Lodato, G. 2007, MNRAS, 381, 1597

Pinilla, P., Benisty, M., & Birnstiel, T. 2012a, A&A, 545, A81

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012b, A&A, 538, A114

Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95

Pinilla, P., Quiroga-Nuñez, L. H., Benisty, M., et al. 2017, ApJ, 846, 70

Pinilla, P., Natta, A., Manara, C. F., et al. 2018, A&A, 615, A95

Pinilla, P., Kurtovic, N. T., Benisty, M., et al. 2021, A&A, 649, A122

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62

Price, D. J. 2007, PASA, 24, 159

Price, D. J., Cuello, N., Pinte, C., et al. 2018a, MNRAS, 477, 1270

Price, D. J., Wurster, J., Tricco, T. S., et al. 2018b, PASA, 35, e031

Rabago, I., & Zhu, Z. 2021, MNRAS, 502, 5325

Ragusa, E., Dipierro, G., Lodato, G., Laibe, G., & Price, D. J. 2017, MNRAS, 464, 1449

Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87

Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

Ricci, L., Testi, L., Natta, A., et al. 2014, ApJ, 791, 20

Rilinger, A. M., Espaillat, C. C., & Macías, E. 2019, ApJ, 878, 103

Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790

Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114

Schaefer, G. H., Dutrey, A., Guilloteau, S., Simon, M., & White, R. J. 2009, ApJ, 701, 698

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tazzari, M., Clarke, C. J., Testi, L., et al. 2021, MNRAS, 506, 2804

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339

Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111

Toci, C., Lodato, G., Christiaens, V., et al. 2020a, MNRAS, 499, 2015

Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020b, ApJ, 888, L4

Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19

Ubeira Gabellini, M. G., Miotello, A., Facchini, S., et al. 2019, MNRAS, 486, 4638

Veronesi, B., Ragusa, E., Lodato, G., et al. 2020, MNRAS, 495, 1913

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2021, A&A, 653, A46

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Warren, S. G., & Brandt, R. E. 2008, J. Geophys. Res. (Atmos.), 113, D14220

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Williams, J. P., Cieza, L., Hales, A., et al. 2019, ApJ, 875, L9

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321
## Appendix A: Sound speed power law exponent from the SED

<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

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 (Beckwith et al. 1990). In the frequency range $k_B T_{out} \ll h\nu \ll 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{d\log_1(\lambda F_\lambda)}{d\log_1(\lambda)} = -4 + \frac{1}{q}. \qquad (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_{10}(\lambda F_\lambda)$ with respect to $\log_{10}(\lambda)$ in the wavelength range 2.5 – 100 µm (Fig. A.1), obtaining a slope of $\approx -0.62$. Using the (A.1), we derive that $q \approx 0.3$.

## 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^4$ 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 \times 10^4$ 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

<figure><img src="image_15.png" /><figcaption>**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 M$_{Jup}$.</figcaption></figure>

can calculate the planet orbital eccentricity:

$$ e = \frac{r_a - r_p}{r_a + r_p}, \qquad (B.1) $$

obtaining a value $e \approx 0.01$ after $8 \times 10^4$ yrs.

## 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 $\ll 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.0 mJy. In our synthetic images, without an inner disk, we consider the flux contained
<figure><img src="image_16.png" /><figcaption>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 M<sub>Jup</sub> after $4 \times 10^4$ yrs.</figcaption></figure>

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.

## 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 \times 10^4$ yrs. Compared to the results after $4 \times 10^4$ 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<sub>Jup</sub> 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<sub>Jup</sub> and $M_P = 2.6$ M<sub>Jup</sub> 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.

## 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 $^{12}$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 $^{12}$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 $^{13}$CO $J = 3-2$ emission. CIDA 1 is located in the subregion L1495 W. Since in this area the $^{13}$CO data are missing, we consider the nearby subregion L1595 E for our estimates.

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

$$ \frac{\int T_A^*(^{12}\text{CO}) dv}{\int T_A^*(^{13}\text{CO}) dv} = \frac{1 - e^{-\tau_{12}}}{1 - e^{-\tau_{13}}} \quad (E.1) $$

where $\tau_{13} = \tau_{12}/70$ is the $^{13}$CO optical depth. We integrate the antenna temperatures along velocity intervals with the same centers and width of the observed $^{12}$CO channel maps (top panels in Fig. 6). The resulting $^{12}$CO optical depth is presented in Fig. E.3. The morphology of the optical depth spectrum follows approximately the cloud absorption behavior on the $^{12}$CO observation, being strongest at ~6.0–7.0 km s<sup>-1</sup> while fading at ~4.0 km s<sup>-1</sup> and ~8.5 km s<sup>-1</sup>, 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 $^{12}$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.
<figure><img src="image_17.png" /><figcaption>Band 4 (2.1 mm)</figcaption></figure>

<figure><img src="image_18.png" /></figure>

**Fig. D.1.** As in Fig. 5, but for the case of simulations evolved for 8 × 10⁴ yrs.
<figure><img src="image_19.png" /><figcaption>**Fig. E.1.** Map of the Taurus molecular cloud obtained by the James Clerk Maxwell Telescope using the ¹²CO (J = 3 – 2) emission (integrated between 2.0 and 11.0 km s⁻¹). Areas also mapped with ¹³CO (J = 3 – 2) emission (integrated between 4.9 and 9.1 km s⁻¹) 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).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. E.3.** ¹²CO (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 ¹²CO channel maps.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. E.2.** ¹²CO (J = 3 – 2) and ¹³CO (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).</figcaption></figure>