# Mapping the Planetary Wake in HD 163296 with Kinematics

Josh Calcino, Thomas Hilder, Daniel J. Price, Christophe Pinte, Francesco Bollati, Giuseppe Lodato, Brodie J. Norfolk

► To cite this version:

Josh Calcino, Thomas Hilder, Daniel J. Price, Christophe Pinte, Francesco Bollati, et al.. Mapping the Planetary Wake in HD 163296 with Kinematics. The Astrophysical journal letters, 2022, 929, 10.3847/2041-8213/ac64a7. insu-03705358

HAL Id: insu-03705358

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

Submitted on 27 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# Mapping the Planetary Wake in HD 163296 with Kinematics

Josh Calcino<fnref n="1" /> ID, Thomas Hilder<fnref n="2" /> ID, Daniel J. Price<fnref n="2" /> ID, Christophe Pinte<fnref n="2" />,<fnref n="3" /> ID, Francesco Bollati<fnref n="4" />, Giuseppe Lodato<fnref n="5" /> ID, and Brodie J. Norfolk<fnref n="6" /> ID

## Abstract

We map the planetary wake associated with the embedded protoplanet creating the CO kink in the disk of HD 163296. We show that the wake can be traced by a series of correlated perturbations in the peak velocity map. The sign change of the perturbations across the disk's major axis confirms that the wake induces predominantly radial motion, as predicted by models of planet-disk interaction. These results provide the first direct confirmation of planet wakes generated by Lindblad resonances. Mapping the wake provides a constraint on the disk aspect ratio, which is required to measure the mass of the planet.

*Unified Astronomy Thesaurus concepts: Planet formation (1241); Protoplanetary disks (1300); Hydrodynamical simulations (767); Radiative transfer (1335)*

## 1. Introduction

The disk around HD 163296 displays the kinematic signatures associated with an embedded planet (Pinte et al. 2018a), which was first hypothesized by Grady et al. (2000). Higher-spatial-resolution observations have revealed additional substructures in millimeter-continuum and CO line observations from the Disk Substructures At High Angular Resolution Project (DSHARP) (Andrews et al. 2018; Huang et al. 2018; Pinte et al. 2020; Izquierdo et al. 2022) and more recently the Molecules with Atacama Large Millimeter/submillimeter Array (ALMA) at Planet forming Scales (MAPS) (Oberg et al. 2021; Teague et al. 2021) large programs using ALMA. Teague et al. (2021) pointed out spiral features both in channel maps and velocity residuals. Velocity kinks arise in the CO channel maps due to perturbations in the motion of gas, which modifies the projected velocity. We define a “velocity kink” as a localized distortion of the isovelocity curve, caused by deviations from Keplerian motion, which shift the line emission into an adjacent channel. They can be induced by, for example, spiral arms, gravitating bodies, or gravitational instability in the disk. In particular, velocity kinks detected in the CO emission were associated with an embedded planet (Pinte et al. 2018a).

In Bollati et al. (2021, hereafter B21), we recently developed a semianalytic theory of velocity kinks caused by planet-disk interaction, based on the mathematical theory of planet-induced density waves launched at Lindblad resonances developed by Goldreich & Tremaine (1979), Goldreich & Tremaine (1980), Goodman & Rafikov (2001), Ogilvie & Lubow (2002), and Rafikov (2002). Applying this model to the kinematic detection of a planet in HD 163296 made by Pinte et al. (2018a) allows one to interpret velocity kinks as occurring whenever the spiral wake from the planet intersects a velocity channel in the data.

More puzzling was our finding in B21 that such kinks should therefore extend throughout the disk, including far from the planet (see Figure 6 of B21). Hence, nonlocalized “secondary” kinks away from the planet should be detectable. In B21 we tried to explain away the apparent conflict between the nondetection of these secondary kinks in the data—perhaps being washed out by the finite beam, or suppressed due to viscous or other forms of damping occurring in the disk. With the MAPS data, the former hypothesis is no longer tenable.

