11institutetext: Dipartimento di Fisica, Università degli Studi di Milano, via
Celoria 16, 20133 Milano, Italy
11email<EMAIL_ADDRESS>22institutetext: European Southern
Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany
33institutetext: Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300
RA Leiden, The Netherlands 44institutetext: INAF – Osservatorio Astrofisico di
Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy 55institutetext: School of
Cosmic Physics, Dublin Institute for Advanced Studies, 31 Fitzwilliams Place,
Dublin 2, Ireland 66institutetext: Max-Planck-Institut für Astronomie,
Königstuhl 17, 69117, Heidelberg, Germany 77institutetext: 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 88institutetext: Univ. Grenoble Alpes, CNRS, IPAG, 38000
Grenoble, France 99institutetext: Institute of Astronomy, University of
Cambridge, Madingley Road, CB3 0HA Cambridge, UK 1010institutetext: INAF –
Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate (LC), Italy
1111institutetext: Freie Universität Berlin, Institute of Geological Sciences,
Malteserstr. 74-100, 12249 Berlin, Germany 1212institutetext: Department of
Physics and Astronomy, California State University Northridge, 18111 Nordhoff
Street, Northridge, CA 91130, USA

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

P. Curone 11 A. F. Izquierdo 22 3 3 L. Testi 22 4 4 G. Lodato 11 S. Facchini
11 2 2 A. Natta 55 P. Pinilla 66 N. T. Kurtovic 66
C. Toci 11 M. Benisty 7788 M. Tazzari 99 F. Borsa 1010 M. Lombardi 11 C. F.
Manara 22
E. Sanchis 111122 L. Ricci 1212

