$\langle$reception date$\rangle$ $\langle$acception date$\rangle$
$\langle$publication date$\rangle$

protoplanetary disk – planet-disk interaction – submillimeter: planetary
systems

# Kinematic signatures of a  low-mass planet with a moderately inclined orbit
in a protoplanetary disk

Kazuhiro D. Kanagawa11affiliation: College of Science, Ibaraki University,
2-1-1 Bunkyo, Mito, Ibaraki 310-8512, Japan Tomohiro Ono22affiliation: School
of Natural Sciences, Institute for Advanced Study, Princeton, NJ 08544, USA
Munetake Momose11affiliationmark<EMAIL_ADDRESS>
###### Abstract

A planet embedded in a protoplanetary disk produces a gap by disk–planet
interaction. It also generates velocity perturbation of gas, which can also be
observed as deviations from the Keplerian rotation in the channel map of
molecular line emission, called kinematic planetary features. These observed
signatures provide clues to determine the mass of the planet. We investigated
the features induced by the planet with an inclined orbit through three-
dimensional hydrodynamic simulations. We found that a smaller planet, with the
inclination being $\sim 10^{\circ}$ – $20^{\circ}$, can produce kinematic
features as prominent as those induced by the massive coplanar planet. Despite
the kinematic features being similar, the gap is shallower and narrower as
compared with the case in which the kinematic features are formed by the
coplanar planet. We also found that the kinematic features induced by the
inclined planet were fainter for rarer CO isotopologues because the velocity
perturbation is weaker at the position closer to the midplane, which was
different in the case with a coplanar massive planet. This dependence on the
isotopologues is distinguished if the planet has the inclined orbit. We
discussed two observed kinematic features in the disk of HD 163296. We
concluded that the  kink observed at 220 au can be induced by the inclined
planet, while the kink at 67 au is consistent to that induced by the coplanar
planet.

##  1 Introduction

A planet is formed in a protoplanetary disk and gravitationally interacts with
the surrounding gas via spiral waves (Lin & Papaloizou, 1979; Goldreich &
Tremaine, 1980). When a planet is sufficiently large, a gas depleted region is
formed along with the planet orbit, called a planetary gap (e.g.,Lin &
Papaloizou (1986), Kley (1999), Crida et al. (2006), Duffell & MacFadyen
(2013), Fung et al. (2014), Kanagawa et al. (2015)). Recent observations with
Atacama Large Millimeter/submillimeter Array (ALMA) have revealed a number of
protoplanetary disks with gap structures (e.g., ALMA Partnership et al.
(2015), Tsukagoshi et al. (2016), Long et al. (2018), Andrews et al. (2018),
van der Marel et al. (2019),  Cieza et al. (2019)). The planetary gap is one
of the convincing origins to explain these observed gap structures, although
other processes (e.g., the dust evolution in the vicinity of the snow line
(Zhang et al., 2015; Okuzumi et al., 2016), secular gravitational
instabilities (Takahashi & Inutsuka, 2014, 2016), disk-wind (Takahashi & Muto,
2018)) are proposed. Moreover, deviations from the Keplerian rotation have
been detected in the channel maps of 12CO emissions in some protoplanetary
disks (e.g., Pinte et al. (2018), Teague et al. (2018), Casassus & Pérez
(2019), Teague et al. (2019), Pinte et al. (2019) , Pinte et al. (2020), Bae
et al. (2022), Izquierdo et al. (2022); Calcino et al. (2022), and see also
Pinte et al. (2022)). These anomalies from the Keplerian rotation are thought
to be originated from the velocity perturbation induced by the embedded planet
(Perez et al., 2015; Pinte et al., 2018) and density perturbation (spiral and
gap) in addition to the velocity perturbations (Pérez et al., 2018; Rabago &
Zhu, 2021).

Particularly, Pinte et al. (2020) found the kinematic planetary features (or
called velocity kinks) that can be generated from the planet within the
observed gap, in nine protoplanetary disks in the DSHARP sample (Andrews et
al., 2018), which strongly suggests that the observed gap is induced by the
planet. However, they also highlighted that the planet mass suggested by the
kinematic planetary feature was 4–10 times larger than that estimated by the
gap shape. To address this discrepancy of the masses predicted by the gap and
kinematic planetary feature, we consider a planet with an inclined orbit in
this paper.

Exoplanets with misalignment of the angular momentum vector with the stellar
spin axis have been discovered (Winn & Fabrycky, 2015). This misalignment can
indicate that the planet has an inclined orbit against the protoplanetary disk
due to, for instance, the planet– scattering in the protoplanetary disk
(e.g.,Chatterjee et al. (2008)). The disk–planet interaction for the low-mass
planet with low-inclination has been investigated by using a linear theory
(Tanaka & Ward, 2004) and hydrodynamic simulations (e.g., Arzamasskiy et al.
(2018)). In the case with the massive planet, it has been investigated by
hydrodynamic simulations (e.g., Cresswell et al. (2007), Marzari & Nelson
(2009), Bitsch & Kley (2011), Xiang-Gruess & Papaloizou (2013), Kloster &
Flock (2019), Zhu (2019)). With the orbital inclination, a larger mass of a
planet is required to form a gap (Bitsch et al., 2013). Moreover, the planet
could disturb gas at a high altitude compared with the case with the coplanar
planet, which may appear in the structure of the channel map as the kinematic
planetary features.

In this paper, we investigate the kinematic planetary features formed by the
low-mass planet with moderate orbital inclination using three-dimensional (3D)
hydrodynamic simulations and radiative transfer calculations. We present our
model and numerical method in Section 2. In Section 3, we present our results
and discuss the implications of our results to observations and possible
scenarios to maintain the orbital inclination in the disk phase in Section 4.
Finally, the summary of our results is provided in Section 5.

##  2 Our model and method

###  2.1 Basic equations

We considered 3D hydrodynamics of a protoplanetary disk with an inclined-
orbital planet, assuming a non-self-gravitational and locally isothermal disk.
Hence, the basic equations are expressed as

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\bm{v}\right)=0$ (1)
$\displaystyle\frac{\partial\rho\bm{v}}{\partial
t}+\nabla\cdot\left(\rho\bm{v}+P{\rm\bm{I}}\right)=-\rho\nabla\Phi+\nabla\cdot\bm{\Pi},$
(2)

where $\rho$ and $\bm{v}$ are the gas density and velocity, respectively. We
used the spherical coordinates, $r$ is the spherical radius, $\theta$ and
$\phi$ are the polar and azimuth angles, respectively, and hence $\bm{r}$ =
(r,$\theta$,$\phi$) and $\bm{v}$=($v_{r}$,$v_{\theta}$,$v_{\phi}$). The origin
of the coordinate is set to be the position of the central star. For
simplicity, note that we also used the cylindrical coordinates (R, $\phi$, z),
where $R=r\sin(\theta)$ and $z=r\cos(\theta)$. The unit tensor is denoted by
${\rm\bm{I}}$. Since we adopted the simple locally isothermal equation of
state, the gas pressure $P$ is expressed as

$\displaystyle P$ $\displaystyle=c_{s}^{2}\rho,$ (3)

where $c_{s}$ is the sonic speed, which is given by a function of $R$. The
gravitational potential $\Phi$ is expressed as