In this Letter, we demonstrate that not only are secondary kinks observed in HD 163296, but they can also be used to map the wake generated by the embedded planet found by Pinte et al. (2018a) and Teague et al. (2018) over a large fraction of the disk. In particular, we show that the spiral wake generated by the planet can be easily traced in a peak velocity map given sufficient spatial and spectral resolutions. By determining the shape of the planet wake we constrain the aspect ratio of the disk, which is needed in order to measure the planet mass from observed velocity kinks (B21).

## 2. Methods

### 2.1. Observations

We used ¹²CO (2–1) robust = 0.5 line emission observations of HD 163296 from the MAPS large program (2018.1.01055.L, Czekala et al. 2021; Oberg et al. 2021).<fnref n="7" /> The observations have a JvM correction (Jorsater & van Moorsel 1995), which accounts for a mismatch between the CLEAN and dirty beam sizes when a CLEAN model image is combined with the residuals. The observations have a beam size of approximately 0.′′14 × 0.′′11 with a position angle of 104°, a spectral resolution of 92 m s⁻<fnref n="1" />, and a spacing of 200 m s⁻<fnref n="1" />. These are the same observations used by Teague et al. (2021) in their kinematic analysis of HD 163296. We assumed a line-of-sight velocity of $v_{\text{los}} = 5.76$ km s⁻<fnref n="1" />, as found by

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA; jcalcino@lanl.gov</footnote>

<footnote n="2">² School of Physics and Astronomy, Monash University, Clayton, Vic 3800, Australia</footnote>

<footnote n="3">³ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4">⁴ Dipartimento di Scienza e Alta Tecnologia, Università degli Studi dell'Insubria, Via Valleggio 11, I-22100, Como, Italy</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16 I-20133 Milano, Italy</footnote>

<footnote n="6">⁶ Centre for Astrophysics and Supercomputing (CAS), Swinburne University of Technology, Hawthorn, Victoria 3122, Australia</footnote>

<footnote>Received 2021 November 10; revised 2022 March 29; accepted 2022 April 6; published 2022 April 22</footnote>

<footnote n="7">⁷ The data are available for download at http://alma-maps.info/.</footnote>
Teague et al. (2021). The emitting layer of CO was estimated by Law et al. (2021) using the method developed by Pinte et al. (2018b). The resulting surface layer was fitted by a flared disk structure with an exponential taper, with the following functional form:

$$z(r) = z_0 \left(\frac{r}{1^{\prime\prime}}\right)^{\phi} \exp\left(-\left[\frac{r}{r_{\text{taper}}}\right]^{\psi}\right). \quad (1)$$

We adopted the best-fit values of HD 163296 from Law et al. (2021) of $z_0 = 0.388$, $\phi = 1.851$, $r_{\text{taper}} = 2^{\prime\prime}362$, and $\psi = 1.182$ when computing the CO-emitting layer.

## 2.2. Hydrodynamical Simulations and Radiative Transfer

We used the smoothed particle hydrodynamics (SPH) code PHANTOM (Price et al. 2018) to simulate the interaction of a 3 $M_J$ planet with a gas-only disk in 3D. We employed $2 \times 10^6$ SPH particles to model the gas disk with mass $10^{-2} M_\odot$ orbiting a 1.9 $M_\odot$ central star. Both the planet and central star were modeled using sink particles (Bate et al. 1995). The central star was given an accretion radius of 5 au, while the planet has a radius of 5.56 au, or 0.25 of the Hill radius. The SPH gas particles were initialized with a pressure-corrected Keplerian velocity, assuming a power-law distribution with an exponential taper given by

$$\Sigma(r) \propto (r/r_{\text{p,i}})^{-p} \times \exp[-(r/r_c)^{2-p}], \quad (2)$$