(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\,\mathrm{M}_{\odot}\lesssim
M_{\star}\lesssim 0.3\,\mathrm{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\text{\,}\mathrm{mm}$) and Band 4
($2.1\text{\,}\mathrm{mm}$), along with its 12CO ($J=3{-}2$) and 13CO
($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\,\mathrm{M}_{\mathrm{Jup}}$ and initially located at a distance of
between 9 and $11\text{\,}\mathrm{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
${{\sim}}1.4\,\mathrm{M}_{\mathrm{Jup}}$ orbiting at a distance of
${{\sim}}9{-}$10\text{\,}\mathrm{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
${{\sim}}10^{-2}$. We are able to reproduce the low spectral index
(${{\sim}}2$) observed where the dust ring is detected, with a
${{\sim}}40{-}50\%$ fraction of optically thick emission. Assuming a 12CO
abundance of $5\times 10^{-5}$ and a 13CO abundance 70 times lower, our
synthetic images reproduce the morphology of the 12CO ($J=3{-}2$) and 13CO
($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\,\mathrm{M}_{\odot}$ and a systemic velocity
$\varv_{\mathrm{sys}}=$6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{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
${{\sim}}4{-}$8\text{\,}\mathrm{M}_{\mathrm{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 $\lesssim 0.3\,\mathrm{M}_{\odot}$ but still above the hydrogen-
burning limit of ${{\sim}}0.08\,\mathrm{M}_{\odot}$, below which the brown
dwarf regime begins (Liebert & Probst 1987). Particularly fascinating are the
cases of TRAPPIST-1, a ${\sim}0.085\,\mathrm{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
${\sim}0.12\,\mathrm{M}_{\odot}$ and hosts a terrestrial planet, a super-Earth
candidate, and a newly discovered sub-Earth candidate (Anglada-Escudé et al.
2016; Damasso et al. 2020, Faria et al. 2022). In general, data from the
Kepler spacecraft show that the occurrence rate of small planets (with radii
of $1.0{-}2.8\,\mathrm{R}_{\oplus}$) 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 ${\sim}5\,\mathrm{M}_{\mathrm{Jup}}$ planet around a
${\sim}25\,\mathrm{M}_{\mathrm{Jup}}$ young brown dwarf, and Morales et al.
(2019) discovered a planet with a minimum mass of
$0.46\,\mathrm{M}_{\mathrm{Jup}}$ orbiting a ${\sim}0.12\,\mathrm{M}_{\odot}$
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
(${\sim}0.1^{\prime\prime}$ at $0.87\text{\,}\mathrm{m}\mathrm{m}$, 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\arcsec$ (corresponding to ${\sim}$55\text{\,}\mathrm{au}$$ at
$137.5\text{\,}\mathrm{pc}$, the estimated distance to CIDA 1111Throughout 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 et al. 2021) and derived a distance of
$133.4^{+0.7}_{-0.9}\,\mathrm{pc}$. Such a distance differs by ${\approx}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 ${\sim}6\%$ in the estimated value of the disk mass.), Ricci et al. (2014)
showed a dust disk and detected the CO ($J=3{-}2$) line emission.

An inner cavity with a radius of ${\sim}$ 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\arcsec\times 0.12\arcsec$
(${\sim}29\times$17\text{\,}\mathrm{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^{-4}$ and a
stellar mass $M_{\star}=0.1\,\mathrm{M}_{\odot}$. They found that a planet
with a minimum mass of ${\sim}0.3\,\mathrm{M}_{\mathrm{Jup}}$ orbiting at a
distance of $15\text{\,}\mathrm{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 ${\sim}0.050^{\prime\prime}\times
0.034^{\prime\prime}$ (${\sim}7\times$5\text{\,}\mathrm{au}$$). These images
clearly reveal a dust ring, whose emission peaks at ${\sim}$ 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 12CO
($J=3{-}2$) and 13CO ($J=3{-}2$). The authors used 1D dust evolution
simulations and excluded that a planet more massive than
$M_{\star}=0.5\,\mathrm{M}_{\mathrm{Jup}}$ could cause the observed dust gap.

Figure 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\arcsec\times 0.034\arcsec$, indicated by the white ellipse in the
bottom-left corner.

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

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

##  2 Methods

###  2.1 Gas and dust hydrodynamical simulation

We performed a set of five hydrodynamical simulations with the smoothed
particle hydrodynamics (SPH) code Phantom (Price et al. 2018b). This code has
been extensively used for the modeling of protoplanetary disks, with the
intent of reproducing the observed structures in both the dust and gas
emission (e.g., Dipierro et al. 2015, Ragusa et al. 2017, Pinte et al. 2018,
Price et al. 2018a, Cuello et al. 2019, Ubeira Gabellini et al. 2019, Toci et
al. 2020a, Veronesi et al. 2020, Toci et al. 2020b). 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
$\mathrm{St}\propto\rho_{0}a_{\mathrm{d}}/\Sigma_{\mathrm{g}}<1$, where
$\rho_{0}$ is the dust grain intrinsic density, $a_{\mathrm{d}}$ the grain
size, and $\Sigma_{\mathrm{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_{\mathrm{SPH}}=10^{6}$ SPH particles. The disk has a gas
mass of $M_{\mathrm{gas}}=1.5\times 10^{-3}\,\mathrm{M}_{\odot}$ and initially
extends from $R_{\mathrm{in}}=$1\text{\,}\,\mathrm{au}$$ to
$R_{\mathrm{out}}=$100\text{\,}\,\mathrm{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_{\mathrm{g}}(R)=\Sigma_{\mathrm{c}}\left(\frac{R}{R_{\mathrm{in}}}\right)^{-p}\exp\left[-\left(\frac{R}{R_{\mathrm{c}}}\right)^{2-p}\right]\,\left(1-\sqrt{\frac{R_{\mathrm{in}}}{R}}\,\right)\,,$
(1)

where $\Sigma_{\mathrm{c}}$ is a normalization constant depending on the total
disk mass, $R_{\mathrm{c}}=$80\text{\,}\,\mathrm{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_{\mathrm{g}}c_{\mathrm{s}}^{2}$, with the following sound speed radial
profile:

$c_{\mathrm{s}}(R)=c_{\mathrm{s,in}}\biggl{(}\frac{R}{R_{\mathrm{in}}}\biggr{)}^{-q},$
(2)

where $\rho_{\mathrm{g}}$ is the gas volume density, $c_{\mathrm{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_{\mathrm{s}}(R)}{\varv_{\mathrm{K}}(R)}=\frac{H(R_{\mathrm{in}})}{R_{\mathrm{in}}}\,\left(\frac{R}{R_{\mathrm{in}}}\right)^{1/2-q}\,,$
(3)

where $\varv_{\mathrm{K}}$ is the Keplerian velocity and the aspect ratio at
the inner radius is $H(R_{\mathrm{in}})/R_{\mathrm{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_{\mathrm{AV}}$
introduced by Lodato & Price (2010). We set the Shakura & Sunyaev (1973)
viscosity $\alpha_{\mathrm{SS}}=5\times 10^{-3}$, corresponding to
$\alpha_{\mathrm{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).

Table 1: Initial parameters used for the simulations with Phantom.

$\begin{array}[]{p{0.3\linewidth}l}\hline\cr\hline\cr\vskip 3.0pt plus 1.0pt
minus 1.0pt\cr Parameters&\text{Value}\\\ \vskip 3.0pt plus 1.0pt minus
1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr$N_{\text{SPH}}$&10^{6}\\\ $M_{\star}[\mathrm{M}_{\odot}]$&$0.2$\\\
$R_{\star\,\text{acc}}$ [au]&1\\\ $R_{\text{in}}$ [au]&1\\\ $R_{\text{out}}$
[au]&100\\\ $R_{\text{c}}$ [au]&80\\\ $p$&1.5\\\ $q$&0.3\\\
$H(R_{\mathrm{in}})/R_{\mathrm{in}}$&0.08\\\ $\alpha_{\text{SS}}$&5\times
10^{-3}\\\ $\alpha_{\text{AV}}$&0.36\\\ $\langle h\rangle/H$&0.14\\\
$M_{\text{gas}}$ [$\mathrm{M}_{\odot}$]&1.5\times 10^{-3}\\\
$M_{\text{dust}}/M_{\text{gas}}$&10^{-2}\\\ $a_{\text{d}}$
[cm]&[a_{\text{min}}=1.6\times 10^{-5},\,a_{\text{max}}=0.1,N=11]\\\
$\xi$&3.5\\\ $\rho_{\text{0}}\,[\text{g cm}^{-3}]$&2\\\ \vskip 3.0pt plus
1.0pt minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus
1.0pt\cr$R_{\text{P}\,\text{acc}}$ [au]&1/4\text{ of Hill radius}\\\
$d_{\text{P},0}$ [au]&10.0\\\ $M_{\text{P},0}$
[$\mathrm{M}_{\mathrm{Jup}}$]&0.1,0.5,1.0,1.5,2.0\\\ \vskip 3.0pt plus 1.0pt
minus 1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\end{array}$ Notes.
Each entry is commented on in the main text.

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\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $1\text{\,}\mathrm{mm}$. The
grain size distribution follows the power law
$\mathrm{d}N(a_{\mathrm{d}})/\mathrm{d}a_{\mathrm{d}}\propto
a_{\mathrm{d}}^{-\xi}$, truncated at the minimum and maximum grain size, where
$a_{\mathrm{d}}$ indicates the grain size and $N(a_{\mathrm{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\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$$.

Figure 2: Face-on gas and dust surface density, in units of
$\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-2}$, from the Phantom simulation
with an initial planet mass of $1.5\,\mathrm{M}_{\mathrm{Jup}}$ and a final
planet mass of $2.0\,\mathrm{M}_{\mathrm{Jup}}$ after
${\approx}$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$$ (${\approx}560$
planet orbits). Gas density is displayed in the top-left corner, and the
remaining panels present the surface density for the 11 simulated dust grain
sizes. The bigger central green dot represents the star, and the smaller green
dot indicates the planet position. Images were produced using the tool Splash
(Price 2007).

####  2.1.2 Properties of the central star and the embedded planet

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

$R_{\mathrm{H}}=\left(\frac{1}{3}\frac{M_{\mathrm{P}}}{M_{\star}}\right)^{1/3}d_{\mathrm{P}}\,,$
(4)

where $d_{\mathrm{P}}$ is the distance of the planet from the central star.

We performed a set of SPH simulations, each with an embedded planet starting
at a distance $d_{\mathrm{P},0}=$10.0\text{\,}\mathrm{au}$$ from the star. We
varied the initial planet mass $M_{\mathrm{P},0}$: the chosen values are
$0.1$, $0.5$, $1.0$, $1.5$, and $2.0$ $\mathrm{M}_{\mathrm{Jup}}$. We let the
simulations evolve until they reach a stable configuration. We analyzed the
simulation results at $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$,
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\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, 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_{\mathrm{P},0}=1.5\,\mathrm{M}_{\mathrm{Jup}}$, evaluated after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. 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 ${\sim}0.1$ to $1\text{\,}\mathrm{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_{\mathrm{P},0}=0.5\,\mathrm{M}_{\mathrm{Jup}}$
and $d_{\mathrm{P},0}=$12.0\text{\,}\mathrm{au}$$, the second with
$M_{\mathrm{P},0}=1.0\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P},0}=$11.0\text{\,}\mathrm{au}$$. The maximum mass of the planet
is assessed in Sect. 4.5, considering a simulation with
$M_{\mathrm{P},0}=4.0\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P},0}=$9.0\text{\,}\mathrm{au}$$.

####  2.1.3 Limitations of the modeling

Figure 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_{\mathrm{d}}$ values associated with the two smallest simulated
grain sizes to obtain densities for grains down to $5\text{\,}\mathrm{nm}$
(blue line). Then, we use a cubic spline interpolation on the remaining
$\log\rho_{\mathrm{d}}$ values (orange line). In the  right panel, with the
linear fit, we obtain ten dust densities related to grain sizes between
$5\text{\,}\mathrm{nm}$ and $0.37\text{\,}\mathrm{\SIUnitSymbolMicro 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\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $1\text{\,}\mathrm{mm}$ (to
the right of the vertical dashed gray line).

High-resolution ALMA images of the dust continuum emission from CIDA 1 (Fig.
1) show a dust ring peaking at ${\sim}$20\text{\,}\mathrm{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_{\star\,\mathrm{acc}}=$1\text{\,}\mathrm{au}$$;
that is, we are not simulating the dynamics of the material within a radius of
$1\text{\,}\mathrm{au}$ from the central star. Adopting a smaller
$R_{\mathrm{acc}}$ rapidly causes a substantial increase in computational time
due to a higher dynamic range. SPH particles closer to the central star move
faster, so smaller time steps and a greater computational cost are needed to
simulate their motion.

Furthermore, an intrinsic consequence of the artificial viscosity
implementation in Phantom is that $\alpha_{\mathrm{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 sf  3 dmodels 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\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$1\text{\,}\mathrm{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.,
$\mathrm{d}N(a_{\mathrm{d}})/\mathrm{d}a_{\mathrm{d}}\propto
a_{\mathrm{d}}^{-3.5}$, which can be rearranged as $\rho_{\mathrm{d}}\propto
a_{\mathrm{d}}^{0.5}$, with $\rho_{\mathrm{d}}$ indicating the dust mass
density), we fitted the values of $\log_{10}\rho_{\mathrm{d}}$ associated with
the two smallest simulated grain sizes
($0.16\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$) using a straight line with slope
${\sim}0.5$. With this relation, we sampled 10 density values logarithmically
spaced between $5\text{\,}\mathrm{nm}$ and
$0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$. For larger dust grains, we
performed a cubic spline interpolation for $\log_{10}\rho_{\mathrm{d}}$ values
related to the remaining simulated grain sizes. Thanks to this interpolation,
we sampled 50 density values, logarithmically spaced between
$0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $1\text{\,}\mathrm{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\text{\,}\mathrm{nm}$ and
$0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and 50 related to grain sizes
between $0.37\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$1\text{\,}\mathrm{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{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$). We
randomly distributed $1/5\,N_{\mathrm{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_{\star}=0.15\pm 0.03\,\mathrm{L}_{\sun}$. 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_{\mathrm{eff}}=$3050\text{\,}\mathrm{K}$$ and a radius
$R_{\star}=1.3\,\mathrm{R}_{\sun}$, 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\text{\,}\mathrm{pc}$. From the analysis of Pinilla et al. (2021), we
used an inclination angle $i=37.5^{\circ}$ and a position angle
$\mathrm{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\text{\,}\mathrm{mm}$)
and Band 4 ($2.1\text{\,}\mathrm{mm}$). The gas line emission is calculated
for 12CO and 13CO 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\text{\,}\mathrm{km}\text{\,}{\mathrm{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\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. We included freeze-out and
photo-dissociation for CO following the parameterizations in Williams & Best
(2014). The condition $\mathrm{T}_{\mathrm{gas}}<$20\text{\,}\mathrm{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_{\mathrm{dissoc}}=$1.3\text{\times}{10}^{21}\text{\,}\mathrm{c}\mathrm{m}^{-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 12CO abundance relative to $\mathrm{H}_{2}$
$[^{12}\mathrm{CO}/\mathrm{H}_{2}]=5\times 10^{-5}$, and the isotopologue
abundance ratio $[^{12}\mathrm{CO}/^{13}\mathrm{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\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ for 12CO and
$1.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ for 13CO, we first
computed the simulated channel maps at a higher velocity resolution of
$0.045\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. Then, to obtain the
final spectrally convoluted synthetic channels, we averaged 11 contiguous
simulated 12CO channels, centered at the same velocities of the observed
channels. We repeated the same procedure for the computed 13CO 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,\varv)$ points of the ALMA observations, in Band 7 and Band 4,
respectively. After that, we calculated the residuals between the observations
and the simulations. Finally, we performed the imaging and acquired the
synthetic observations using the tclean algorithm from the software casa,
version 5.7 (McMullin et al. 2007). We applied the same parameters used to
obtain the dust continuum Band 7 and Band 4 images in Pinilla et al. (2021):
the Briggs weighting scheme and a robust parameter of 0.5.

At this point, we can compare the flux in our synthetic images to the one from
the observations. We consider the total dust mass of our models as a free
parameter that can be rescaled by a constant factor in the radiative transfer
simulations to match the flux emitted in Band 7 ($0.9\text{\,}\mathrm{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 12CO and 13CO 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\arcsec$ to recover the same beam size of the observations.

##  3 Results

Table 2: Resulting parameters for our simulations with different initial
planet masses.

$M_{\text{P},0}$ | $t_{\mathrm{Sim}}$ | $M_{\text{P}}$ | $d_{\text{P}}$ | $M_{\text{dust}}$ | $M_{\text{dust}}/M_{\text{gas}}$
---|---|---|---|---|---
[$\mathrm{M}_{\mathrm{Jup}}$] | [yrs] | [$\mathrm{M}_{\mathrm{Jup}}$] | [au] | [$\mathrm{M}_{\odot}$] |
$0.1$ | $4\times 10^{4}$ | 0.2 | 9.4 | $6.0\times 10^{-6}$ | $6.3\times 10^{-3}$
$8\times 10^{4}$ | 0.2 | 9.2 | $4.3\times 10^{-6}$ | $5.6\times 10^{-3}$
$0.5$ | $4\times 10^{4}$ | 0.9 | 9.3 | $4.8\times 10^{-6}$ | $5.8\times 10^{-3}$
$8\times 10^{4}$ | 0.9 | 9.4 | $3.3\times 10^{-6}$ | $5.1\times 10^{-3}$
$1.0$ | $4\times 10^{4}$ | 1.4 | 9.4 | $7.0\times 10^{-6}$ | $8.7\times 10^{-3}$
$8\times 10^{4}$ | 1.5 | 9.4 | $4.0\times 10^{-6}$ | $6.4\times 10^{-3}$
$1.5$ | $4\times 10^{4}$ | 2.0 | 9.4 | $7.0\times 10^{-6}$ | $8.9\times 10^{-3}$
$8\times 10^{4}$ | 2.0 | 9.5 | $5.4\times 10^{-6}$ | $8.8\times 10^{-3}$
$2.0$ | $4\times 10^{4}$ | 2.5 | 9.4 | $7.3\times 10^{-6}$ | $9.4\times 10^{-3}$
$8\times 10^{4}$ | 2.6 | 9.6 | $6.4\times 10^{-6}$ | $1.1\times 10^{-2}$

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

Figure 4: Band 7 ($0.9\text{\,}\mathrm{mm}$, top panels) and Band 4
($2.1\text{\,}\mathrm{m}\mathrm{m}\mathrm{,}$ 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\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, along with their residuals. The
synthesized beam ($0.050\arcsec\times 0.034\arcsec$ for Band 7,
$0.048\arcsec\times 0.032\arcsec$ 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.
Figure 5: Comparison between radial profiles of azimuthally averaged
intensities from the observations, simulations, and corresponding residuals,
after the simulations have evolved for
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. 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\,\mathrm{M}_{\mathrm{Jup}}$ (top row) to
$2.5\,\mathrm{M}_{\mathrm{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.

Figure 6: Comparison between the ALMA observation (top panels) and the
synthetic images (bottom panels) of CIDA 1 12CO ($J=3{-}2$) channel maps.
Color scales are the same in the top and bottom panels. Synthetic images are
obtained from the simulation with a final planet mass of
$2.0\;\mathrm{M}_{\mathrm{Jup}}$ after an evolution time of
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. We fixed the systemic velocity
at $6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The contour level
traces the 15$\sigma$ emission from the Band 7 continuum image. The velocity
resolution is $0.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, and the
central velocities of each channel are indicated on the top of each panel. The
synthesized beam ($0.101\arcsec\times 0.075\arcsec$, FWHM) is shown in the
lower-left corner of each panel. Figure 7: Contour levels (left panel in each
subplot) at $5,\,10,\,$15\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$$, and
intensity profile comparisons (right panel in each subplot) between the 12CO
observation and simulation with a $2.0\;\mathrm{M}_{\mathrm{Jup}}$ planet mass
after $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, in the six channels not
completely obscured by the cloud
($\varv=4.0,\,4.5,\,5.0,\,8.0,\,8.5,\,$9.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$).
Intensity profiles are taken along the $\mathrm{PA}=11^{\circ}$ 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\sigma$
noise level ($1.6\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$).

Figure 8: Comparison between the ALMA observation (top panels) and the
synthetic images (bottom panels) of CIDA 1 13CO ($J=3{-}2$) channel maps.
Color scales are the same in the top and bottom panels. Synthetic images are
obtained from the simulation with a final planet mass of
$2.0\;\mathrm{M}_{\mathrm{Jup}}$ after an evolution time of
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. We fixed the systemic velocity
at $6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The contour level
traces the 15$\sigma$ emission from the Band 7 continuum image. The velocity
resolution is $1\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, and the
central velocities of each channel are indicated on the top of each panel. The
synthesized beam ($0.104\arcsec\times 0.076\arcsec$, FWHM) is shown in the
lower-left corner of each panel. Figure 9: Contour levels (left panel in each
subplot) at $3,\,6,\,$9\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$$, and
intensity profile comparisons (right panel in each subplot) between the 13CO
observation and simulation with a $2.0\;\mathrm{M}_{\mathrm{Jup}}$ planet mass
after $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, in the four channels not
completely obscured by the cloud
($\varv=4.5,\,5.5,\,7.5,\,$8.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$).
Intensity profiles are taken along the $\mathrm{PA}=11^{\circ}$ 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\sigma$
noise level ($1.2\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$).

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_{\mathrm{Sim}}=$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$$ and
$t_{\mathrm{Sim}}=$8\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$$. We note
little difference in planet mass ($M_{\mathrm{P}}$) and distance from the
central star ($d_{\mathrm{P}}$) in each simulation at the two times, meaning
that $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ 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\;\mathrm{M}_{\mathrm{Jup}}$. Here, we note a fast
transient in the first ${10}^{4}\text{\,}\mathrm{yrs}$, when the planet is
carving the gap, followed by a more stable evolution. At this later stage, the
system has adapted from the arbitrary initial conditions and reached a quasi-
equilibrium state. If the planet is massive enough, such a situation
corresponds to the case in which the planet has already carved an internal
cavity both in the gas and in the dust, as in the case illustrated in Fig. 2.

###  3.1 Dust continuum emission

The first aim of our work is reproducing the dust thermal emission of the dust
ring in CIDA 1 observed by ALMA in Band 7 ($0.9\text{\,}\mathrm{mm}$) and Band
4 ($2.1\text{\,}\mathrm{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\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, 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\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ 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\,\mathrm{M}_{\mathrm{Jup}}$, does not produce any major effect on the
disk morphology. The $0.9\,\mathrm{M}_{\mathrm{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\,\mathrm{M}_{\mathrm{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\,\mathrm{M}_{\mathrm{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\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, 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}\,\mathrm{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\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ are presented in
Appendix D. The synthetic profiles for simulations with a minimum planet mass
of $2.0\,\mathrm{M}_{\mathrm{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\text{\times}{10}^{4}\text{\,}$ and
$8\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ 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.

###  3.2 Gas line emission

The dust continuum modeling shown in Fig. 5 proves that, while obtaining a
good match in the case of a $1.4\,\mathrm{M}_{\mathrm{Jup}}$ planet, the
minimum planet mass able to best replicate the observed dust emission is
$2.0\,\mathrm{M}_{\mathrm{Jup}}$. Therefore, we decided to use the model with
such a planet to compute the 12CO ($J=3{-}2$) and 13CO ($J=3{-}2$) channel
maps.

In Fig. 6 we compare the 12CO channel maps observed by ALMA to our synthetic
gas emission images. In our model, we set a systemic velocity of
$6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. Strong cloud absorption
is evident in the central channels of the observation, from ${\sim}5.5$ to
$7.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. We analyze the cloud
contamination effects in Appendix E, using previous observations of the gas
emission from the Taurus molecular cloud by Davis et al. (2010).

Our simulation reproduces appropriately the spatial extent of the channels
centered at the velocities
$\varv=4.0,\,4.5,\,5.0,\,8.0,\,8.5,\,$9.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$.
This is more evident in Fig. 7, which shows for these channels the comparison
between the observed and simulated 12CO emission in terms of intensity contour
levels at $5,\,10,\,$15\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$$ and
intensity profiles along the $\mathrm{PA}=11^{\circ}$. The channels
$\varv=4.5,\,5.0,\,$8.0\text{\,}\mathrm{km}\text{\,}{\mathrm{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
$\varv=4.0,\,8.5,\,$9.0\text{\,}\mathrm{km}\text{\,}{\mathrm{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 $\varv=$3.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and
$\varv=9.5{-}$10.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$. 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 13CO emission. Figure 8
presents the comparison between the 13CO channel maps observed by ALMA and the
synthetic ones computed from our simulation. Here, we have a lower velocity
resolution of $1.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. In the
observations, we note that the
$6.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ channel is completely
affected by cloud absorption, which also impacts the channel at
$\varv=$5.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$. In Fig. 9, we
compare the contour levels at
$3,\,6,\,$9\text{\,}\mathrm{mJy}\,\mathrm{beam}^{-1}$$ and the intensity
profiles at $\mathrm{PA}=11^{\circ}$ from the observed and simulated 13CO
emission. Even in this case, our models are able to reproduce the spatial
distribution of the observed emission while presenting higher intensity
values, except for the $8.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$
channel. In particular, we note a strong cloud absorption for the observed
intensity profile of the $5.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$
channel, but the corresponding contour plot shows a good agreement in the
extent of the emission between observation and simulation. As for the 12CO,
there is some observed inner emission for the high-velocity channels at $3.5$
and $9.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, 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\,\mathrm{M}_{\mathrm{Jup}}$ is able to globally reproduce
both the dust and gas emission observed in CIDA 1. In particular, this case
perfectly recovers the external dust ring as observed both in Band 7 and Band
4, along with replicating the spatial distribution of the gas emission as
traced by 12CO ($J=3{-}2$) and 13CO ($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 12CO and 13CO channel maps are strongly
affected by cloud absorption in the channels near the systemic velocity (Sect.
3.2 and Appendix E). Hence, we compare our synthetic images to the
observations only in those channels probing higher velocities in the disk
(Fig. 7 and 9). Models are consistent with observations, being able to
reproduce the spatial distribution of the emission. The stellar mass
$M_{\star}=0.2\,\mathrm{M}_{\odot}$ and systemic velocity
$\varv_{\mathrm{sys}}=$6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
assumed in our simulations are therefore reasonable estimates of these
properties in CIDA 1. This stellar mass is compatible with the known range
${\sim}0.1-0.2\,\mathrm{M}_{\odot}$ (Pinilla et al. 2021) and is in very good
agreement with the estimate of $0.19\,\mathrm{M}_{\odot}$ by Kurtovic et al.
(2021). We also simulated the case with $M_{\star}=0.1\,\mathrm{M}_{\odot}$,
but the morphology of the resulting gas channel maps differed significantly
from the observed ones; therefore, we decided not to include these models in
the simulation sample.

###  4.2 Spectral index and optical depth

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\;\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$.

| Observation | Simulation
---|---|---
| (external ring) | $M_{\mathrm{P}}=2.0\,\mathrm{M}_{\mathrm{Jup}}$
$F_{B7}$ [mJy] | $29.0\pm 2.9$ | $28.6\pm 2.9$
$F_{B4}$ [mJy] | $4.7\pm 0.5$ | $4.0\pm 0.4$
Spectral index $\alpha_{\mathrm{B7},\,\mathrm{B4}}$ | $2.08\pm 0.16$ | $2.26\pm 0.16$

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.

Figure 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\;\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, 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 $\tau=1$ level. The region
where the emission is optically thick (optical depth $\tau\geq 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\;\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. 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 ($\tau\geq 1$) in both Band 7 and Band 4.

Having dust emission data at two different wavelengths allows us to compare
the observed and simulated spectral index. Given the fluxes and frequencies
corresponding to Band 7 and Band 4, respectively $F_{\mathrm{B7}}$,
$\nu_{\mathrm{B7}}$, and $F_{\mathrm{B4}}$, $\nu_{\mathrm{B4}}$, the spectral
index $\alpha_{\mathrm{B7},\,\mathrm{B4}}$ is calculated as

$\alpha_{\mathrm{B7},\,\mathrm{B4}}=\frac{\log_{10}\left[F_{\mathrm{B7}}/F_{\mathrm{B4}}\right]}{{\log_{10}\left[\nu_{\mathrm{B7}}/\nu_{\mathrm{B4}}\right]}}\,.$
(5)

Assuming that opacity follows a power-law $\kappa_{\nu}\propto\nu^{\beta}$,
the emission is optically thin and within the Rayleigh-Jeans regime, then the
spectral index can be approximated as
$\alpha_{\mathrm{B7},\,\mathrm{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_{\mathrm{B7},\,\mathrm{B4}}<3$) could be indicative of the
presence of grains with size $\gtrsim$0.1\text{\,}\mathrm{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\;\mathrm{M}_{\mathrm{Jup}}$ 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 $\tau$ is above 1. We calculate that ${\approx}51\%$ of the
emission in Band 7 is optically thick, whereas this fraction decreases to
${\approx}38\%$ in Band 4. These results are consistent with the findings in
Tazzari et al. (2021), showing similar optically thick fractions in a sample
of disks located in the Lupus star-forming region. Such a mixture of optically
thin and thick emission in our model leads us to interpret with caution the
spectral index in CIDA 1. Moreover, the dust temperature predicted in our
simulation (${\sim}$20\text{\,}\mathrm{K}$$ in the disk midplane, where most
of the settled dust is located) lies at the limit of the Rayleigh-Jeans
domain. In this spectral region, and for a temperature of
$20\text{\,}\mathrm{K}$, dropping the Rayleigh-Jeans approximation but
retaining the optically thin one, we calculate that the spectral index is
$\alpha_{\mathrm{B7},\,\mathrm{B4}}\approx 1.7+\beta$.

Figure 11: Synthetic images of the 13CO ($J=3{-}2$) emission from our
simulation with $M_{\mathrm{P}}=2.0\,\mathrm{M}_{\mathrm{Jup}}$ after an
evolution time of $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. We compute
the images at two different channels:
$\varv=$5.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ (top panels) and
$\varv=$7.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ (bottom panels),
with a velocity resolution of
$0.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. We use three different
circular beam sizes: $0.03\arcsec$, $0.02\arcsec$, and $0.01\arcsec$ (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 $\phi=145^{\circ}$ (relative to the angle $\phi=0^{\circ}$, which is
located at the disk major axis indicated by the dashed white line in the top-
left panel).

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\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. 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\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. Our model best reproduces the
observed spectral index profile between ${\sim}$10\text{\,}\mathrm{au}$$ and
${\sim}$30\text{\,}\mathrm{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\text{\,}\mathrm{au}$ from the central
star. Then, we note that the spectral indices start diverging beyond
${\sim}$30\text{\,}\mathrm{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
${\sim}$14\text{\,}\mathrm{au}$$ and ${\sim}$20\text{\,}\mathrm{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
${\sim}$20\text{\,}\mathrm{au}$$ and ${\sim}$30\text{\,}\mathrm{au}$$, still
showing a spectral index around ${\sim}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 12CO ($J=3{-}2$) and 13CO
($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.

Figure 12: As in Fig. 5, but for the cases of simulations with final planet
masses $0.9$ and $1.4\,\mathrm{M}_{\mathrm{Jup}}$ with different initial
distances from the star. For the simulation in the first row,
$M_{\mathrm{P},0}=0.5\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P},0}=$12.0\text{\,}\mathrm{au}$$, and, after an evolution time of
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$,
$M_{\mathrm{P}}=0.9\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P}}=$11.0\text{\,}\mathrm{au}$$. Second row:
$M_{\mathrm{P},0}=1.0\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P},0}=$11.0\text{\,}\mathrm{au}$$, and, after the same evolution
time, $M_{\mathrm{P}}=1.4\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P}}=$10.3\text{\,}\mathrm{au}$$.

We consider our simulation with a final planet mass of
$2.0\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. 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 12CO emission, there is no
evidence of kinks in any channels. On the other hand, a local perturbation
nearby the planet is visible in the 13CO emission. Figure 11 compares the
simulated 13CO emission at three different angular resolutions: $0.03\arcsec$,
$0.02\arcsec$, and $0.01\arcsec$. We show the channels at
$\varv=$5.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and
$\varv=$7.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$. These velocities,
equidistant from the systemic velocity, are just outside the window where the
cloud totally absorbs 13CO emission. We only show the results with a channel
width of $0.5\text{\,}\mathrm{km}\text{\,}{\mathrm{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\arcsec$ in the
$5.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ channel, becoming more
evident at $0.01\arcsec$. Instead, no kink is detected in the
$7.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ channel. Therefore, this
analysis proves that a planet-induced kink would be visible only in 13CO
emission and at a very high resolution of at least $0.02\arcsec$, 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 13CO 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 13CO and not in 12CO in our simulations. First, the emission of
12CO comes from regions at high altitudes, whose dynamics could be less
affected by the planet orbiting on the midplane (Rabago & Zhu 2021), whereas
13CO is emitted closer to the midplane. Second, the spatial resolution of our
models naturally becomes coarser at higher altitudes, where 12CO 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 12CO (Pinte et al. 2018,
2019, 2020). The explanation is that 12CO has a bright line emission, allowing
for very high signal-to-noise ratios even at high spectral and spatial
resolution. Kinks in 13CO 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\;\mathrm{M}_{\mathrm{Jup}}$ and $1.4\;\mathrm{M}_{\mathrm{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_{\mathrm{P},0}$ from the central star of
$10.0\text{\,}\mathrm{au}$, now, we aim to better constrain the minimum planet
mass able to explain the observed dust ring by modifying the parameter
$d_{\mathrm{P},0}$. Using the same initial planet masses that eventually
resulted in $0.9\;\mathrm{M}_{\mathrm{Jup}}$ and
$1.4\;\mathrm{M}_{\mathrm{Jup}}$, that is, $M_{\mathrm{P},0}=0.5$ and
$1.0\;\mathrm{M}_{\mathrm{Jup}}$, respectively, we explore different values of
$d_{\mathrm{P},0}$. The best cases are presented in Fig. 12, after an
evolution time of $4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$.

The planet with $M_{\mathrm{P}}=0.9\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P}}=$11.0\text{\,}\mathrm{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 ${\sim}4$ and ${\sim}$10\text{\,}\mathrm{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_{\mathrm{P}}=1.4\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P}}=$10.3\text{\,}\mathrm{au}$$ in both Band 7 and Band 4.
Therefore, we conclude that a minimum mass of
${\sim}1.4\,\mathrm{M}_{\mathrm{Jup}}$ can explain the observed dust ring.

Such constraints on the planet mass have been obtained by varying the initial
planet mass and location but always fixing the initial value of the disk
viscosity at $\alpha_{\mathrm{SS}}=5\times 10^{-3}$. 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
$\alpha_{\mathrm{SS}}=10^{-3}$ to $10^{-4}$), 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

Figure 13: As in Fig. 5, but for the case of a simulation with higher planet
mass. The planet started with
$M_{\mathrm{P},0}=4.0\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P},0}=$9.0\text{\,}\mathrm{au}$$, and, after an evolution time of
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$, it reached
$M_{\mathrm{P}}=4.5\,\mathrm{M}_{\mathrm{Jup}}$ and
$d_{\mathrm{P}}=$8.5\text{\,}\mathrm{au}$$. Figure 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\text{\,}\mathrm{au}$ from the star..

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\,\mathrm{M}_{\mathrm{Jup}}$,
starting at a distance of $10.0\text{\,}\mathrm{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_{\mathrm{P},0}=4.0\,\mathrm{M}_{\mathrm{Jup}}$. Choosing an initial
location $d_{\mathrm{P},0}=$9.0\text{\,}\mathrm{au}$$, after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ the planet reached
$d_{\mathrm{P}}=$8.5\text{\,}\mathrm{au}$$ with
$M_{\mathrm{P}}=4.5\,\mathrm{M}_{\mathrm{Jup}}$, and also this model
reproduces the observed dust ring in both Band 7 and Band 4 (see Fig. 13).

Thus, we followed another approach to give an estimate of the possible maximum
planet mass. Lodato et al. (2019) report an empirical relation according to
which the width of an observed dust gap scales with the planet Hill radius
(see Eq. 4) in case of low disk viscosity ($\alpha_{\mathrm{SS}}\lesssim
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_{\mathrm{H}}\,,$ (6)

where $\Delta$ is defined as the distance between the minimum intensity in the
gap and the peak intensity in the ring. Hence, we take the value of $2\Delta$
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\text{\,}\mathrm{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\text{\,}\mathrm{M}_{\odot}$ and adopted a planet
location between $9.5\text{\,}\mathrm{au}$, consistent with the values in
Table 2, and $7.5\text{\,}\mathrm{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\,\mathrm{M}_{\mathrm{Jup}}$
planet (Fig. 13). Employing Eq. 6 with these values, we obtained the total
dust gap width $2\Delta$ 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
${\sim}4{-}8\,\mathrm{M}_{\mathrm{Jup}}$. A more massive planet would carve a
wider gap, either generating the dust ring at a greater distance from the star
or depleting the inner disk.

###  4.6 Comparison with previous estimates of the planet mass

Assuming a stellar mass of $0.1\,\mathrm{M}_{\odot}$, a disk viscosity
$\alpha_{\mathrm{SS}}=10^{-3}$, and that the gap is carved by a single planet
located at the minimum of the gap (${\sim}$8\text{\,}\mathrm{au}$$), Pinilla
et al. (2021) employed the Crida et al. (2006) criterion and estimated that
the minimum planet mass to open a gap in the gas is
${\sim}0.4\,\mathrm{M_{\mathrm{Jup}}}$. However, such a planet would generate
a pressure maximum at ${\sim}$13\text{\,}\mathrm{au}$$. So, either a more
massive planet or a higher distance from the star is needed to retrieve the
observed dust ring at ${\sim}$20\text{\,}\mathrm{au}$$. Following this
argument, the authors computed dust evolution simulations considering two
cases: a $2.4\,\mathrm{M}_{\mathrm{Jup}}$ planet located at
$8\text{\,}\mathrm{au}$ and a $0.5\,\mathrm{M}_{\mathrm{Jup}}$ planet at
$12\text{\,}\mathrm{au}$, assuming a gap shape obtained analytically. Their 1D
models evolved for ${10}^{6}\text{\,}\mathrm{yrs}$, 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\,\mathrm{M}_{\mathrm{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\,\mathrm{M}_{\mathrm{Jup}}$ planet case, instead, was able to reproduce
the contrast between the inner disk and the ring after
${10}^{5}\text{\,}\mathrm{yrs}$ 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}^{6}\text{\,}\mathrm{yrs}$. 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 (${\sim}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\,\mathrm{M}_{\odot}$, whereas Pinilla et al. (2021) assumed a
$0.1\,\mathrm{M}_{\odot}$ central star. This difference in the stellar mass,
along with a disk viscosity $\alpha_{\mathrm{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
${\sim}0.5\,\mathrm{M}_{\mathrm{Jup}}$ while our models require a minimum
planet mass of ${\sim}1.4\,\mathrm{M}_{\mathrm{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\,\mathrm{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}^{6}\text{\,}\mathrm{yrs}$$, whereas planet formation via
gravitational instabilities should occur at earlier times (Kratter & Lodato
2016). The properties of disks in the early Class 0/I protostars is an active
area of research and the properties of such disks are highly uncertain and
debated at the moment, both on observational and modeling grounds (Maury et
al. 2019; Tychoniec et al. 2020; Lebreuilly et al. 2021). Nevertheless, it
seems likely that at least some young disks will be prone to undergo
gravitational instabilities and possibly form massive planets early on through
this path.

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

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

##  5 Conclusions

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

  * •

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

  * •

Our numerical models cannot directly constrain the maximum planet mass
consistent with observations. Thus, we used an empirical relation between the
observed dust gap width and the Hill radius of the planet and estimate a
maximum planet mass of ${\sim}4{-}8\,\mathrm{M}_{\mathrm{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
${\sim}$7\text{\times}{10}^{-6}\text{\,}\mathrm{M}_{\odot}$$, corresponding to
a dust-to-gas mass ratio of ${\sim}${10}^{-2}\text{\,}$$.

  * •

The observed 12CO and 13CO channel maps are strongly affected by cloud
absorption. Nonetheless, our model with a final planet mass of
$2.0\,\mathrm{M}_{\mathrm{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
${\sim}{10^{-2}}$, a 12CO abundance of
$[^{12}\mathrm{CO}/\mathrm{H}_{2}]=5\times 10^{-5}$, and an isotopologue
abundance ratio $[^{12}\mathrm{CO}/^{13}\mathrm{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\,\mathrm{M}_{\odot}$ and systemic velocity of
$6.25\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ are solid estimates for
the actual properties of CIDA 1.

  * •

In the case of our simulation with a final planet mass of
$2.0\,\mathrm{M}_{\mathrm{Jup}}$, we are able to reproduce the low spectral
index (${\sim}2$) observed at the location of the dust ring
($10{-}30\,\mathrm{au}$). The simulated dust emission in both Band 7 and Band
4 is optically thick in the range $14{-}$20\text{\,}\mathrm{au}$$, with a
total fraction of optically thick emission corresponding to ${\sim}40{-}50\%$.
However, the emission from our model is optically thin between $20$ and
$30\text{\,}\mathrm{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\,\mathrm{M}_{\mathrm{Jup}}$ planet creates a velocity perturbation
that is predicted to be extremely arduous to observe. It would require 13CO
channel maps at a very high angular resolution of at least $0.02\arcsec$, with
a velocity resolution of $0.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$.

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 C$52$I$13000140001$), by the Deutsche Forschungs-gemeinschaft
(DFG, German Research Foundation) - Ref no. 325594231 FOR $2634$/$1$ TE
$1024$/$1$-$1$, and by the DFG cluster of excellence Origins (www.origins-
cluster.de). This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No. 823823 (RISE DUSTBUSTERS project) and from the
European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant
855130). This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2018.1.00536.S, ADS/JAO.ALMA#2015.1.00934.S, and
ADS/JAO.ALMA#2016.1.01511.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. M.T. has been supported by the UK Science and Technology research
Council (STFC) via the consolidated grant ST/S000623/1, and by the European
Union’s Horizon 2020 research and innovation programme under the Marie
Sklodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project).

## References

  * Anglada-Escudé et al. (2016) Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Ballabio et al. (2018) Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766
  * Barenfeld et al. (2016) Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142
  * Beckwith et al. (1990) Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924
  * Bergin et al. (2007) Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 751
  * Bitsch et al. (2015) Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112
  * Bohlin et al. (1978) Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132
  * Bollati et al. (2021) Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, MNRAS, 504, 5444
  * Briceno et al. (1993) Briceno, C., Calvet, N., Gomez, M., et al. 1993, PASP, 105, 686
  * Bruderer (2013) Bruderer, S. 2013, A&A, 559, A46
  * Bruderer et al. (2012) Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91
  * Chauvin et al. (2005) Chauvin, G., Lagrange, A. M., Dumas, C., et al. 2005, A&A, 438, L25
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698
  * Clarke et al. (2018) Clarke, C. J., Tazzari, M., Juhasz, A., et al. 2018, ApJ, 866, L6
  * Comeron et al. (1998) Comeron, F., Rieke, G. H., Claes, P., Torra, J., & Laureijs, R. J. 1998, A&A, 335, 522
  * Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
  * Cuello et al. (2019) Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114
  * Damasso et al. (2020) Damasso, M., Del Sordo, F., Anglada-Escudé, G., et al. 2020, Science Advances, 6, eaax7467
  * Davis et al. (2010) Davis, C. J., Chrysostomou, A., Hatchell, J., et al. 2010, MNRAS, 405, 759
  * Dipierro et al. (2015) Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73
  * Dodson-Robinson & Salyk (2011) Dodson-Robinson, S. E. & Salyk, C. 2011, ApJ, 738, 131
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1026
  * Draine (2006) Draine, B. T. 2006, ApJ, 636, 1114
  * Facchini et al. (2018) Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104
  * Faria et al. (2022) Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115
  * Fung & Chiang (2016) Fung, J. & Chiang, E. 2016, ApJ, 832, 105
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456
  * Haworth et al. (2020) Haworth, T. J., Cadman, J., Meru, F., et al. 2020, MNRAS, 494, 4130
  * Hendler et al. (2017) Hendler, N. P., Mulders, G. D., Pascucci, I., et al. 2017, ApJ, 841, 116
  * Hutchison et al. (2018) Hutchison, M., Price, D. J., & Laibe, G. 2018, MNRAS, 476, 2186
  * Izquierdo et al. (2018) Izquierdo, A. F., Galván-Madrid, R., Maud, L. T., et al. 2018, MNRAS, 478, 2505
  * Izquierdo et al. (2021a) Izquierdo, A. F., Smith, R. J., Glover, S. C. O., et al. 2021a, MNRAS, 500, 5268
  * Izquierdo et al. (2021b) Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021b, A&A, 650, A179
  * Kratter & Lodato (2016) Kratter, K. & Lodato, G. 2016, ARA&A, 54, 271
  * Kurtovic et al. (2021) Kurtovic, N. T., Pinilla, P., Long, F., et al. 2021, A&A, 645, A139
  * Laibe & Price (2014) Laibe, G. & Price, D. J. 2014, MNRAS, 440, 2136
  * Lambrechts & Johansen (2012) Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32
  * Lebreuilly et al. (2021) Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10
  * Liebert & Probst (1987) Liebert, J. & Probst, R. G. 1987, ARA&A, 25, 473
  * Liu et al. (2020) Liu, B., Lambrechts, M., Johansen, A., Pascucci, I., & Henning, T. 2020, A&A, 638, A88
  * Liu et al. (2019) Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75
  * Lodato et al. (2005) Lodato, G., Delgado-Donate, E., & Clarke, C. J. 2005, MNRAS, 364, L91
  * Lodato et al. (2019) Lodato, G., Dipierro, G., Ragusa, E., et al. 2019, MNRAS, 486, 453
  * Lodato & Price (2010) Lodato, G. & Price, D. J. 2010, MNRAS, 405, 1212
  * Manara et al. (2018) Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Maury et al. (2019) Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Mercer & Stamatellos (2020) Mercer, A. & Stamatellos, D. 2020, A&A, 633, A116
  * Miguel et al. (2020) Miguel, Y., Cridland, A., Ormel, C. W., Fortney, J. J., & Ida, S. 2020, MNRAS, 491, 1998
  * Mohanty et al. (2005) Mohanty, S., Jayawardhana, R., & Basri, G. 2005, ApJ, 626, 498
  * Morales et al. (2019) Morales, J. C., Mustill, A. J., Ribas, I., et al. 2019, Science, 365, 1441
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130
  * Mulders et al. (2021) Mulders, G. D., Pascucci, I., Ciesla, F. J., & Fernandes, R. B. 2021, ApJ, 920, 66
  * Natta & Testi (2001) Natta, A. & Testi, L. 2001, A&A, 376, L22
  * Natta et al. (2002) Natta, A., Testi, L., Comerón, F., et al. 2002, A&A, 393, 597
  * Okabe et al. (2000) 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 et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125
  * Payne & Lodato (2007) Payne, M. J. & Lodato, G. 2007, MNRAS, 381, 1597
  * Pinilla et al. (2012a) Pinilla, P., Benisty, M., & Birnstiel, T. 2012a, A&A, 545, A81
  * Pinilla et al. (2013) Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95
  * Pinilla et al. (2012b) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012b, A&A, 538, A114
  * Pinilla et al. (2021) Pinilla, P., Kurtovic, N. T., Benisty, M., et al. 2021, A&A, 649, A122
  * Pinilla et al. (2018) Pinilla, P., Natta, A., Manara, C. F., et al. 2018, A&A, 615, A95
  * Pinilla et al. (2017) Pinilla, P., Quiroga-Nuñez, L. H., Benisty, M., et al. 2017, ApJ, 846, 70
  * Pinte et al. (2020) Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Price (2007) Price, D. J. 2007, PASA, 24, 159
  * Price et al. (2018a) Price, D. J., Cuello, N., Pinte, C., et al. 2018a, MNRAS, 477, 1270
  * Price et al. (2018b) Price, D. J., Wurster, J., Tricco, T. S., et al. 2018b, PASA, 35, e031
  * Rabago & Zhu (2021) Rabago, I. & Zhu, Z. 2021, MNRAS, 502, 5325
  * Ragusa et al. (2017) Ragusa, E., Dipierro, G., Lodato, G., Laibe, G., & Price, D. J. 2017, MNRAS, 464, 1449
  * Reissl et al. (2016) Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87
  * Ricci et al. (2010) Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15
  * Ricci et al. (2014) Ricci, L., Testi, L., Natta, A., et al. 2014, ApJ, 791, 20
  * Rilinger et al. (2019) Rilinger, A. M., Espaillat, C. C., & Macías, E. 2019, ApJ, 878, 103
  * Rosotti et al. (2016) Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790
  * Sanchis et al. (2020) Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114
  * Schaefer et al. (2009) Schaefer, G. H., Dutrey, A., Guilloteau, S., Simon, M., & White, R. J. 2009, ApJ, 701, 698
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Tazzari et al. (2018) Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527
  * Tazzari et al. (2021) Tazzari, M., Clarke, C. J., Testi, L., et al. 2021, MNRAS, 506, 2804
  * Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339
  * Testi et al. (2016) Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111
  * Toci et al. (2020a) Toci, C., Lodato, G., Christiaens, V., et al. 2020a, MNRAS, 499, 2015
  * Toci et al. (2020b) Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020b, ApJ, 888, L4
  * Tychoniec et al. (2020) Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19
  * Ubeira Gabellini et al. (2019) Ubeira Gabellini, M. G., Miotello, A., Facchini, S., et al. 2019, MNRAS, 486, 4638
  * Veronesi et al. (2020) Veronesi, B., Ragusa, E., Lodato, G., et al. 2020, MNRAS, 495, 1913
  * Villenave et al. (2021) Villenave, M., Ménard, F., Dent, W. R. F., et al. 2021, A&A, 653, A46
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323
  * Warren & Brandt (2008) Warren, S. G. & Brandt, R. E. 2008, Journal of Geophysical Research (Atmospheres), 113, D14220
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Williams & Best (2014) Williams, J. P. & Best, W. M. J. 2014, ApJ, 788, 59
  * Williams et al. (2019) Williams, J. P., Cieza, L., Hales, A., et al. 2019, ApJ, 875, L9
  * Wilson & Rood (1994) Wilson, T. L. & Rood, R. 1994, ARA&A, 32, 191
  * Zhu et al. (2018) Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850
  * Zubko et al. (1996) Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

##  Appendix A Sound speed power law exponent from the SED

Figure 15: 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\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, between the two vertical
dashed gray lines.

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_{\mathrm{B}}T_{\mathrm{out}}\ll h\nu\ll k_{\mathrm{B}}T_{\mathrm{in}}$
(where $k_{\mathrm{B}}$ is the Boltzmann constant, $T_{\mathrm{in}}$ and
$T_{\mathrm{out}}$ are the disk temperatures at the inner and outer radii of
the disk, respectively), usually corresponding to a wavelength range roughly
between $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, the slope of the SED is connected
to the exponent $q$ of the disk sound speed profile by this relation:

$\frac{\mathrm{d}\mathrm{\log_{10}}(\lambda\,F_{\lambda})}{\mathrm{d}\mathrm{\log_{10}(\lambda)}}=-4+\frac{1}{q}.$
(7)

The measured SED of CIDA 1 is reported in Hendler et al. (2017). We performed
a linear fit of the value of $\mathrm{\log_{10}}(\lambda\,F_{\lambda})$ with
respect to $\log_{10}(\lambda)$ in the wavelength range
$2.5{-}$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (Fig. 15), obtaining a
slope of ${\approx}-0.62$. Using the (7), we derive that $q\approx 0.3$.

##  Appendix B Planet accretion and migration

Figure 16: Time evolution of planet-star distance (top panel) and planet mass
(bottom panel) from the simulation with an initial planet mass of
$1.5\,\mathrm{M}_{\mathrm{Jup}}$.

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. 16, 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\,\mathrm{M}_{\mathrm{Jup}}$. In both
plots, a fast transient is evident in the first
${10}^{4}\text{\,}\mathrm{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\text{\times}{10}^{4}\text{\,}\mathrm{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. 16, a spread in the planet-star distance is
visible, indicating an elliptical orbit of the planet in our simulation.
Knowing the periastron $r_{\mathrm{p}}$ and apoastron $r_{\mathrm{a}}$, we can
calculate the planet orbital eccentricity:

$e=\frac{r_{\mathrm{a}}-r_{\mathrm{p}}}{r_{\mathrm{a}}+r_{\mathrm{p}}}\,,$ (8)

obtaining a value $e\approx 0.01$ after
$8\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$.

##  Appendix C Dust mass rescaling

Figure 17: 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\,\mathrm{M}_{\mathrm{Jup}}$ after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$.

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. 17. The outer ellipse is
positioned in the image center, has a semimajor axis of $0.27\arcsec$, a
semiminor axis of $0.22\arcsec$, and an inclination of $11^{\circ}$; for the
inner ellipse, we only change the semiaxis length to $0.13\arcsec$ and
$0.10\arcsec$. The flux contained in this elliptical annulus is
$29.0\text{\,}\mathrm{mJy}$. In our synthetic images, without an inner disk,
we consider the flux contained within the red dashed ellipse in the bottom
panel of Fig. 17, 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\text{\,}\mathrm{mJy}$.

##  Appendix D Dust emission profiles at longer evolution times

Fig. 18 shows the dust emission radial profiles in Band 7 and Band 4 of each
simulation after an evolution time of
$8\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$. Compared to the results after
$4\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$ (Fig. 5), dust grains had more
time to undergo radial drift and possibly accrete onto the star.

Figure 18: As in Fig. 5, but for the case of simulations evolved for
$8\text{\times}{10}^{4}\text{\,}\mathrm{yrs}$.

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

##  Appendix E Cloud absorption

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

We aim here to qualitatively estimate the cloud absorption affecting the
observed gas channel maps. Davis et al. (2010) mapped the 12CO 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 19 shows the 12CO 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 13CO $J=3{-}2$
emission. CIDA 1 is located in the subregion L1495 W. Since in this area the
13CO data are missing, we consider the nearby subregion L1595 E for our
estimates.

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

${\dfrac{\int T^{*}_{\mathrm{A}}(^{12}\mathrm{CO})\,\mathrm{d}\varv}{\int
T^{*}_{\mathrm{A}}(^{13}\mathrm{CO})\,\mathrm{d}\varv}=\dfrac{1-\mathrm{e^{-\tau_{12}}}}{1-\mathrm{e^{-\tau_{13}}}}}\,,$
(9)

where $\tau_{13}=\tau_{12}/70$ is the 13CO optical depth. We integrate the
antenna temperatures along velocity intervals with the same centers and width
of the observed 12CO channel maps (top panels in Fig. 6). The resulting 12CO
optical depth is presented in Fig. 21. The morphology of the optical depth
spectrum follows approximately the cloud absorption behavior on the 12CO
observation, being strongest at
${\sim}6.0{-}$7.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ while fading
at ${\sim}$4.0\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$ and
${\sim}$8.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$, 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 12CO 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 20: 12CO ($J=3{-}2$) and 13CO ($J=3{-}2$) spectra with respect to the
local standard of rest (LSR) in the Taurus subregion L1495 E. The spectra are
averaged over the entire extent of the subregion. Data are from Davis et al.
(2010).  Figure 21: 12CO ($J=3{-}2$) optical depth in the Taurus subregion
L1495 E, calculated from the spectra in Fig 20 using Eq. 9. We focus only on
the velocity interval where line emission is detected. Bins replicate the
central velocity and resolution of the observed 12CO channel maps.