$\displaystyle\Phi$
$\displaystyle=-\frac{GM_{\ast}}{\left|\bm{r}\right|}-\frac{GM_{\rm
p}}{\sqrt{\left(\left|\bm{r}-\bm{r}_{\rm
p}\right|\right)^{2}+\epsilon^{2}}}+\frac{GM_{\rm p}}{\left|\bm{r}_{\rm
p}\right|^{3}}\bm{r}\cdot\bm{r}_{\rm p},$ (4)

where $G$ is the gravitational constant and $M_{\ast}$ and $M_{\rm p}$ are the
masses of the central star and planet, respectively. The plant is located at
$\bm{r}=\bm{r}_{\rm p}$. In the following, subscript p denotes the quantities
of the planet. The smoothing parameter of planet gravity $\epsilon$ is set as
$0.1r_{H}$, where $r_{H}$ is the hill radius of the planet defined by $a_{\rm
p}\left[M_{\rm p}/(3M_{\ast})\right]^{1/3}$ in which $a_{\rm p}$ is the
semimajor axis of the planet. The viscous stress tensor $\bm{\Pi}$ is

$\displaystyle\bm{\Pi}$
$\displaystyle=\rho\nu\left[\nabla\bm{v}+\left(\nabla\bm{v}\right)^{{\dagger}}-\frac{2}{3}{\rm\bm{I}}\nabla\cdot\bm{v}\right],$
(5)

where $\nu$ is the kinetic viscosity. We adopted the $\alpha$-viscosity
prescription (Shakura & Sunyaev (1973)). Hence,

$\displaystyle\nu$ $\displaystyle=\alpha c_{s}h,$ (6)

where $h$ is the scale height ($=c_{s}/\Omega_{K}$) and $\Omega_{K}$ is the
Keplerian angular velocity defined by $\sqrt{GM_{\ast}/R^{3}}$, and we use
$\alpha=10^{-3}$ in this paper . For convenience, we also defined
$v_{K}=R\Omega_{K}$ as the Keplerian rotation velocity.

###  2.2 Initial and boundary conditions

We assumed a power-law distribution of the initial surface density $\Sigma$,
which is expressed as

$\displaystyle\Sigma$
$\displaystyle=\Sigma_{0}\left(\frac{R}{R_{0}}\right)^{p_{s}},$ (7)

where $R_{0}$ is the unit of the radius. We adopted $p_{s}=-1/2$ and
$\Sigma_{0}=1$ in the following. Note that the choice of $\Sigma_{0}$ does not
affect the result of the hydrodynamic simulation, except for the absolute
value of the density, since we do not consider the self-gravity of the gas
disk. We also used $M_{\ast}$ as the unit of the mass and
$\Omega_{K}^{-1}(R_{0})$ as the unit of time; hence, $\Sigma$ and $\rho$ are
normalized by $M_{\ast}/R_{0}^{2}$ and $M_{\ast}/R_{0}^{3}$, respectively. We
also assumed the power-law distribution of the sound speed and vertically
isothermal distribution, which is expressed as

$\displaystyle c_{s}$
$\displaystyle=c_{s}{}_{0}\left(\frac{R}{R_{0}}\right)^{q},$ (8)

where we set $q=-1/2$ in this paper. The aspect ratio of the disk, $H$, is
written by $c_{s}/v_{K}$, and it is constant as
$H_{0}=c_{s}{}_{0}/v_{K}(R_{0})$ throughout the computational domain. We
investigated the cases with $H_{0}=0.07$ and $H_{0}=0.1$, referring to the
recent observations of the outer region of the planetary disks (e.g., Rich et
al. (2021), Öberg et al. (2021)). Assuming hydrostatic equilibrium, we adopted
the Gaussian profile in vertical direction, which is expressed as

$\displaystyle\rho$
$\displaystyle=\frac{\Sigma_{0}}{\sqrt{2\pi}h{}_{0}}\left(\frac{R}{R_{0}}\right)^{p}\exp\left[\frac{GM_{\ast}}{c_{s}^{2}}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$
(9)

where $h{}_{0}$ is the scale height at $R_{0}$ and $p=p_{s}-(2q+3)/2$ is equal
as $-1.5$ in this paper.

The initial azimuthal velocity is given by hydrodynamic equilibrium in $R$
direction; that is,

$\displaystyle\Omega$
$\displaystyle=\Omega_{K}\left[\left(1+q\right)+\left(p+q\right)H^{2}-q\frac{R}{r}\right].$
(10)

The viscous drift velocity is given by

$\displaystyle v_{\rm vis}$
$\displaystyle=-\frac{3\nu}{R}\left[p+\frac{4}{3}q+\frac{10q+9}{3}\frac{GM_{\ast}}{c_{s}^{2}}\left(\frac{1}{r}-\frac{1}{R}\right)\right].$
(11)

We assumed vertical hydrostatic equilibrium $v_{z}=0$. Therefore,
$v_{r}=v_{\rm vis}\sin(\theta)$ and $v_{\theta}=v_{\rm vis}\cos(\theta)$ in
the initial profile.

At the inner and outer boundaries in the $r$ direction, we fixed the physical
quantities on the initial values. At the upper and lower boundaries in the
$\theta$ direction, we adopted a vertically outgoing boundary for $v_{r}$ and
$v_{\theta}$. That is, the velocities are copied just inside of the boundary
if the gas flows outside of the computational domain ($v_{z}>0$ and $z>0$, or
$v_{z}<0$ and $z<0$); otherwise, $v_{r}$ and $v_{\theta}$ are set to zero. The
azimuthal velocity is copied from that of the cell just next to the boundary.
Note that we applied no wave-damping zone for simplicity.

###  2.3 Numerical methods

####  2.3.1 Hydrodynamic simulations

We used Athena++111https://www.athena-astro.app (Stone et al., 2020) to solve
basic equations (1) and (2). Athena++ is a grid-based magnetohydrodynamic code
using a higher-order Godunov scheme with the orbital advection scheme (Masset
(2000), see Ono et al. in preparation) and static/adaptive mesh refinement.

The computational domain runs from $r=0.4R_{0}$ to $r=3.0R_{0}$ in
$r$-direction and covers the region of $-5.6H_{0}<\theta-\pi/2<5.6H_{0}$ and
$-\pi<\phi<\pi$. We logarithmically divided the computational domain into 192
meshes in the $r$-direction and equally divided it by 80 meshes in the
$\theta$-direction and 640 meshes in the $\phi$ direction. In addition, we
used one-level static mesh refinement in the range of $0.6R_{0}<r<1.4R_{0}$,
and the root grid was divided into two in the refinement region. The grid size
is approximately $5\times 10^{-3}$ in the region where the refinement is
adopted. When $H_{0}=0.07$, which is our fiducial value, the scale height is
divided into 14 meshes in the refinement region, while it is divided in 7
meshes in the other region.

We adopted the second-order piecewise linear method for the spatial
reconstruction method (it is the default of ATHENA++) and third-order accurate
Strong Stability Preserving (SSP) variant method for time-integration with the
CFL coefficient being 0.6. Note that we confirmed that the result does not
change if we chose the second-order van Leer scheme (the default of ATHENA++)
with CFL coefficient being 0.3.

####  2.3.2 Planet

We considered a planet with an orbital inclination. For simplicity, we ignored
a torque from the surrounding gas and planetary migration. We only considered
the gravity of the central star when calculating the planetary orbit. The
initial location of the planet was set by $\bm{r}_{\rm
p}=(R_{0},\pi/2-\sin\left(i_{\rm orb}\right)\sin\left(\pi/4\right),0)$, and
the initial velocity was set by $\bm{v}_{p}=(0,\sin\left(i_{\rm
orb}\right)\sin\left(\pi/4\right),v_{K}(R_{0}))$, where $i_{\rm orb}$ is the
orbital inclination of the planet222 The initial condition is a bit different
from the steady state condition, and hence the planet move slight from the
initial radius. However, the planet reaches steady state after just after a
few orbits. Hence, it could not affect our result significantly. . The time
integration for the orbital calculation was used by the second-order leapfrog
method to stably calculate the planetary orbit long time.

We increased the planet mass gradually in the beginning using

$\displaystyle M_{\rm p}$ $\displaystyle=M_{\rm
p}{}_{\rm,fin}\sin^{2}\left(\frac{\pi}{2}\frac{t}{\tau_{m}}\right),$ (12)

until $t=\tau_{m}$, where $\tau_{m}$ is the timescale of the mass growth and
we adopted $\tau_{m}=2\pi/\Omega_{K}(R_{0})$. After $t=\tau_{m}$, we fixed the
planetary mass on $M_{\rm p}{}_{\rm,fin}$. In the following, we referred
$M_{\rm p}{}_{\rm,fin}$ as $M_{\rm p}$.

####  2.3.3 Radiative transfer simulations

We used the density and velocity profiles given by the 3D hydrodynamic
simulations described above as an input model, but is slightly modified (see
below), to RADMC3D333https://www.ita.uni-heidelberg.de/
dullemond/software/radmc-3d (Dullemond et al., 2012) that is a publicly
available code for radiative transfer calculations for dust and molecules. In
the following, we used the results of hydrodynamic simulations at $t=100$
orbits, corresponding to $0.1$ Myr when $R_{0}$=100 au.

First, we calculated the gas temperature by carrying out thermal Monte Carlo
simulation by assuming that the gas temperature is equal to the dust
temperature. Then, we carried out molecular line transfer simulations to
produce synthetic channel maps using the temperature produced by the thermal
Monte Carlo. Note that the temperature adopted in the radiative transfer
simulations is not consistent to that adopted in the 3D hydrodynamic
simulations with the locally isothermal equation of state. This temperature
inconsistency may not affect our main result described in Section 3 (see more
discussion in Appendix A).

We assumed a T-Tauri star that is similar to a protostar as a central star
with $M_{\ast}=1M_{\odot}$, and its radius is $2.5R_{\odot}$. The effective
temperature of the star was set to be 4500 K; hence, its bolometric luminosity
was approximately $1.7L_{\odot}$. The star is located in the origin of the
coordinate the same as in the hydrodynamic simulation.

We assumed that the dust grains are well-mixed to gas; hence, for simplicity,
the scale height of the dust grains is the same as that of the gas. The size
distribution of the dust grains is assumed to be $s^{-3.5}$, where $s$ is the
size of the dust grain. The minimum and maximum sizes of the dust grain are
$0.1$ $\mu$m and $1$ mm, respectively. We used the same compositions of the
dust grains to that adopted in DSHARP (Birnstiel et al., 2018), and the
absorption and scattering coefficients are averaged by the size distribution.
We assumed that the gas-to-dust mass ratio is 100, and it is constant in the
entire domain.

In this paper, we considered the rotation transition in $J=2-1$ of 12CO and
its isotopologues 13CO and C18O. The fractional abundance of 12CO relative to
H2 was set to $9.83\times 10^{-5}$, and the abundances of 13CO and C18O were
$1.49\times 10^{-6}$ and $2.26\times 10^{-7}$, respectively. The line transfer
was done in local thermal equilibrium using molecular data from the LAMBDA
database444https://home.strw.leidenuniv.nl/ moldata/ (Schöier et al., 2005;
van der Tak et al., 2020).