where $r_{\text{p,i}} = 280$ au is the initial location of the planet, $p = 1$, and $r_c = 100$ au. The temperature profile was assumed to be vertically isothermal with the sound speed $c_s \propto (r/r_{\text{p,i}})^{-q}$, where we set $q = 0.35$. The preceding parameters were selected following the model from Pinte et al. (2018a), except for the initial planet location, which was adjusted to allow for migration. We assumed a disk aspect ratio $(h/r)_{\text{p,i}} = 0.1$ following de Gregorio-Monsalvo et al. (2013) and Pinte et al. (2018a). The disk was allowed to evolve for 50 orbits of the planetary companion, at which point it had migrated to a radius of $r_p = 256$ au from the central star where $h_p/r_p \sim 0.09$ and has a mass of $\sim 4 M_J$. This difference in mass compared to the analytics is not significant, being smaller than the measurement error from kinematics at present (e.g., Pinte et al. 2018a).

The output of our SPH simulation was used to generate a Voronoi mesh for input into the Monte Carlo radiative transfer code MCFOST (Pinte et al. 2006, 2009). Because we did not include dust in our simulation, we assumed that the dust followed the gas particles with a power-law grain size distribution $dn/ds \propto s^{-3.5}$ for $0.03 \ \mu m \le s \le 1 \ mm$ with a gas-to-dust ratio of 100. We assumed spherical and homogeneous grains composed of astronomical silicate (Weingartner & Draine 2001). The central star is set with an effective temperature $T_{\text{eff}} = 9250$ K and a radius $R_* = 1.6 R_\odot$, giving a luminosity of $\sim 17 L_\odot$ ( Setterholm et al. 2018).

The temperature profile and specific intensities were generated using $10^8$ photon packets. Images were produced by ray-tracing the computed source function. Our CO isotopologue observations were generated assuming $T_{\text{gas}} = T_{\text{dust}}$ and all molecules were assumed to be in Local Thermodynamical Equilibrium (LTE). We assumed an initially constant CO abundance, with an abundance ratio of $^{12}\text{CO}/\text{H}_2 = 5 \times 10^{-5}$. This abundance was then altered due

to CO freeze-out ($T = 20$ K, depletion factor of $10^{-4}$), photodesorption, and photodissociation, following Appendix B of Pinte et al. (2018b).

Channels for the CO emission are created at a separation of $20 \ \text{m s}^{-1}$. We linearly interpolate over five channels and average the interpolated data after weighting by a Hann window function, producing a channel width of $100 \ \text{m s}^{-1}$ with a spacing of $200 \ \text{m s}^{-1}$. The images were then smoothed with a Gaussian beam with size $0^{\prime\prime}15 \times 0^{\prime\prime}15$.

## 2.3. Spiral-wake Mapping

In a Keplerian power-law disk, $c_s(r) = c_{s,p}(r/r_p)^{-q}$, $\Sigma = \Sigma_p(r/r_p)^{-p}$ and $\Omega_K = (GM_*/r^3)^{1/2}$, where $r_p$ is the planet location, $c_s$ is the sound speed, and $\Sigma$ is the column density, the planetary wake generated by a planet lies on the curve given by Rafikov (2002),

$$\varphi_{\text{wake}}(r) = \varphi_{\text{p}} + \text{sign}(r - r_{\text{p}}) \left[ \frac{h_{\text{p}}}{r_{\text{p}}} \right]^{-1} \left[ \frac{(r/r_{\text{p}})^{q-1/2}}{q-1/2} - \frac{(r/r_{\text{p}})^{q+1}}{q+1} - \frac{3}{(2q-1)(q+1)} \right], \quad (3)$$

where $(h_p/r_p)$ is the scale height at the planet location. We initially assumed $q = 0.35$ and $h_p/r_p = 0.09$ as in the simulations but found the outer wake matched better to the prominent kink in the outer disk with a slightly lower $h_p/r_p$ of 0.08.