We adopted $10^{7}$ photons for the thermal Monte Carlo simulations and
$3\times 10^{6}$ photons for imaging of dust continuum and molecular line
emissions, respectively. The density and velocity structures in the radiative
transfer were fed from the 3D hydrodynamic simulation with $R_{0}=100$ au.
However, the profile has a relatively large cavity because of the inner
boundary of the computational domain (its size was 40 au when $R_{0}=100$ au).
This unintentional cavity can slightly affect the temperature and channel map.
With this, we inserted a Keplerian rotation disk profile to extend the initial
profile from the inner boundary of the hydrodynamic simulation to $10$ au
(=$0.1R_{0}$). We also should emphasis that the structure close to the inner
boundary of the simulation, that is $\sim 40$ au, can be affected by the
boundary condition of the simulation. To address the inner structure, e.g.,
inner secondary gap (e.g., Dong et al. (2017),Bae et al. (2017)), a simulation
with the inner boundary at smaller radii is required. These structures are out
of the scope of this paper.

We assumed that the distance from the star was 140 pc and the inclination and
position angles were $45^{\circ}$ and $135^{\circ}$, respectively. For the
initial disk mass, we assumed to be $0.1M_{\odot}$. The planet is located at
($x_{\rm p},y_{\rm p}$) = (-0.06 arcsec, 0.56 arcsec) or (-8.4 au,78.4 au )in
the images.

##  3 Results

###  3.1 CO emission layers

Before showing the results with the planet, we presented a height of the
emission layer of CO isotopologues in the disk with no planet because it
reveals the region that mainly contributes in the determination of the
structure of the channel map.

Figure 1:  Locations of $\tau=1$ for $J=$2-1 emissions of 12CO, 13CO, and
C18Ofor the unperturbed disk (solid lines) with $H_{0}=0.07$ (top) and
$H_{0}=0.1$ (bottom). The dotted lines represent the case with the planet (Run
3 in Table 3.2).

Figure 1 shows the emission layers of 12CO, 13CO, and C18O, for the cases with
$H_{0}=0.07$ and $0.1$. We defined the emission layer as a location of
$\tau=1$ in the face-on disk, which is calculated by RADMC3D, where $\tau$ is
the optical depth. For reference, we plotted the emission layers in the case
with the inclined planet (Run 3 in Table 3.2) in the figure. There is no
significant difference between the cases with/without planet because the
planet only forms a shallow gap (see next subsection). The emission layers
shift down close to the midplane when the planet forms a deep gap (Appendix
B).

The ratio of the height of the emission layer $z_{\tau=1}$ to $R$ is
approximately constant in the simulation domain (Figure 1). In the case with
$H_{0}=0.07$ and $R_{0}=100$ au, $z_{\tau=1}/R_{0}=0.29\simeq 4H_{0}$ for
12CO, $z_{\tau=1}/R_{0}=0.21\simeq 3H_{0}$ for 13CO, and
$z_{\tau=1}/R_{0}=0.11\simeq 1.6H_{0}$ for C18O. When $H_{0}=0.1$, the
emission layers for 12CO and 13CO are $4H_{0}$ and $3H_{0}$, respectively,
following the similar scaling as in the case with $H_{0}=0.07$. However, for
C18O, $z_{\tau=1}/R_{0}=2.3H_{0}$ when $H_{0}=0.1$, and it is slightly larger
than that in the case of $H_{0}=0.07$. In all cases, the emission layers are
always contained within the computational domain of hydrodynamic and radiative
transfer simulations.

###  3.2 Kinematic features in channel maps with varying planetary orbital
inclination

We simulated the disk structure with a planet with planet mass $M_{\rm p}$ and
orbital inclination $i_{\rm orb}$ (Table 3.2). We chose $H_{0}=0.07$ as our
fiducial value (Runs 1–6) and demonstrated the cases with $H_{0}=0.1$ for
$M_{\rm p}/M_{\ast}=5\times 10^{-4}$ (Runs 7-10) to see a dependence on the
aspect ratio. For reference, we also carried out the simulation with a massive
planet ($M_{\rm p}/M_{\ast}=2\times 10^{-3}$) with the coplanar orbit (Run
11). We used $\alpha=10^{-3}$ in all the cases.

Parameters  Run $M_{\rm p}/M_{\ast}$ $H_{0}$ $\alpha$ $i_{\rm orb}$ (degree) 1
$5\times 10^{-4}$ ($0.5M_{J}$)∗ $0.07$ $10^{-3}$ $0.0$ 2 $5\times 10^{-4}$
$0.07$ $10^{-3}$ $7.1$ 3 $5\times 10^{-4}$ $0.07$ $10^{-3}$ $14.1$ 4 $5\times
10^{-4}$ $0.07$ $10^{-3}$ $21.1$ 5 $3\times 10^{-4}$ ($0.3M_{J}$) $0.07$
$10^{-3}$ $14.1$ 6 $1\times 10^{-4}$ ($0.1M_{J}$) $0.07$ $10^{-3}$ $14.1$ 7
$5\times 10^{-4}$ ($0.5M_{J}$) $0.10$ $10^{-3}$ $0.0$ 8 $5\times 10^{-4}$
$0.10$ $10^{-3}$ $14.1$ 9 $5\times 10^{-4}$ $0.10$ $10^{-3}$ $21.1$ 10
$5\times 10^{-4}$ $0.10$ $10^{-3}$ $29.4$ 11 $2\times 10^{-3}$ ($2.0M_{J}$)
$0.07$ $10^{-3}$ $0.0$ ∗ Mass of the planet in Jupiter mass when
$M_{\ast}=1M_{\odot}$.

Figure 2:  The dust continuum images at 230 GHz (left) and the channel maps of
12CO J=2-1 emission for various $i_{\rm orb}$ in the cases with $M_{\rm
p}/M_{\ast}=5\times 10^{-4}$ (Runs 1–4 in Table 3.2). The uppermost panels are
the case of the simple Keplerian disk with no planet, and the subsequent lines
correspond to the cases with the planet with $i_{\rm orb}=0^{\circ}$,
$7.1^{\circ}$, $14.1^{\circ}$, and $21.1^{\circ}$. The white cross in the dust
continuum image represents the location of the planet. The disk rotation is
clockwise and the white arrow in the lower right corner of the continuum
images indicate the direction of the rotation. In the following figures, the
direction of the disk rotation is the same.

Figure 2 shows the dust continuum images at 230 GHz and the channel maps of
12CO J=2-1 emission when $M_{\rm p}/M_{\ast}=5\times 10^{-4}$ with various
$i_{\rm orb}$. The meridian angle of the planet varies with time when the
orbital inclination is nonzero. In the figure, the planet is located at the
midplane ($z_{\rm p}$=0). For reference, we show the case with no planet,
demonstrating the pattern of the purely Keplerian rotation disk (the effect of
the pressure gradient is included, as expressed in Equation 10), in which
there is no signature due to, for instance, vertical shear instability as
discussed by Barraza-Alfaro et al. (2021) because of the relatively large
viscosity we adopted and our coarse resolution in vertical direction (Flores-
Rivera et al., 2020) . When $i_{\rm orb}=0^{\circ}$, the relatively deep gap
is observed in the dust continuum, and the channel map is similar to that for
the Keplerian rotation disk with no planet, although the small deviations can
be observed in some points, particularly in $\Delta v=1$ km/s map. The gap
becomes shallower and narrower as $i_{\rm orb}$ increases. Moreover, different
from the case with $i_{\rm orb}=0^{\circ}$, significant deviations from the
Keplerian rotation are observed around the planet in the channel map in the
cases of $i_{\rm orb}=7^{\circ}$ and $i_{\rm orb}=14^{\circ}$. When the
inclination is large and $i_{\rm orb}=21^{\circ}$, the deviation in the
channel map becomes as weak as that in the case of $i_{\rm orb}=0^{\circ}$,
whereas almost no gap can be observed in the dust continuum.

Figure 3:  Density ($\rho$) and perturbed velocity ($v_{r}$, $v_{\theta}$,
$v_{\phi}-v_{\phi}{}_{,\rm Kep}$) distributions that are averaged in the range
of $\phi=\phi_{\rm p}\pm R_{H}/R_{0}$, which is the same cases shown in Figure
2, except for the Keplerian case with no planet. The horizontal dashed lines
indicate the emission layer of 12CO (see Figure 1). Note that the front layer
of Figure 1 corresponds the layer of $z<0$ since we adopted $i_{\rm
disk}=135^{\circ}$.

In Figure 3, we show the velocities that are averaged in the vicinity of the
planet location (in the range of $\phi=\phi_{\rm p}\pm R_{H}/R_{0}$) in the
cases shown in Figure 2. When $i_{\rm orb}=0^{\circ}$, the gas velocity around
the emission layer is small because it is distant from the planet. Similarly,
the velocity around the emission layer is not significant when $i_{\rm
orb}=21^{\circ}$. Meanwhile, with a moderate inclination, such as $i_{\rm
orb}=7^{\circ}$ and $14^{\circ}$, the velocities at a high altitude are
significantly disturbed by the planet. Particularly, the perturbation of
$v_{r}$ is prominent, and its amplitude is approximately 30 % of the Keplerian
rotation velocity. The disturbance of $v_{r}$ causes the deviation from the
Keplerian rotation appeared in the channel map shown in Figure 2 (see also
Appendix C). Note that the velocity distribution and relevant kinematic
features induced by the inclined planet vary in the orbital period of the
planet, which is further discussed in Section 3.5

Figure 4:  Azimuthal distribution of density and perturbed velocities around
$R=R_{p}$ (averaged in the range of $R_{p}\pm R_{H}$), in the same case as
shown in Figure 3.

Figure 4 shows the azimuthal distributions of density and perturbed velocities
around $R=R_{p}$, in the same case of Figure 3. As shown in Figures 3 and 4,
one finds corn-like shape of $\rho$ and $v_{r}$, $v_{\phi}$, which corresponds
to the mach corn induced by the planet. Hence, we conclude that the kinematic
features shown in Figure 2 are reflected from the planetary mach corn. The
corn is thinner with larger inclination angle and higher vertical velocity of
the planet, and it is almost vanished in the case of $i_{\rm
orb}=21^{{}^{\circ}}$.

###  3.3 Mass dependence and variations of kinematic features with CO isotopes

Figure 5:  Dust continuum at 230 GHz (left) and the channel maps for $\Delta
v=1$ km/s of 12CO, 13CO, and C18O for various planet masses. Except for the
panels in the bottom line, the orbital inclination is $14^{\circ}$. For
reference, the bottom line shows the case with the coplanar massive planet
($M_{\rm p}/M_{\ast}=2\times 10^{-3}$ and $i_{\rm orb}=0^{\circ}$). In all the
panels, the planet locates at the midplane.

Figure 5 shows the parts of the channel maps for the CO isotopologues (12CO,
13CO, and C18O) with various planet mass with $i_{\rm orb}=14^{\circ}$ (Runs
4, 5, and 6), except for the panels in the bottom line. In the bottom line of
the figure, we show the channel maps in the case with the massive planet in
the coplanar orbit ($M_{\rm p}/M_{\ast}=2\times 10^{-3}$ and $i_{\rm
orb}=0^{\circ}$, Run 11), for reference. Note that the planet locates at the
midplane in all cases shown in the figure, as the same as Figure 2