We followed the prescription outlined in B21 to compute the velocity perturbations created by the planet wake in the disk midplane using nonlinear evolution theory, with the same disk and planet parameters as above. We generated these perturbations on a $500 \times 500$ Cartesian grid and added them to an unperturbed background given by Takeuchi & Lin (2002):

$$v_{\phi}(r) = \left( \frac{GM_{*}}{r} \right)^{\frac{1}{2}} \left[ 1 - (p+q) \left( \frac{h_{p}}{r_{p}} \right)^{2} \left( \frac{r}{r_{p}} \right)^{1-2q} \right]^{\frac{1}{2}}, \quad (4)$$

where the deviation from Keplerian rotation is due to the gas pressure support.

We projected the wake and associated velocity fields (Equation (3)) onto the CO-emitting layer (Equation (1)) assuming no vertical dependence in the velocity perturbations (see Section 4.3) and then onto the plane of the sky to derive a line-of-sight velocity map of the disk upper surface (using a disk inclination of $i = 46^\circ$ and position angle $PA = 313^\circ$). In this projection, we implicitly assume identical radial temperature profiles in both the midplane and the emitting layer. In reality, the temperature is vertically stratified, which may result in spirals with a larger pitch angle (see Section 4.3).

## 3. Results

Figure 1 displays a selection of CO channel maps of the disk around HD 163296 (subplot A) along with analogous channels from our SPH model (subplot B). We have annotated several features of interest. First, in the top-left panel, one may observe the velocity kink associated with an embedded planet, first reported by Pinte et al. (2018a). In the third panel of the top row, we highlight a secondary kink on the far side of the disk, which is at a larger radial separation from the central star than the embedded protoplanet.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Secondary kinks seen in channel maps of $^{12}$CO (2–1) line emission in HD 163296 (Oberg et al. 2021) (subplot A), along with the same channels from our SPH simulation (subplot B). Labels indicate the location of the embedded planet, as well as secondary kinks caused by the spiral wake, denoted by the cyan, dashed line (Equation (3)); projected to the CO-emitting layer). The predicted location of the wake explains the main spiral features seen in the channel maps. The magenta line in the second panel on the bottom of each subplot is the wake projected onto the lower surface. The opaque crosses mark the locations of additional kinks that do not appear associated with the planetary wake.</figcaption></figure>

We find the same structure in the SPH model (third panel from left in the top half of subplot B). This is the outer wake of the planet. Interestingly, this kink was predicted by the original SPH simulations by Pinte et al. (2018a) (see their Figure 5) and in retrospect can be seen in the original data used by Pinte et al. (2018a) to detect the embedded protoplanet.

In subplot A of Figure 1, we have labeled a faint kink seen on the lower surface of the disk (bottom row, second from left). An analogous feature is also seen in our SPH model. This is the outer wake of the planet seen on the lower surface of the disk. We plot the projected location of the outer planet wake on the lower surface to demonstrate this (magenta dashed line). Our SPH model suggests we should also see the spiral kink on the

upper surface of the disk. Although there is a hint of this in the observations, the contrast is lower in this location with respect to the midplane and the lower surface.

Opaque crosses in Figure 1 denote several additional kinks visible in the channel maps. Not all of these kinks are associated with the planet wake generated by the embedded protoplanet (see Section 4.1). The height of the CO-emitting layer is lower in our model than in the observations. Increasing the disk aspect ratio would correct this; however, the resulting planet wake would not be tightly wound enough as in the observations. The most likely reason for this discrepancy is a slight difference between the temperature profile assumed in the SPH model (which sets the wake propagation) and the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Analytic prediction of the velocity perturbation generated by the planetary wake of a 3 $M_J$ planet located at 250 au. The left and right panels show radial and azimuthal velocity perturbations, respectively. Note that positive radial velocities indicate motion away from the central star, and positive azimuthal velocities indicate super-Keplerian flow.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Channel map in $^{12}$CO (2–1) line emission at $\Delta v = 0.5$ km s$^{-1}$ from systemic (left panel). We define positive velocity as moving away from the observer (redshifted) while negative velocity is toward the observer (blueshifted). Labels highlight redshifted and blueshifted material caused by the crossing of the planetary wake. The sign flip over the semimajor axis indicates that the dominant velocity perturbation is radial, as expected for planet-induced spiral density waves (right panel).</figcaption></figure>