In the case of $M_{\rm p}/M_{\ast}=10^{-4}$, almost no gap is formed in the
dust continuum, and no feature of the planet is observed in all the channel
maps of the isotopologues. As the mass increases, the dust continuum gap is
gradually prominent, and some features associated to the planet can be
observed in the 12CO channel map. However, there is no significant feature in
the channel maps of other isotopologues, including 13CO and C18O. The emission
layers of 13CO and C18O are located closer to the midplane as compared to that
of 12CO, as presented in Section 3.1. Considering the velocity perturbation
shown in Figure 3, the gas disturbance is weaker as it is close to the
midplane. No significant deviation appears in the channel maps of those
molecules since the motion of gas near the emission layers of 13CO and C18O is
approximately Keplerian motion.

In contrast, in the case with the coplanar massive planet shown in the bottom
line of Figure 5, the deep and wide gap is observed in the dust continuum. In
the channel map of 12CO, the significant anomaly from the Keplerian rotation
can be found. In the channel map of 13CO, the emission in the vicinity of the
planet becomes a bit fainter because of gap formation, and the neck-like
structure can be observed. In the C18O map, the emission around the planet is
much fainter than other parts, and then the gap structure is formed. The neck
and gap structures in the channel map also appear in the near and far sides
that the planet is located. These structures are produced by the density
perturbation, rather than velocity perturbation, as discussed in Appendix C.
The appearance of the substructures in the rare isotopologues is one of
distinguishing features between the coplanar massive planet and the small
planet with an inclined orbit.

Note that one may see the inner secondary gap in the case of $M_{\rm
p}/M_{\ast}=2\times 10^{-3}$ (the bottom line of Figure 5). However, it
locates close to the inner boundary of the 3D simulation. Hence it could be
just the artifact of the simulation. A simulation with the inner boundary at
smaller radii is needed to investigate this issue.

###  3.4 Case with a larger disk scale height ($H_{0}=0.1$)

Figure 6:  The same as in Figure 5, but for the case with $M_{\rm
p}/M_{\ast}=5\times 10^{-4}$, $H_{0}=0.1$ and various $i_{\rm orb}$ (Run
7–10).

We carried out simulations with higher $H_{0}$ (=0.1) and various orbital
inclination (Runs 7–10) to determine the dependence on the disk aspect ratio.
The results are shown in Figure 6. There is no visible gap in the dust
continuum in all the cases shown in the figure. However, in the channel maps,
a significant anomaly from the Keplerian disk is observed when $i_{\rm
orb}=14^{\circ}$ and $21^{\circ}$. In the case with $H_{0}=0.1$, the emission
layer lies at a higher altitude compared to the case of $H_{0}=0.07$ (Figure
1). Considering this, the deviation from the Keplerian rotation in the channel
map of 12CO can be observed with a larger orbital inclination. Similar to the
case with $H_{0}=0.07$, the deviation becomes weaker and vanishes for rarer
isotopologues.

Considering the results with $H_{0}=0.07$ (Figure 5) and $H_{0}=0.1$ (Figure
6), we may find conditions to see the anomaly originated from the planet in
the channel map: the planet can approach the emission layer (i.e., $z_{\rm
em}\sim\sin\left(i_{\rm orb}\right)R_{0}$, where $z_{\rm em}$ is the height of
the emission layer). Since the minimum mass to see the anomaly is $M_{\rm
p}/M_{\ast}=3\times 10^{-4}$ when $H_{0}=0.07$ and $M_{\rm p}/M_{\ast}=5\times
10^{-4}$ when $H_{0}=0.1$, the minimum mass may be given by $M_{\rm
p}/M_{\ast}\sim 0.5H_{0}^{3}$.

###  3.5 Variation of kinematic features in the orbital period

Figure 7:  The variation of the channel map of $\Delta v=1.0$ km/s of 12CO
emission in Run 3 by the meridian location of the planet. In the cases labeled
by $z_{\rm p}=-z_{\max}$ and $z_{\max}$, the planet is located at the highest
altitude of its orbit in $z<0$ or $z>0$ planes (its front/back layer in the
case shown in the figure), where $z_{\rm max}/R_{0}=0.24$ when $i_{\rm
orb}=14^{\circ}$. The labels $z_{\rm p}\rightarrow-0$ and $z_{\rm
p}\rightarrow+0$ indicate the cases that the planet is located at the midplane
($z_{p}=0$) from the $z<0$ and $z>0$ layers.

Figure 8:  Distribution of $v_{r}$ averaged in the range of $\phi=\phi_{\rm
p}\pm R_{H}/R_{0}$ in the case shown in Figure 7. The black cross represents
the location of the planet, and the horizontal line represents the emission
layer of 12CO.

The velocity distribution varies on the meridian location of the planet;
hence, the kinetic features in the channel map also change in the orbital
period of the planet. Figures 7 and 8 show the channel map of 12CO for $\Delta
v=1$ km/s and the distribution of $v_{r}$ averaged in the range of
$\phi=\phi_{\rm p}\pm R_{H}/R_{0}$, respectively. In the left panels of the
figures (the cases labeled by $z_{\rm p}=\pm z_{\max}$), the planet is located
at the highest altitude of its orbit in the front ($z<0$)/back ($z>0$ layers).
Moreover, in the right panels of the figures, the planet is located at the
midplane, but it comes from the front/back layer in the upper/lower panel.
Significant deviations from the Keplerian rotation are observed when the
planet is located in the front layer (in the top panels in Figures 7 and 8).
However, the appearances of the deviation are different depending on the
location of the planet because the velocity distribution is different (Figure
8). Meanwhile, the deviation is much weaker when the planet is located in the
back layer because the emission from the back layer that is disturbed by the
planet is hidden by the front layer. Note that the velocity maps of the left
top and left bottom panels in Figure 8 are not symmetric against the midplane,
but this asymmetry comes from that the plotting timings are slightly off. The
figure that the plotting timing is adjusted is shown in Appendix D .

###  3.6 Comparison with the case with the coplanar planet

Figure 9:  Series of the channel maps of 12CO J=2–1 emission. The upper panels
show the channel maps in the case of the inclined planet with $M_{\rm
p}/M_{\ast}=5\times 10^{-4}$ (Run 3), and the lower panels show these in the
case of the coplanar massive planet with $M_{\rm p}/M_{\ast}=2\times 10^{-3}$
(Run 11).

Kinematic planetary features can appear in a series of channel maps, although
we focused on that for $\Delta v=1$ km/s in aforementioned subsections. In
this section, the series of the channel maps for CO isotopologues are shown,
comparing the cases with the inclined and coplanar planets. Figure 9 shows a
series of channel maps of 12CO emissions in the case of the inclined planet
(Run 3) and the massive coplanar planet (Run 11). In the case of the massive
coplanar planet (lower panels), there are some small deviations from the
Keplerian rotation in all the channels compared with those in the case of the
inclined planet. However, prominent deviations similar with that in the
$\Delta v=1$ km/s appear at the similar location in the same channels, such as
$\Delta v=0$, 0.5, and 1.5 km/s, in both cases.

Figure 10:  The same of the upper panels of Figure 9, but in the case of
$M_{\rm p}/M_{\ast}=3\times 10^{-4}$ (Run 5).

As an intriguing feature only appearing in the case of the inclined planet,
the tail from the planet location in the $\Delta v=-0.5$ km/s channel can be
found. However, this tail can be observed only when the planet is located near
the midplane. Moreover, the tail is much fainter with the smaller mass of the
planet even if the planet locates around the midplane (Figure 10). The series
of the channel maps in the case of $M_{\rm p}/M_{\ast}=3\times 10^{-4}$ (Run
5) is presented in Figure 10. Except for the tail, the structures are similar
to these in the case of $M_{\rm p}/M_{\ast}=5\times 10^{-4}$, which is also
highlighted in Section 3.3. Nevertheless, the tail could be evidence of the
inclined planet and constrained the mass of the planet, if it is observed.
However, this tail originates from the gas very close to the planet, with a
size comparable to the circumplanetary disk, but our simulation does not
resolve the circumplanetary disk (its size is thought to be on the order of
$0.1R_{H}$). Therefore, this tail structure can be the artifact caused by the
unresolved circumplanetary disk.

Figure 11:  The same as Figure 9, but for the C18O J=2–1 emissions.

Figure 11 shows the series of channel maps of C18O J=2–1 emissions, which is
the same cases shown in Figure 9. In the case of the inclined planet shown in
the upper panels of the figure, there is almost no deviation from the pattern
of the Keplerian rotation, as discussed in Section 3.3 and shown in Figure 5.
In contrast, the prominent gap is observed in all the channels in the case of
the massive coplanar planet, as shown in the lower panels of Figure 11. The
difference in the channel map of C18O is much distinct as compared with the
difference in that of 12CO. Hence, the observations of C18O can be decisive to
determine the orbital inclination of the planet.

Figure 12:  Moment 1 maps for 12CO (upper panels) and C18O (lower panels) with
$M_{p}/M_{\ast}=5\times 10^{-4}$ (Run 3, left panels) and
$M_{p}/M_{\ast}=2\times 10^{-3}$ (Run 11, right panels). The contours denotes
the lines of the moment 1 being equal to $0$ km/s, $\pm 0.5$ km/s, $\pm 1.0$
km/s, and $\pm 1.5$ km/s .

The trends appearing in the channel maps shown above also appear in the moment
1 maps shown in Figure 12. As similar to Figure 9, one finds the deviation
from the Keplerian rotation in the 12CO maps (the upper panels of Figure 12).
On the other hands, in the C18O maps, only a weak deviation can be found in
the case of Run 3, while the case of Run 11 shows the deviation as strong as
that shown in the 12CO map. This is because the inclined planet does not
perturb around the midplane and the perturbation is weaker as closer to the
midplane, while the massive co-planer planet affects both the midplane and
upper layers, as shown above.

##  4 Discussion

###  4.1 Implication to observations

When the planet has an inclined orbit, the kinematic features in the channel
map has several differences from those induced by the coplanar planet, as
presented in Section 3. First, the kinematic features induced by the planet
with an inclined orbit is significant in common isotopologues, such as 12CO ,
because the planet perturbs stronger in a higher altitude (Figure 3). For the
rarer isotopologues, such as 13CO and C18O, only much weaker features are
observed (Figure 5). In contrast, in the case of the coplanar planet, a gap
feature can be seen and is clear for rarer isotopologues, because they trace a
closer location to the midplane with more decreasing gas density. The
dependence of kinematic features on isotopologues can distinguish the cases
with the planet coplanar and inclined orbits. The gap is narrower and
shallower when these are formed by the inclined planet, as compared with those
induced by the coplanar planet.

Moreover, we pointed out that, in the outer region of the disk with the disk
aspect ratio, the inclined planet can make significant kinematic planetary
features in the channel map of 12CO with no visible gap in dust emission
(Figure 6). If only the kinematic planetary signature is observed with no dust
gap, it may indicate a relatively small planet in the outer region.

The appearance of kinematic planetary features varies relevant to the vertical
location of the planet (Figure 7). Kinematic planetary features can be changed
because of the azimuthal rotation of the planet. If it is formed by the
inclined planet, an additional time variation due to the variation of the
vertical location of the planet should be observed. Moreover, the planet in
the front layer from the observer is likely to be observed because the planet
in the back layer produces the weak features in the channel map (Figure 7).
That is, half of the planets are hardly detectable. This bias should be
considered when the planetary population from the observations of the
kinematic planetary feature is discussed.

We should note that in this paper we have fixed the CO/H2 ratio, the dust/gas
ratio, and the H2 gas density itself, and we have also ignored the freezing of
CO at very low temperatures, but these quantities determine the height of the
emitting layer and thus also affect the appearance of kinetic features. In a
less massive disk than our model in this paper, for instance, the emission
layer of 12CO lies closer to the midplane, and the neck structure can be seen
in 12CO channel map in this case, and planetary features of CO emission are
fainter when the freezing of CO is significant. Moreover, the degree of the
dust-settling and size distribution of the dust grains can affect the disk
temperature (e.g., Zhang et al. (2021b)), which also affects the height of the
emission layer. For direct comparison with the observation and numerical
simulations, we should consider non-isothermal equation of state because it
changes the shape of the match cone formed by the planet. With such the
uncertainly, however, comparing the channel maps of different isotopologues,
we may distinguish the features induced by the inclined planet and co-planer
planet. We discuss the case of HD 163296 in detail as an example, below.

Pinte et al. (2020) has reported nine disks with significant anomaly from the
Keplerian rotation in the channel map of 12CO emission, among the DSHARP
sample. According to the results of hydrodynamic simulations with the coplanar
planet, they suggested that a massive planet with $\gtrsim 1M_{J}$ is required
to form the observed anomalies from the Keplerian rotation, which is referred
as velocity kinks. However, they also reported that the planet mass estimated
from the gap shape (i.e., width and depth) is $\sim 0.1$ times smaller than
that suggested by the velocity kink. This discrepancy of the planet masses
estimated from the gap and velocity kink may be resolved when the inclined
planet is considered, as discussed in Section 3. If the observations of other
isotopologues or other transitions are provided, we could determine which is
more suitable to the observation between the coplanar and the inclined
planets.

Figure 13:  A part of the channel map of 12CO (left), 13CO (middle), and C18O
(right) for the disk of HD 163296 is provided by MAPS (Öberg et al., 2021).
The upper and lower panels show the channel map with $\Delta v=1$ km/s and
$\Delta v=-2.4$ km/s, respectively. The white circles indicate the locations
of the velocity kink reported by Pinte et al. (2020). The white arrows in the
lower middle and right panels represent the neck and gap structures in the
opposite side of the rotation pattern against the planet location. The system
velocity of HD 163296 is 5.8 km/s (Zhang et al., 2021a).

Figure 13 shows the part of the channel maps for 12CO , 13CO, and C18O, for
the disks of HD 163296, which is provided by MAPS (Öberg et al., 2021). In the
disk of HD 163296, Pinte et al. (2020) found two localized velocity kinks: one
is seen in the $\Delta v=1$ km/s map at 2.2 arcsec (220 au) (called the outer
kink), and the other is seen in the $\Delta v=-2.4$ km/s map at 0.67 arcsec
(67 au) (called the inner kink). For the outer kink, the significant kink in
the 12CO map is found, while there is no significant kink feature in the same
location of 13CO kink in the maps of 13CO and C18O. As discussed before, such
a dependence on the isotopologues is consistent with that induced by the
inclined planet, rather than the coplanar one. In contrast, for the inner
kink, the neck-like feature and the gap can be observed in the 13CO and C18O
maps. The neck and gap structures can also be observed in the opposite side of
the rotation pattern (represented by the arrows in the figure). These
structures are consistent with that induced by the coplanar planet, as shown
in the bottom panels of Figure 5. We pointed out that the point-like emission
within the white circle in the C18O map with $\Delta v=-2.4$ km/s, as shown in
the lower right panel of Figure 13). The location of the point-like emission
may be consistent to the location of the planet estimated by Pinte et al.
(2020). If the planet at 220 au has the inclined orbit, the density waves
induced by the planet is weaker than that estimated with the co-planer planet,
which could explain that Speedie & Dong (2022) explored the density waves in
dust continuum but found no clear detection of any spiral structure around 220
au of HD 163296.