temperature profile determined self-consistently in the radiative transfer by MCFOST (which sets the CO-emitting layer). Proper treatment of the disk thermal structure in the SPH simulation is required to better understand the propagation of the wake and the disk scale height. We discuss the implications of this more in Section 4.3.

Figure 2 shows the velocity perturbations predicted for a 3 $M_J$ planet in the B21 analytic model. The left and middle panels show the velocity contribution in the radial and azimuthal directions, respectively. Comparing the left and middle panels, we see that near the planet the azimuthal components dominate, as also seen in the supplementary material of Pinte et al. (2019), while far from the planet, the radial motions dominate. That $\Delta v_{\phi} \ll \Delta v_r$ in the wake was mathematically proved by Rafikov (2002, Equations (A4) and (A8), respectively) so this is no surprise. The azimuthal contribution is only important within

~10 au of the planet. Figure 2 (right panel) shows the projection of the velocity perturbations onto the surface of the CO-emitting layer in HD 163296. Because the perturbations are mostly radial, their projected amplitude decreases substantially along the semimajor axis of the disk where these perturbations are perpendicular to the line of sight.

Figure 3 highlights the perturbations caused by the planet wake crossing the $\Delta v = 0.2$ km s$^{-1}$ channel in the observations. Labels in the figure highlight the sign flip over the semimajor axis, as expected for radial-dominated motion (see above). More subtle but visible nonetheless is the N-wave structure transverse to the wake predicted from the analytics (Rafikov 2002, B21).

Figure 4 compares the peak velocity map of HD 163296 to the semianalytic model (top right) and to our SPH + radiative transfer model (bottom right). Perturbations from Keplerian
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Peak velocity map of the $^{12}$CO (2–1) emission of HD 163296 from MAPS (Oberg et al. 2021) (panel A), the analytic model prediction (panel B), and the prediction from our SPH + radiative transfer model (panel C). The projected wake in the simulation and analytic models is shown by the solid line and can be seen to trace a series of kinks mapping the wake over the whole disk. The same series of kinks is visible in the observations (panel A), revealing the wake from the planet in the data, almost perfectly corresponding to the projected wake from the models.</figcaption></figure>

rotation, corresponding to “kinks” in the channel maps, are evident in the observations (panel A). In particular, the channelization of the data (in 200 m s⁻¹ bins) highlights that the isovelocity levels are distorted from the typical Keplerian butterfly pattern along the planetary wake. These perturbations coincide with those predicted from the semianalytic theory in the right panel of Figure 2 and in the predicted velocity map in panel B. Similar perturbations are also seen in the SPH model (panel C), in agreement with the analytic model. Plotting the corresponding spiral wake from Equation (3) in all three plots (solid line) reveals the extended planetary wake in the observational data.

A final bonus is that the spiral arm is not just seen as velocity perturbations in the observations, but also in intensity. Figure 5 demonstrates this, showing selected channels from the data. Beginning with the top-right panel, one may first observe the velocity kink associated with the spiral. The label indicates where the spiral arm can also be seen in intensity. The spiral continues into the preceding channels in the bottom-left and -right panels. The spiral is colocated with the predicted location of the outer planet wake generated by the embedded protoplanet, as shown by the dotted cyan line.

## 4. Discussion

We have demonstrated that several of the recently reported spiral structures in HD 163296 (Teague et al. 2021) are due to planet wakes generated by the embedded protoplanet (Pinte

et al. 2018a). These results provide the first direct confirmation of planetary spirals generated by Lindblad resonances (Goldreich & Tremaine 1979, 1980; Goodman & Rafikov 2001; Ogilvie & Lubow 2002; Rafikov 2002). We also show that the velocity-damping prescription adopted by B21 is not required. That is, the planetary wake indeed induces secondary kinks in velocity channels that extend far from the planet’s location.