###  4.2 Possible mechanisms to maintain the orbital inclination

The orbital inclination is damped due to disk–planet interaction; however, we
did not consider this in our paper. In the linear regime, the orbital
inclination exponentially decays with time and the damping timescale is much
shorter than the orbital migration time (Tanaka & Ward (2004)). Numerically,
Cresswell et al. (2007) and Bitsch & Kley (2011) confirmed that the orbital
inclination exponentially decays for the small planet with a small
inclination, as predicted by the linear theory, except in the large
inclination limit as discussed below. In this case, the inclination of $\sim
10^{\circ}$ can decay within several ten orbits. For a massive planet with a
deep gap, the damping rate can be longer than that estimated by the linear
theory (Bitsch et al., 2013; Zhu, 2019). However, it is not appropriate with
the cases that we treated in this paper because only a shallow gap is opened
by the inclined planet. Hence, if the observed features are formed by the
inclined planet, the mechanism to maintain the inclination is required.

The inclination may be maintained by the balance between damping and
excitation, which are due to the disk-planet interaction and planet– planet
interaction, respectively, when the planets are captured by the mean-motion
resonance (e.g., Sotiriadis et al. (2017)). Almost no visible gap and no
kinetic features are observed when the other planet of the pair has a large
inclination, as shown in Figure 2. Rein (2012) investigated the damping
timescale of the inclination for a highly inclined planet by the dynamical
friction theory. In addition, he showed that the damping because of the
disk–planet interaction is inefficient as the inclination increases. Moreover,
the inward migration is slower with a larger inclination. When the inner
planet has a large inclination, the inclination of the outer planet can be
excited up to $\sim 10^{\circ}$, and the inward migration can also slow down
due to the interaction with the inner planet, which is consistent with the
fact that the kinematic features are frequently observed in the outer region
of the disk. It is difficult to directly detect the inner planet with a large
inclination because it makes no visible gap and kinematic features. However,
it may be revealed by the inclination and eccentricity of the outer planet
that makes the observed kinematic features because the inclination and
eccentricity should be rapidly damped without the inner planet.

##  5 Summary

In this paper, we investigated the kinematic features in the channel maps of
CO isotopologues induced by the inclined planet through 3D hydrodynamic
simulations. We demonstrated that the appearance of the kinematic features
induced by the inclined planet is different from that induced by the coplanar
planet at several points as follows:

  1. 1.

When the planet has an inclination of $10^{\circ}$ – $20^{\circ}$, the smaller
planet ($M_{\rm p}/M_{\ast}=3\times 10^{-4}$) can make the kinematic features
as clear as that induced by the coplanar massive planet ($M_{\rm
p}/M_{\ast}=2\times 10^{-3}$) in the channel map of 12CO (Figure 5).

  2. 2.

In the case of the inclined planet, the kinematic features is fainter for
rarer isotopologues because the velocity perturbation is weaker at the
position closer to the midplane (Figure 3 shows that the emission layer of the
rarer isotopologues lies closer to the midplane (Figure 1). In contrast, in
the case of the coplanar massive planet (the bottom panels of Figure 5), the
neck and gap structures are observed in the channel maps of the rarer
isotopologues.

  3. 3.

The inclined planet can make the significant kinematic planetary feature with
no prominent gap in dust emission in the outer region of the disk with the
large disk aspect ratio (Figure 6).

  4. 4.

The appearance of the kinematic planetary feature induced by the inclined
planet varies with time, which is relevant to the vertical location of the
planet (Figure 7). Hence, the time variation of the feature that is not
related to the azimuthal rotation would be observed if it is induced by the
inclined planet. Moreover, only half the planets are detectable, which is
important to consider planetary population, since the planet in the back layer
from the observer produces much weaker features as compared to the planet in
the front layer.

As discussed in Section 4.1, compared to the case of the coplanar planet, the
fact that the inclined planet can make a shallower gap and prominent kinematic
features may resolve the discrepancy, which is pointed by the previous works
(e.g.,Pinte et al. (2020)), in which the planet mass required to explain the
observed kink is $4$–$10$ times larger than that estimated from the observed
gap shape. Moreover, the channel maps of 13CO and C18O can provide a clue to
distinguish between the cases with the coplanar massive and inclined planets.
From the structures of the channel maps of 12CO, 13CO, and C18O in in the case
of HD 163296, we concluded that the outer velocity kink at $2.2$ arcsec is
consistent with that induced by the inclined planet. In contrast, the inner
velocity kink at $0.67$ arcsec has the neck and gap structures in the maps of
13CO and C18O, respectively, which is consistent with the case with the
coplanar planet.

The orbital inclination of the planet significantly affects the gap shape and
kinematic features in the channel map. Hence, in the future, a more
sophisticated model including the planet inclination is required to confirm
the properties of the embedded planet.

We would like to thank the anonymous referee for his/her careful review and
helpful comments. This work was supported by the JSPS KAKENHI (Grant Nos.
17H01103, 18H05441 and 19K14779). Numerical computations were carried out on
the Cray XC50 at the Center for Computational Astrophysics, National
Astronomical Observatory of Japan. This paper uses the following ALMA data:
ADS/JAO.ALMA#2018.1.01055.L. ALMA is a partnership of ESO (representing its
member states), NSF (U.S.), and NINS (Japan), together with NRC (Canada), NSC
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The joint ALMA observatory is operated by ESO, auI/NRAO,
and NAOJ.

## Software

ATHENA++ (Stone et al. (2020), https://www.athena-astro.app), RADMC3D version
2.0 (Dullemond & Dominik, 2005) and radmc3dPy python package
(https://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d), Matplotlib
(Hunter (2007), http://matplotlib.org), NumPy (Harris et al. (2020),
http://www.numpy.org).

##  Appendix A Channel maps in the vertically isothermal disk

Figure 14:  The series of the channel maps of 12CO (upper panels) and C18Oin
the same case are shown in the upper panels of Figure 9, but the temperature
distribution given by Equation (8) is adopted.

Figure 14 shows the channel map in the same case shown in the upper panel of
Figure 9. However, the temperature calculated from Equation (8), vertically
isothermal, is adopted. Comparing the case that the temperature calculated by
the thermal-Monte Carlo calculations (the upper panels of Figures 9 and 11),
the intensities becomes much fainter. The temperature at an altitude is the
same as the midplane temperature in the case shown in Figure 14, while the
temperature at an altitude is higher than the midplane temperature when it is
given by the thermal Monte Carlo calculation. The difference of the intensity
is originated from the difference of the temperature around the emission
layer. However, except for the magnitude of the intensity, the channel maps is
considerably similar to each other. Although our simulations contain a
temperature contradiction adopted in hydrodynamic simulations and radiative
transfer simulations, this contradiction is not expected to significantly
affect the main results of this paper.

##  Appendix B Disk structure in the case of the massive coplanar planet

Figure 15:  Distributions of densities and velocities in the case of Run 11.
The dashed, dotted, and solid horizontal lines represent the emission layers
of 12CO, 13CO, and C18O , respectively, in the unperturbed disk (shown in
Figure 1). The actual heights of the emission layers are affected by gap
formation, and these are shown in Figure 16.

Figure 16:  The same as Figure reffig:tau1surf, but for the case of Run 11.

The distributions of densities and velocities in the case of Run 11 is shown
in Figure 15, and the heights of the emission layers of CO isotopologues are
shown in Figure 16.

##  Appendix C Velocity perturbation and the kinematic planetary feature

Figure 17:  The channel map of 12CO in $\Delta v=1$ km/s channel in the case
of Run 3. In the panel labeled as ’full’, all physical quantities are fed from
the result of hydrodynamic simulation. In the panel labeled as ’all velocities
are initial’ (upper left), all velocities that are given by the initial
profile (Keplerian rotation) are adopted instead of those produced by the
hydrodynamic simulation, and the lower right and left panels, $v_{r}$ and
$v_{\phi}$, respectively, are set as the initial profiles.

In this section, we discuss which components of the velocity mainly contribute
to the appearance of the kinematic planetary features. Considering this, we
demonstrated additional radiative transfer simulations in the following three
setups: first, the initial profiles of $v_{r}$, $v_{\theta}$, and $v_{\phi}$
were adopted instead of calculating through hydrodynamic simulations. Hence,
in this case, only the density was fed from the hydrodynamic simulation. In
the other two setups, we adopted the initial $v_{r}$ and $v_{\phi}$,
respectively. The channel maps with these setups are presented in Figure 17,
when $M_{\rm p}/M_{\ast}=5\times 10^{-4}$ and $i_{\rm orb}=14^{\circ}$ (Run
3). Based on the figure, when the initial $v_{\phi}$ is applied, the deviation
from the Keplerian pattern is observed to be as clear as that in the case
where all velocities are fed from the simulation, whereas there is no
deviation in the other two cases. Therefore, the perturbation of $v_{r}$
mostly makes the kinematic planetary feature in this case.

Figure 18:  Similar as Figure 18, but for the case of Run 11 (the case of the
massive coplanar planet).

Figure 18 shows the comparison similar to that presented in Figure 17 in the
case of the massive coplanar planet (Run 11). The kink feature around the
planet is originated from the perturbation of the $v_{r}$, which is similar to
Figure 17. However, even in the case where all velocities are set to be the
initial profiles, deviations from the Keplerian pattern can be found, as
opposed to the case presented in Figure 17. Therefore, the perturbation of
$\rho$ and $v_{\phi}$ marginally contribute to make the kinematic planetary
features and the perturbation of $v_{r}$.

Figure 19:  The same as Figure 18, but for C18O emissions.

Figure 19 is the same as in Figure 18, but for C18O emission. In this case,
even the case where all velocities are set as initial profiles, a similar
image is produced. Hence, in this case, the deviation from the Keplerian
pattern is mostly originated from the density perturbation.

##  Appendix D Time correction of Figure 8

Figure 20:  The same as Figure 8, but the plotting timing of the left top
panel is adjusted to match that of the left bottom panel.

The plotting time of the left top panel in Figure 8 is later in
$\pi/\Omega_{K}/6$ that that corresponding to the plotting time of the left
bottom panel, and it is fixed in Figure 20.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41
  * Arzamasskiy et al. (2018) Arzamasskiy, L., Zhu, Z., & Stone, J. M. 2018, MNRAS, 475, 3201
  * Bae et al. (2017) Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201
  * Bae et al. (2022) Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20
  * Barraza-Alfaro et al. (2021) Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45
  * Bitsch et al. (2013) Bitsch, B., Crida, A., Libert, A. S., & Lega, E. 2013, A&A, 555, A124
  * Bitsch & Kley (2011) Bitsch, B., & Kley, W. 2011, A&A, 530, A41
  * Calcino et al. (2022) Calcino, J., Hilder, T., Price, D. J., et al. 2022, ApJ, 929, L25
  * Casassus & Pérez (2019) Casassus, S., & Pérez, S. 2019, ApJ, 883, L41
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698
  * Cresswell et al. (2007) Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329
  * Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, icarus, 181, 587
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127
  * Duffell & MacFadyen (2013) Duffell, P. C., & MacFadyen, A. I. 2013, ApJ, 769, 41
  * Dullemond & Dominik (2005) Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool
  * Flores-Rivera et al. (2020) Flores-Rivera, L., Flock, M., & Nakatani, R. 2020, A&A, 644, A50
  * Fung et al. (2014) Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88
  * Goldreich & Tremaine (1980) Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Izquierdo et al. (2022) Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2
  * Kanagawa et al. (2015) Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994
  * Kley (1999) Kley, W. 1999, MNRAS, 303, 696
  * Kloster & Flock (2019) Kloster, D. L., & Flock, M. 2019, MNRAS, 487, 5372
  * Lin & Papaloizou (1979) Lin, D. N. C., & Papaloizou, J. 1979, MNRAS, 186, 799
  * Lin & Papaloizou (1986) —. 1986, ApJ, 307, 395
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
  * Marzari & Nelson (2009) Marzari, F., & Nelson, A. F. 2009, ApJ, 705, 1575
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1
  * Okuzumi et al. (2016) Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82
  * Pérez et al. (2018) Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018, MNRAS, 480, L12
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJ, 811, L5
  * Pinte et al. (2022) Pinte, C., Teague, R., Flaherty, K., et al. 2022, arXiv e-prints, arXiv:2203.09528
  * 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
  * Pinte et al. (2020) Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9
  * Rabago & Zhu (2021) Rabago, I., & Zhu, Z. 2021, MNRAS, 502, 5325
  * Rein (2012) Rein, H. 2012, MNRAS, 422, 3611
  * Rich et al. (2021) Rich, E. A., Teague, R., Monnier, J. D., et al. 2021, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 312.01
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Sotiriadis et al. (2017) Sotiriadis, S., Libert, A.-S., Bitsch, B., & Crida, A. 2017, A&A, 598, A70
  * Speedie & Dong (2022) Speedie, J., & Dong, R. 2022, ApJ, 940, L43
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4
  * Takahashi & Inutsuka (2014) Takahashi, S. Z., & Inutsuka, S.-i. 2014, ApJ, 794, 55
  * Takahashi & Inutsuka (2016) —. 2016, AJ, 152, 184
  * Takahashi & Muto (2018) Takahashi, S. Z., & Muto, T. 2018, ApJ, 865, 102
  * Tanaka & Ward (2004) Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12
  * Teague et al. (2019) Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019, ApJ, 884, L56
  * Tsukagoshi et al. (2016) Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJ, 829, L35
  * van der Marel et al. (2019) van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112
  * van der Tak et al. (2020) van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15
  * Winn & Fabrycky (2015) Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409
  * Xiang-Gruess & Papaloizou (2013) Xiang-Gruess, M., & Papaloizou, J. C. B. 2013, MNRAS, 431, 1320
  * Zhang et al. (2015) Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7
  * Zhang et al. (2021a) Zhang, K., Booth, A. S., Law, C. J., et al. 2021a, ApJS, 257, 5
  * Zhang et al. (2021b) Zhang, S., Hu, X., Zhu, Z., & Bae, J. 2021b, ApJ, 923, 70
  * Zhu (2019) Zhu, Z. 2019, MNRAS, 483, 4221