### 4.1. Density Waves or Buoyancy Spirals?

We have demonstrated that kinematic perturbations from planetary companions are more extended than previously assumed (B21; Teague et al. 2021). The extended velocity perturbation associated with the embedded planet was also reported by Teague et al. (2021) after the subtraction of a flared Keplerian disk model. In their work, they suggest this is evidence of buoyancy spirals (Zhu et al. 2012; Bae et al. 2021) being excited by the embedded planet. Figure 2 and panel B of Figure 4 demonstrate that buoyancy spirals are not required to explain the secondary kinks (B21).

However, several additional substructures are seen in the channels of Figure 1 that do not appear to be associated with the planetary wake of the embedded protoplanet. These features were also reported by Teague et al. (2021). We labeled these features with semiopaque crosses in Figure 1. These additional features could possibly be attributed to buoyancy resonances from the embedded protoplanet. Buoyancy spirals are expected to produce the greatest perturbations close to the planet (Bae
<figure><img src="image_8.png" /><figcaption>**Figure 5.** A closer view of a selection of channels containing a faint spiral arm in CO emission. The location of the spiral is consistent with the projected location of the outer wake of the embedded protoplanet (dotted cyan line).</figcaption></figure>

et al. (2021) and with a smaller opening angle. An alternative possibility is that these additional features arise from a second companion at a smaller radial separation, as suggested by Teague et al. (2018). Pinte et al. (2020) reported the presence of a localized kink within the outermost gap of the continuum ring and interpreted this as a possible signature from an embedded planet. Teague et al. (2021) also found evidence of a velocity perturbation in this region. If such a planet does exist, we may be witnessing the outer wake of this planet.

Figure 3 confirms our hypothesis that the spirals are density waves rather than buoyancy spirals. Density waves are dominated by radial motions (Rafikov 2002), whereas buoyancy waves produce mostly vertical motions (Bae et al. 2021). As seen in Figure 2 (right panel) radial-velocity perturbations change signs along the semimajor axis, while vertical motions would not. Figure 3 shows that the sign flip expected for radial-dominated motions is present in the observational data. Another approach would be to look for this sign change in residual maps (Izquierdo et al. 2022; Teague et al. 2021), but we found that the sign flip measured this way is easily confused even in residuals from simulations, as it depends sensitively on the background disk model subtracted (see Appendix).

## 4.2. Estimating the Planet Mass

The observed planet wake in HD 163296 suggests a more robust way to estimate the properties of the embedded planet and the protoplanetary disk. As suggested by B21, the

amplitude of the velocity kink generated by the spiral as it crosses a channel is related to the perturber mass. The problem is that the velocity perturbation, and thus the kink amplitude, is related to the embedded planet mass in units of thermal mass as seen in Equation (9) of B21. The thermal mass is given by Goodman & Rafikov (2001),

$$M_{\text{th}} = \frac{2}{3} \left( \frac{h_p}{r_p} \right)^3 M_{\star}, \qquad (5)$$

and depends sensitively on the disk thermal structure. Thus, to measure the absolute planet mass one must break the degeneracy between the perturbation amplitude and the disk thermal structure (i.e., $h_p/r_p$) (see Bae & Zhu 2018a).

While good constraints on the disk temperature structure are possible from the CO emission (e.g., Pinte et al. 2018b), this does not trace the midplane gas temperature or the pressure scale height at the planet location, which is needed to measure the thermal mass. The CO emission between the upper and lower surface may allow for direct measurement of the midplane gas temperature (Dullemond et al. 2020), but the assumption of optically thick emission (for which $T_b = T_{\text{gas}}$) is not guaranteed to be valid due to the freeze-out of the molecules on the dust grains. Mapping the spiral wake provides a novel way to break this degeneracy and determine the perturber mass using only the disk kinematics. As seen in Equation (3), the propagation of the spiral wake depends on the
disk aspect ratio. Thus, by measuring both the amplitude of the kinks close to the planet (to determine the mass in units of the thermal mass) and fitting the wake shape to the location of the spiral-induced velocity kinks (to measure $h/r$ as a function of $r$), one can in principle directly determine the mass of the embedded planet using nonlinear wake propagation theory. Rabago & Zhu (2021) have also suggested constraining the planet mass using the maximum kink velocity.

We have shown that the shape of the planet wake in the data can be used to constrain the disk aspect ratio. This is needed for an accurate estimate of the thermal mass. Fitting the planet mass itself would require the simultaneous fitting of the kink amplitude, the wake shape, and the emitting surface, which is a complex procedure and beyond the scope of this Letter.

### 4.3. Limitations of the Models

Zhu et al. (2015) and Bae & Zhu (2018b) found that the wake predicted by Equation (3) underestimates the pitch angle of the spiral arm for large planet masses. Cimerman & Rafikov (2021) also found that the nonlinear planet wake propagation theory overpredicts the decay of the density wave. The perturbation amplitude near the planet may also be wrong for planet masses exceeding the thermal mass (Rafikov 2002). This is the reason we compared to both simulations and analytics, as the simulations do not have this limitation. Because we find similar results in both cases, this suggests the effects on the observed wake are limited.

In both the simulations and analytic models, we assumed a vertically isothermal temperature profile, which is not strictly true in the disk around HD 163296 (de Gregorio-Monsalvo et al. 2013; Rosenfeld et al. 2013). Temperature stratification would mean that the sound wave launched by the planet would propagate faster in the (hotter) upper layers of the disk and result in spirals with a larger pitch angle (Juhász & Rosotti 2018). Additional spirals due to buoyancy can also arise when there is a vertical temperature gradient (Bae et al. 2021). However, as seen in Figure 10 of Law et al. (2021), the temperature stratification is not significant in the outer regions of the disk up to several scale heights, which is also the region we are principally concerned with in this Letter. Additionally, the supplementary material of Pinte et al. (2019) shows similar morphology for the wake generated in both isothermal and stratified disks, although they only compared the regions close to the planet.

Additionally, wake propagation in 3D deviates from the linear wave theory of Equation (3). Zhu et al. (2015) showed that the pitch angle of the outer planetary wake is smaller (i.e., the spiral is more tightly wound) with increasing height above the disk midplane. This is taken into account in our SPH models but is not when we project the wake onto the CO-emitting layer. This effect will somewhat counteract the faster sound speed of the CO-emitting layer relative to the midplane. However, given that thermal stratification is not significant in the outer disk, this may be negligible. Additionally, shock heating also changes the wake shape and propagation (Zhu et al. 2015). A more robust treatment of the disk thermodynamics (e.g., accounting for radiative and shock heating) is required to explore this issue but is beyond the scope of this Letter.

## 5. Summary

We demonstrate that the recently reported spiral structure in HD 163296 (Teague et al. 2021) is the outer planetary wake generated by an embedded protoplanet (Pinte et al. 2018a), providing the first confirmation of planetary wakes generated by Lindblad resonances.

We thank Richard Teague and Andrés Izquierdo for useful discussion. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the US Department of Energy National Nuclear Security Administration under contract No. 89233218CNA000001. J.C. acknowledges the support from LANL/LDRD program (approved for public release as LA-UR-21-31084). D.J.P. and C.P. acknowledge funding from the Australian Research Council via FT130100034, DP180104235, and FT170100040. B.J.N. is supported by an Australian Government Research Training Program (RTP) Scholarship. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement 823823 (DUSTBUSTERS). We used PLONK (Mentiplay 2019), which utilizes functions written in SPLASH (Price 2007).

# Appendix
## Residuals from Keplerian Rotation

In Figure 6, we plot the Keplerian velocity residuals from our SPH before and after radiative transfer (left and middle), along with the velocity residuals from Teague et al. (2021). We used the code BETTERMOMENTS to generate the velocity map of our SPH+RT model and used the quadratic method (Teague & Foreman-Mackey 2018). We follow the fitting procedure outlined in Section 3.1 of Teague et al. (2021) to fit our radiative transfer model and obtain a best-fit CO-emitting layer, stellar mass, and disk position angle. The left panel of Figure 6 shows the deviations from Keplerian rotation of the SPH particles close to this layer (accounting for their height above the disk midplane).

An issue with this approach is that EDDY (Teague 2019), when fed the synthetic moment map from our SPH simulation (middle panel in Figure 6), fits a stellar mass of $1.84 M_\odot$, whereas the mass of the central object in our model was $1.9 M_\odot$. This difference is likely due to the gas pressure support not being accounted for in the fitting model. However, the resulting error in the residuals from the incorrect stellar mass is $v_{res} = \sqrt{G\Delta M_*/r_p} \approx 460 \text{ m s}^{-1}$ with $\Delta M_*=0.06 M_\odot$ and $r_p=250 \text{ au}$—on the same order as the residuals of $\approx 300 \text{ m s}^{-1}$ from the planetary wake (Figure 2). These rather minor errors in the stellar mass determination produce residuals that do not reflect those in the SPH model. A similar effect could also be generated by small changes in the CO-emitting layer. Because these small differences in model parameters produce velocity deviations on the order of the velocity deviations induced by the kink, they can produce a residual map that hides the expected kink signature seen in the left panel of Figure 6 and the right panel of Figure 3. This also explains why residual maps produced from either our simulations or the observations (Teague et al. 2021) do not obviously show the sign flip across the major axis expected for radial motions.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Residuals from Keplerian rotation predicted from our SPH model (left panel), from our SPH model after radiative transfer (middle panel), and from the observations presented in Teague et al. (2021). Our SPH model shows radial, vertical, and azimuthal velocity deviations from circular Keplerian rotation. We used EDDY to subtract a flared Keplerian disk model from our SPH+RT model. EDDY does not recover the true residuals in our simulation because of errors in the stellar mass and CO-emitting layer determination, but the residual map contains all the major features of the residual map produced in the same manner from the observations (shown in Figure 1(b) of Teague et al. 2021).</figcaption></figure>

## ORCID iDs

Josh Calcino https://orcid.org/0000-0001-7764-3627

Thomas Hilder https://orcid.org/0000-0001-7641-5235

Daniel J. Price https://orcid.org/0000-0002-4716-4235

Christophe Pinte https://orcid.org/0000-0001-5907-5179

Giuseppe Lodato https://orcid.org/0000-0002-2357-7692

Brodie J. Norfolk https://orcid.org/0000-0001-5898-2420

## References

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

Bae, J., & Zhu, Z. 2018a, ApJ, 859, 119

Bae, J., & Zhu, Z. 2018b, ApJ, 859, 118

Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362

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

Cimerman, N. P., & Rafikov, R. R. 2021, MNRAS, 508, 2329

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApjS, 257, 2

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2

Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Juhász, A., & Rosotti, G. P. 2018, MNRAS, 474, L32

Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4

Mentiplay, D. 2019, JOSS, 4, 1884

Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

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

Pinte, C., Ménard, F., Duchêne, G., et al. 2018b, A&A, 609, A47

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

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

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

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

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

Rafikov, R. R. 2002, ApJ, 569, 997

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Setterholm, B. R., Monnier, J. D., Davies, C. L., et al. 2018, ApJ, 869, 164

Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344

Teague, R. 2019, JOSS, 4, 1220

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJL, 860, L12

Teague, R., & Foreman-Mackey, D. 2018, Zenodo, Zenodo doi:10.5281/ zenodo.1419754

Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJL, 758, L42