# Using planet migration and dust drift to weigh protoplanetary discs

Yinhao Wu (吴寅昊)<fnref n="1" />,*, Clément Baruteau <fnref n="2" /> and Sergei Nayakshin <fnref n="1" />

Accepted 2023 June 12. Received 2023 June 9; in original form 2023 January 16

## ABSTRACT

ALMA has spatially resolved over 200 annular structures in protoplanetary discs, many of which are suggestive of the presence of planets. Constraining the mass of these putative planets is quite degenerate for it depends on the disc physical properties, and for simplicity a steady-state is often assumed whereby the planet position is kept fixed and there is a constant source of dust at the outer edge of the disc. Here, we argue against this approach by demonstrating how the planet and dust dynamics can lift degeneracies of such steady-state models. We take main disc parameters from the well-known protoplanetary disc HD 163296 with a suspected planet at $R \approx 86$ au as an example. By running gas and dust hydrodynamical simulations post-processed with dust radiative transfer calculations, we first find steady-state disc and planet parameters that reproduce ALMA continuum observations fairly well. For the same disc mass, but now allowing the planet to migrate in the simulation, we find that the planet undergoes runaway migration and reaches the inner disc in $\sim 0.2$ Myr. Further, decreasing the disc mass slows down planet migration, but it then also increases the dust’s radial drift, thereby depleting the disc dust faster. We find that the opposing constraints of planet migration and dust drift require the disc mass to be at most $0.025 M_\odot$, must less massive than previously estimated, and for the dust to be porous rather than compact. We propose that similar analysis should be extended to other sources with suspected planetary companions.

**Key words:** planets and satellites: formation–planet–disc interactions – protoplanetary discs.

## 1 INTRODUCTION

Protoplanetary discs are the engines that drive formation of young protostars as well as formation and early orbital evolution of planetary systems. Mass is one the most fundamental characteristics of protoplanetary discs, yet it remains very difficult to constrain observationally. This is mainly because $H_2$, by far the most dominant constituent of protoplanetary discs, emits very weakly in the range of temperatures encountered at radii where most of the disc mass is contained (Carmona 2010). Consequently, one needs to rely on two main approaches to estimate the disc mass: emission by minor constituents (gas species or dust), and discs dynamics.

In the first approach, dust particle continuum emission at mm- wavelengths has been used widely to estimate disc mass. Provided that dust emission is optically thin, that dust scattering is negligible, that the dust’s absorption opacity and temperature are uniform and known, disc dust mass can be inferred from the averaged flux of dust continuum emission for a source with known distance to Earth. The total mass of the disc can then be obtained by further assuming the dust-to-gas mass ratio (see Miotello et al. 2022 for a recent review). Detailed modelling shows that uncertainties in the assumptions mentioned above supplemented by modelling freedom in dust particle size distribution, composition, and porosity (e.g. Kataoka et al. 2014; Woitke et al. 2016), result in order of magnitude

uncertainties in the inferred disc masses (e.g. see Andrews et al. 2018a; Zhu et al. 2019; Miotello et al. 2022; Xin et al. 2023).

The emission of minor gas constituents, like CO and its isotopo- logues, is an alternative route to tracing the discs mass, but it also has its own flaws, including optical depth effects and freeze-out on to dust grains (Bergin et al. 2013; Kamp et al. 2017). Because it does not freeze out and its abundance relative to $H_2$ is well known, HD is a promising disc mass tracer despite the degeneracy between its abundance and the gas temperature (Bergin et al. 2016; Trapman et al. 2017; Schwarz et al. 2021), however, only a handful protoplanetary discs have HD emission data currently. All in all, the estimation of disc masses based on dust and/or gas emission is probably uncertain by an order of magnitude (Miotello et al. 2022).

In the second approach based on discs dynamics, several methods have been proposed. First, the disc mass can be crudely estimated via the observed stellar accretion rate ($\dot{M}_*$) and the age of the star ($t_*$) as $M_{\text{disc}} \approx \zeta \dot{M}_* t_*$, where $\zeta$ is a few (Jones, Pringle & Alexander 2012). This assumes a quasi-steady state evolution for the disc gas driven by radial (turbulent) transport of angular momentum, an assumption that is more and more challenged by the emergent picture of vertical transport due to magnetized winds (Suzuki et al. 2016; Lesur 2021; Elbakyan et al. 2022; Tabone et al. 2023; Wu et al. 2023). A second dynamics method is based on dynamics of dust particles rather than that of gas. So long as dust particles have a Stokes number $St \ll 1$, their radial drift velocity due to gas drag is directly proportional to the dust particle size ($s$) and inversely proportional to the gas density (Birnstiel, Klahr & Ercolano 2012). By observing a disc at different wavelengths, one could in principle infer the radial location of the

<footnote n="1">¹School of Physics and Astronomy, University of Leicester, Leicester LE1 7RH, UK</footnote>

<footnote n="2">²IRAP, Université de Toulouse, CNRS, UPS, CNES, Toulouse, France</footnote>

<footnote>*E-mail: yw505@leicester.ac.uk</footnote>
outer edge of the dust in the disc as a function of $s$, which can then be turned into an estimate of the gas disc mass (Powell, Murray-Clay & Schlichting 2017; Powell et al. 2019). Unfortunately, uncertainties in the dust opacity, and the fact that dust porosity strongly affects the radial dust drift rates at high porosity, preclude this method from being sufficiently constraining (Franceschi et al. 2022). Note also that both aforementioned methods assume the age of the system is known, while it is usually poorly constrained for young stellar objects. Another method measures departure of a disc rotation curve from its Keplerian profile, but this is well suited only to very massive self-gravitating discs (Veronesi et al. 2021).

In this paper, we propose another dynamical method to estimate the discs mass. It is based on planetary migration in discs for which dust sub-structures, in particular dark rings in the dust continuum emission, are assumed to have a planetary origin. In less than a decade, spatially resolved observations of protoplanetary discs have uncovered a wealth of sub-structures in both dust and gas emission that much resemble those seen in simulations of disc-planet interactions (e.g. ALMA Partnership et al. 2015; Huang et al. 2018; Long et al. 2018; Zhang et al. 2018; Casassus & Pérez 2019; Pérez et al. 2019; Tsukagoshi et al. 2019; Pinte et al. 2020). The planetary interpretation of these structures can (and should) be debated as other planet-free mechanisms may also form similar structures in the discs emission (Meru et al. 2017; Riols, Lesur & Menard 2020). Still, if discs structures are due to planets, then it is possible to constrain the mass of these planets by comparing synthetic images of the discs emission with real observations. In doing this comparison, synthetic images are most often based on gas and dust hydrodynamical simulations of disc-planet(s) interactions in which planetary migration is discarded for simplicity (e.g. Zhang et al. 2018). This assumption is valid if planet migration time-scale is long in comparison to the characteristic time-scale for evolution of gas and dust in the disc, a lower limit of which being set by dust particles radial drift. Since the speed of planetary migration increases with both planet mass and disc mass (or, actually, with the disc's surface density), assuming that disc structures are caused by planets on fixed orbits may be reasonable for some system parameters but not for others. A recent illustration of this point has been presented in Wafflard-Fernandez & Baruteau (2020). They found that a single planet, with a mass typically between that of Saturn and Jupiter, can undergo rapid intermittent type III migration provided that the disc is massive enough. This unsteady migration pattern may trigger formation of multiple dust rings and gaps which in a fixed planet orbit scenario would require multiple planets. It is also known that migrating planets produce different effects on the structure of the gaps and rings they carve in the disc compared with static planets (Meru et al. 2019; Pérez et al. 2019; Weber et al. 2019). It is possible that dust continuum emission in planet vicinity can serve as a handle on the unknown disc mass.

Furthermore, in this paper we point out the following simple constraint that planet migration and dust dynamics place on the disc mass when taken together. It is well known that even for low mass planets, in the Neptune to Saturn mass range, the migration time-scale may be quite a small fraction of the disc lifetime $t_*$ at high disc masses (Baruteau et al. 2014). On the other hand, to have a decent statistical chance of catching a migrating planet in the process of structuring its dusty disc we would expect its migration time-scale, $t_m$, to be a fair fraction of the disc age. Indeed, $t_m$ is the rough measure of how long it takes for the planet to disappear from the outer disc into the inner unresolved one, so if we happen to image a planet in the outer disc then it is fair to assume that it is not older than $\sim t_m$. But if planets form at random times with respect to the time when

we image the discs then the probability of observing a given planet is $\sim t_m/t_*$. We thus argue that a maximum mass for structured discs can be set in order to avoid too fast planetary migration. While only an upper disc mass can be set with this approach, it may turn out quite useful in conjunction with the other methods discussed above. For example, the dust lines method of Powell et al. (2017, 2019), gives a lower limit on the disc mass. By combining dust and planet dynamics the discs mass could hopefully be constrained from both below and above.

The present study aims at testing these ideas via multifluid hydrodynamical simulations of disc-planet interactions post-processed by dust radiative transfer calculations. Our disc model is inspired by the well-known disc around the $M_* \approx 2 M_\odot$ Herbig Ae star HD 163296, a relatively old protostar with an estimated age $t_* \approx (7-10)$ Myr (e.g. Setterholm et al. 2018; Vioque et al. 2018). The disc has been observed with ALMA as part of the DSHARP Large Programme (Huang et al. 2018; Andrews et al. 2018b). Dust continuum emission from the system features at least 4 dark rings including two prominent ones located around 48 and 86 au (Isella et al. 2018). Localised velocity perturbations based on CO gas emission strongly suggest the existence of several massive planets in the disc (Teague et al. 2018; Pinte et al. 2020), including one inside the dark ring around 86 au, which our study will focus on. Especially relevant to our study is the fact that the mass of the disc in HD 163 296 has been estimated by various techniques: through dust continuum emission ($\approx 0.12 M_\odot$, Powell et al. 2019), dust dynamics ($\approx 0.2 M_\odot$, Powell et al. 2019) and also via chemodynamical modelling of the emission of various CO isotopologues and of HD (between 0.05 and $0.3 M_\odot$, Booth et al. 2019; Kama et al. 2020).

This paper is organized as follows. Section 2 details our disc model inspired from the HD 163 296 disc, and describes the setup of our hydrodynamical simulations and dust radiative transfer calculations. Section 3 presents our results of simulations and dust radiative transfer calculations. A few concluding remarks follow in Section 4.

## 2 DISC MODEL AND NUMERICAL METHODS

As already stated in the introduction, here we focus on HD 163 296 system as a case study. We focus specifically on the dark ring in the dust continuum emission at around 86 au that has been suggested to be carved by a massive planet (e.g. Liu et al. 2018; Garrido-Deutelmoser et al. 2023). Detection of a localized velocity perturbation from the dark ring's vicinity (Teague et al. 2018; Pinte et al. 2020) makes the case for a massive planet especially compelling. We emphasize that we do not aim here to reproduce all the structural features observed in the disc of HD 163296, but just use the disc parameters as an example. For a recent example of a study setting that aim we point the reader to Garrido-Deutelmoser et al. (2023), which came out while this manuscript was in review. The authors use four massive planets held on fixed orbits in a massive gas disc to try and reproduce HD 163 296 observations in detail. Again, our goal here is to interpret the dark ring in HD 163 296 as a dust gap carved by a single massive planet, and to explore what planet and disc masses (among other physical parameters) could account for the ring properties when planet migration is taken into account. To this aim, we have carried out two-dimensional (2D) hydrodynamical simulations of disc-planet interactions to model the gas and dust evolution of a disc whose physical properties are chosen to reflect those in the HD 163 296 disc (see Section 2.1). The dust surface density obtained in the hydrodynamical simulations is then used to calculate synthetic maps of the dust continuum emission in
the sub-millimetre via 3D dust radiative transfer calculations (see Section 2.2).

Our work is split into two stages. In the first one, we assume that the planet remains on a fixed circular orbit and we examine what disc and planet parameters can account for the flux properties of the dark ring of emission near 86 au in the HD 163 296 disc. Having a fixed planet is a common simplifying assumption that has often been used in previous studies to infer what planet mass would be needed to explain the width and depth of dark rings of emission in the sub-mm. In the second stage of our study, we consider the radial migration of the planet, given the disc and planet parameters obtained in the first stage, and we examine how planet migration affects our findings.

## 2.1 Setup of hydrodynamical simulations

We use the public grid-based code FARGO3D (Benítez-Llambay & Masset 2016) in its multifluid version (Benítez-Llambay, Krapp & Pessah 2019) to model the evolution of a protoplanetary disc with an embedded planet. The disc model includes both gas and dust, which interact via linear drag forces. Planet-disc interaction due to gravity is included, but gas self-gravity is neglected. Our simulations are 2D and in the following, we describe the disc by its polar coordinates $R$ and $\phi$. Since the star remains at the frame's origin, the indirect terms arising from the planet and the disc accelerations are fully taken into account.

### 2.1.1 Gas model

We assume that the evolution of the disc gas is driven by radial transport of angular momentum via a turbulent viscosity. We set the turbulent viscosity to the form of the isotropic prescription by Shakura & Sunyaev (1973) with a constant $\alpha$ viscosity parameter. ALMA CO and DCO$^+$ molecular line observations of the HD 163 296 disc only set an upper limit to the disc turbulence parameter $\alpha \lesssim 3 \times 10^{-3}$ (Flaherty et al. 2017). We have therefore experimented with two values of the $\alpha$ parameter: $10^{-4}$ and $10^{-3}$ (Liu et al. 2022).

Because the total mass in the HD 163 296 disc is not constrained accurately neither, we have also varied the disc mass in our simulations, which allowed us to study its impact on the planet migration rate as well as on the dust drift rate. The disc mass $M_{\text{disc}}$ ranges from 0.0125 $M_{\odot}$ to 0.2 $M_{\odot}$. Following Dullemond et al. (2020), the initial gas surface density profile at time $t=0$ in the simulation is given by

$$ \Sigma_g(R, t=0) = \Sigma_{g,0} \left(\frac{R}{R_0}\right)^{-p} \exp\left[-\left(\frac{R}{R_{\text{exp}}}\right)^{2-p}\right], \quad (1) $$

with $R_0 = 100$ au, $R_{\text{exp}} = 150$ au, $p=1$, and $\Sigma_{g,0}$ is a normalization constant that changes with $M_{\text{disc}}$.

A locally isothermal equation of state is adopted for simplicity, meaning that the gas temperature stays stationary. This assumption is most presumably valid given that we model regions in the HD 163 296 disc located around 100 au, where the cooling time-scale is expected to be much shorter than the dynamical time-scale. We set the radial profile of the gas temperature, $T(R)$, as

$$ T(R) \approx 15 \left(\frac{R}{R_0}\right)^{-0.14} K, \quad (2) $$

which features the same power-law exponent as in the mid-plane temperature inferred from Dullemond et al. (2020) by the modelling of CO(2-1) emission in the disc. Our reference temperature at $R_0 = 100$ au is, however, slightly smaller than their mid-plane temperature at the same location (~25 K). This is motivated by preliminary

simulations and dust radiative transfer calculations, which showed that, for the range of planet masses that we have experimented, a smaller disc temperature was required to obtain a narrower gap around the planet's location, closer to that in the ALMA continuum observation (see also Dong et al. 2018). This is similar to the simulations presented in Pérez et al. (2019), where a disc temperature quite smaller than what would be usually assumed was necessary for a gap-opening planet to account for the short radial separation between the three fine bright rings in the outer parts of the disc around HD 169 142 as observed with ALMA. Assuming a Solar-like composition for the disc gas (mean molecular weight equal to 2.4), our temperature profile translates into the following radial profile for the disc's aspect ratio, $h(R)$:

$$ h(R) = \frac{c_s}{v_K}(R) = 0.054 \times \left(\frac{R}{R_0}\right)^{0.43}, \quad (3) $$

where $v_K$ is the Keplerian velocity in which we have assumed the stellar mass to be $M_{\star} = 2.0 M_{\odot}$ (Andrews et al. 2018b).

### 2.1.2 Dust model

Dust is modelled as multiple pressure-less fluids, and our simulations include dust diffusion (Weber et al. 2019) as well as dust feedback on the gas. We have modified the public source files of FARGO3D such that each dust fluid corresponds to a fixed particle size instead of a fixed Stokes number. The simulations presented in this paper use 5 dust fluids with size $s$ spaced logarithmically uniformly between 30 µm and 3 mm. We adopt an initial size distribution in $s^{-3.5}$ and an initial dust-to-gas mass ratio (and surface density ratio) that defaults to 1 per cent. In the Epstein regime, which is relevant to our disc model, the Stokes number can be approximated as

$$ \text{St} \approx \pi/2 \times s \rho_d / \Sigma_g, \quad (4) $$

where $\rho_d$ is the internal mass volume density for each dust fluid and $\Sigma_g$ is the local gas surface density. By default, our simulations use $\rho_d = 1.3 \text{ g cm}^{-3}$, which is also consistent with the dust composition assumed in our dust radiative transfer calculations. Simulations with smaller internal densities will be presented in Section 3.5.

### 2.1.3 Planet setup

Dong et al. (2018) showed how the radial location of the maxima and minima in the dust's surface density around a gap-opening planet depends on the planet mass $M_p$ and the disc's aspect ratio for inviscid discs. Following these authors, and given our disc's aspect ratio, we estimate that $M_p$ should be $\approx 0.26 M_J$ to approximately match the radial width of the dark ring in the mm-emission that is located around 86 au in the HD 163 296 disc. We therefore set the planet-to-star mass ratio to $1.3 \times 10^{-4}$. Given that this value is inferred from inviscid disc models, while our simulations do include viscosity, we also tried lower and higher planet masses (see Appendix A). In order for the disc to accommodate the disturbance generated by the planet smoothly, we let the planet's mass increase from zero to its final mass according to the formula

$$ M_{\text{taper}} = M_p \sin^2 \left( \frac{\pi}{2} \times \frac{t}{t_{\text{taper}}} \right), \quad (5) $$

with $t_{\text{taper}} = 10$ orbits is the time-scale for the mass of planet to increase from 0 to $M_p$.

In our simulations, we either hold the planet on a fixed circular orbit at 86 au or we let it migrate freely in the disc. In the latter case,
we consider two cases: (i) the planet starts migrating at $R_{p, start} = 86$ au after having remained on a fixed orbit for 100 orbits, and (ii) the planet starts migrating at $R_{p, start} = 130$ au after having remained on a fixed orbit for only 10 orbits. Since neither gas nor dust self-gravity is included in our simulations, the total disc's surface density is subtracted from its azimuthal average prior to the calculation of the disc force on the planet (Baruteau & Masset 2008; Benítez-Llambay & Masset 2016), and the force exerted by the gas and dust inside the planet's Hill radius is discarded (Crida et al. 2009). To mimic the effects of finite vertical thickness, the planet's gravitational potential is smoothed over a softening length set to 0.6 pressure scale heights at the planet's orbital radius (Müller, Kley & Meru 2012).

### 2.1.4 Numerical setup

We use 512 grid cells in the radial direction with a linearly uniform radial spacing, and 1536 in the azimuthal direction. The grid extends from 0.3 $R_0$ (30 au) to 2.2 $R_0$ (220 au) in the radial direction, and from 0 to $2\pi$ in the azimuthal direction. Regarding boundary conditions, for the gas we use so-called wave-damping boundary conditions, whereby the gas density and velocity components are damped towards their initial radial profile. For the dust fluids, we use an open boundary condition for their radial velocity at the grid's inner edge, and we do not use wave-damping. This implies that the surface density of the dust fluids at the grid's outer edge is not kept constant over our simulations, which means there is a finite supply of dust from outside of the system.

## 2.2 Setup of dust radiative transfer calculations

Our dust radiative transfer calculations were carried out with the public code RADMC-3D (Dullemond et al. 2012). The post-processing of our FARGO3D simulations as inputs to RADMC-3D, and the calculation of the final beam-convolved synthetic maps of dust continuum emission were done with the public PYTHON code fargo2radmc3d. The setup of the dust radiative calculations is very similar to that described in Baruteau et al. (2021), in particular the fact that the dust temperatures used in RADMC-3D correspond to the temperature used in the hydrodynamical simulations. We therefore only briefly list the main parameters used. The surface density of each of our 5 dust fluids is projected onto a 3D spherical grid assuming vertical hydrostatic equilibrium and dust scale heights set to

$$h_{i, \text{dust}} = h(R) \times \left( \frac{\alpha}{\alpha + \text{St}_i} \right)^{1/2}, \quad (6)$$

where $i = 1..5$ and $\text{St}_i$ is the average Stokes number for the $i$th dust fluid. The spherical grid covers about $\pm 3H$ over the mid-plane over 40 cells with a non-uniform spacing. Henyey-Greenstein anisotropic dust scattering is included. The dust absorption and scattering opacities are calculated according to Mie theory, assuming a default dust composition comprised of 70 per cent water ices and 30 per cent astro silicates, thereby corresponding to an internal mass volume density of about $1.3 \text{ g cm}^{-3}$ (the same value as in the hydrodynamical simulations). Note that the results of hydrodynamical simulations and dust radiative transfer calculations presented in Section 3.5 were obtained with porous dust having an internal density of $0.1 \text{ g cm}^{-3}$, assuming the same composition as above but an additional volume fraction of vacuum set to $\sim 92$ per cent (Baruteau et al. 2019).

The specific intensity of the dust continuum emission is computed at a wavelength of 1.25 mm with $10^8$ photon packages used for ray-tracing. The disc is taken to have a distance of 101 pc, an

inclination of 46.7 deg and a position angle of 133 deg. To produce the final synthetic image of the disc, we first add white noise to the map of specific intensity computed by RADMC-3D (prior to beam convolution) with an rms noise level of $2.3 \times 10^{-5} \text{Jy beam}^{-1}$, and the obtained map is finally convolved with an elliptical beam of $0.048 \times 0.038 \text{ arcsec}^2$ and position angle of 81.7 deg (Isella et al. 2018).

# 3 RESULTS

## 3.1 Simulations table and guide to main results

In this Section, we present results of our hydrodynamical simulations and of dust radiative transfer calculations. For reader's convenience in Table 1 we summarize the parameters and key results of our simulations. The table is broken into four horizontal blocks separated by double lines, with model names reflecting one or more major model assumptions. In the first section, we present the case of a planet on a fixed orbit at 86 au (Section 3.2). These model names start with 'NM' to emphasize that they include no migration for the planet. The purpose of this set of models is to constrain the disc's mass and turbulent viscosity that yield synthetic images in a reasonable agreement with the observed ones. Having constrained in this way viscosity parameter, $\alpha = 10^{-3}$, we use this value for the rest of the paper.

In Section 3.3, we then relax the assumption of a fixed orbit for the planet, allowing it to migrate freely. The names of these models in Table 1 start with 'O' to indicate that the planet in our model starts from the 'observed' position of the dark ring of continuum emission at 86 au (Huang et al. 2018). The models in this section differ by the disc mass, which we vary in steps of a factor of 2 from the low value of $M_{\text{disc}} = 0.0125 \text{ M}_\odot$ in Model OM0125 to the high value of $M_{\text{disc}} = 0.2 \text{ M}_\odot$ in Model OM2. The main purpose of these experiments is to demonstrate how planet migration can constrain the otherwise unknown disc mass $M_{\text{disc}}$; in the case of HD 163 296 we find $M_{\text{disc}} \le 0.025 \text{ M}_\odot$ for our fiducial planet mass. This value for the disc mass challenges the conclusion of the fixed planet simulations in Section 3.2, where we find a good match to the data with $M_{\text{disc}} = 0.05 \text{ M}_\odot$.

For a more comprehensive comparison of our simulations with the dust continuum observations, in Section 3.4 we let the planet to start much further out, at $R_{p, start} = 130$ au and run the simulations until the planet reaches $R_{p, end} = 86$ au. The names of these models in Table 1 start with 'P' to indicate that the planet starts from a 'past' position – it is unlikely that the planet currently observed at 86 au formed there without undergoing migration from further out. Similar to Section 3.2 we perform dust radiative transfer calculations to build synthetic images when the planet reaches 86 au. The comparison of the synthetic images with the observations shows that there is no disc mass that simultaneously satisfies the opposing constraints of dust drift and planet migration. At high disc masses the planet reaches the inner parts of our computational grid too rapidly, whereas at low disc masses the planet migration is reasonably slow but the dust reaches the inner parts of the grid too rapidly. This conclusion therefore holds whether we start the planet at 86 au or let it migrate from further out.

Therefore, in Section 3.5 we propose that the dust in the HD 163 296 disc may be more porous than assumed in Sections 3.2–3.4. More specifically, we set the internal dust density to $\rho_d = 0.1 \text{ g cm}^{-3}$ in both the hydrodynamical simulations and the dust radiative transfer calculations, and the names of these models in Table 1 end with 'p' to reflect the porous dust assumption. We then find that porous dust does mitigate the opposing challenges of the too rapid dust
<table><caption>**Table 1.** List of models adopted in the hydrodynamical simulations. The main parameters are listed in Columns 2–6: the disc mass $M_{\text{disc}}$, the planet mass $M_{\text{p}}$, the internal mass volume density $\rho_{\text{d}}$ for the dust fluid, whether the planet is allowed to migrate or not, and the initial orbital radius of the planet. Column 7 indicates the migration time-scale $t_{\text{m}}$ equation (7) as measured from the simulations. The last column uses a comment to characterize main results of the simulations.</caption><thead><tr><th rowspan="2">Model</th><th colspan="6">Hydrodynamical simulations models: parameters, results, and comments</th></tr><tr><th><i>M</i><sub>disc</sub><br>(M<sub>☉</sub>)</th><th><i>M</i><sub>p</sub><br>(M<sub>☉</sub>)</th><th><i>ρ</i><sub>d</sub><br>(g cm<sup>-3</sup>)</th><th>Migration<br>(Yes/No)</th><th><i>R</i><sub>p, start</sub><br>(au)</th><th><i>t</i><sub>m</sub><br>(kyr)</th><th>Comments</th></tr></thead><tbody><tr><td>NM05</td><td>0.05</td><td>0.26</td><td>1.3</td><td>No</td><td>86</td><td></td><td>Best-fitting model w/o planet migration</td></tr><tr><td>NM05a4</td><td>0.05</td><td>0.26</td><td>1.3</td><td>No</td><td>86</td><td></td><td>α = 10<sup>-4</sup>, Non-axisymmetric rings</td></tr><tr><td>NM025</td><td>0.025</td><td>0.26</td><td>1.3</td><td>No</td><td>86</td><td></td><td>Too wide gap</td></tr><tr><td>NM05H</td><td>0.05</td><td>0.52</td><td>1.3</td><td>No</td><td>86</td><td></td><td>Too wide gap (Appendix A)</td></tr><tr><td>NM05L</td><td>0.05</td><td>0.13</td><td>1.3</td><td>No</td><td>86</td><td></td><td>Good fit model (Appendix A)</td></tr><tr><td>OM0125</td><td>0.0125</td><td>0.26</td><td>1.3</td><td>Yes</td><td>86</td><td>3612</td><td>Slow migration</td></tr><tr><td>OM025</td><td>0.025</td><td>0.26</td><td>1.3</td><td>Yes</td><td>86</td><td>1868</td><td>Slow migration</td></tr><tr><td>OM05</td><td>0.05</td><td>0.26</td><td>1.3</td><td>Yes</td><td>86</td><td>425</td><td>Runaway migration</td></tr><tr><td>OM1</td><td>0.1</td><td>0.26</td><td>1.3</td><td>Yes</td><td>86</td><td>190</td><td>Runaway migration</td></tr><tr><td>OM2</td><td>0.2</td><td>0.26</td><td>1.3</td><td>Yes</td><td>86</td><td>118</td><td>Runaway migration</td></tr><tr><td>PM025</td><td>0.025</td><td>0.26</td><td>1.3</td><td>Yes</td><td>130</td><td>1985</td><td>No outer bright ring</td></tr><tr><td>PM05</td><td>0.05</td><td>0.26</td><td>1.3</td><td>Yes</td><td>130</td><td>869</td><td>Weak outer bright ring</td></tr><tr><td>PM05H</td><td>0.05</td><td>0.52</td><td>1.3</td><td>Yes</td><td>130</td><td>538</td><td>Too wide gap (Appendix A)</td></tr><tr><td>PM05L</td><td>0.05</td><td>0.13</td><td>1.3</td><td>Yes</td><td>130</td><td>1282</td><td>No outer bright ring (Appendix A)</td></tr><tr><td>PM025p</td><td>0.025</td><td>0.26</td><td>0.1</td><td>Yes</td><td>130</td><td>2006</td><td>Insufficient flux of continuum emission</td></tr><tr><td>PM025d10p</td><td>0.025</td><td>0.26</td><td>0.1</td><td>Yes</td><td>130</td><td>2006</td><td>ε = 0.1, Best-fitting model w/o planet migration</td></tr></tbody></table>

drift and planet migration but then results in an insufficient dust continuum flux in the synthetic images. We finally explore the idea that the disc in HD 163 296 may be much more dust rich than usually assumed, by increasing the initial dust-to-gas mass ratio to 0.1 in Model PM025d10p.

## 3.2 Reference case: a planet on a fixed orbit

The disc's turbulent viscosity is a major unknown that affects the depth of a gap opened by a planet of a given mass in a disc with a given aspect ratio. Recall that in our disc models, $h(R)$ is given by equation (3) and our fiducial planet-to-star mass ratio is $q = 1.3 \times 10^{-4}$. Following previous works (e.g. Liu et al. 2018; Zhang et al. 2018), we constrain the disc's $\alpha$ turbulent viscosity by comparing the shape and width of the dark ring in dust continuum emission caused by the planet in our simulations with those observed in HD 163 296 around 86 au. As pointed out in Section 2.1.1, an upper limit to $\alpha$ from molecular line observations is $\sim 3 \times 10^{-3}$. Therefore, we tested two values of the $\alpha$ turbulent viscosity: $10^{-3}$ (models NM05 and NM025) and $10^{-4}$ (Model NM05a4).

For detailed analysis, we present results at 400 orbits at $R_0$ ($\sim 2.8 \times 10^5$ yr with $R_0 = 100$ au). This corresponds to just about 500 orbits at the planet's orbital radius at 86 au. We find that at this time in the simulation the width and depth of the gap in the gas reach a quasi steady-state.<fnref n="1" /> We then performed dust radiation transfer calculations

and computed synthetic images of the dust continuum emission at 1.25 mm. The side-by-side comparison between our synthetic images and the DSHARP band 6 continuum data for HD 163 296 is presented in Fig. 1. Note that our radiative transfer calculations use a 3D spherical grid that extends the 2D polar grid of our hydrodynamical simulations. The grid is therefore truncated in the radial direction, which is why our synthetic maps feature an inner hole about the central star (it is thus not due to a real lack of dust material in our simulations).

Model NM05a4, which is shown in the lower left panel of Fig. 1, appears to be the least promising fit to the observed emission map, since there is a pronounced vortex in the synthetic data. This large-scale vortex, which is located at the outer edge of the planet gap, forms due to the Rossby-Wave Instability setting in at that location (e.g. Li et al. 2001), and corresponds to the outermost asymmetric bright ring in the synthetic image. The low viscosity also leads to some dust still corotating with the planet at the time in the simulation, which is not seen in the observed image.

The synthetic image of Model NM025 is shown in the lower right panel of Fig. 1. For this model, we use $\alpha = 10^{-3}$ and $M_{\text{disc}} = 0.025 M_{\odot}$. No vortices appear in this case, but the gap carved by the planet in this model is much wider than the one seen in the observation. We also note that the disc region interior to the planet is dimmer than in the other models. The reason is two-fold. First, a lower disc mass implies a lower dust mass for a fixed dust-to-gas mass ratio. This contributes to lowering the overall dust emission, which turns out to be mostly optically thin in this model. Also, a lower disc mass implies faster dust drift, in particular for the 1 mm-

orbits, and found no change in our qualitative results. The only change is the amount of dust around L4, L5 points is smaller due to dust diffusion (e.g. Zhang et al. 2018; Wu et al. 2023). Finally, the low viscosity discs are not the focus of our study, but the fact that low viscosity can generate different substructures in disc will help the understanding of our results.

<footnote n="1">¹Some studies, e.g. Kanagawa et al. (2016), find that much longer integration times are needed to reach steady state if the disc viscosity is as low as $\alpha = 10^{-4}$. However, here we find that for such low viscosities synthetic images show significant deviations from the observed one already at early times. Evolving the simulations for longer would only increase the disagreement between synthetic and observed images. Finally, we find that the system evolves after the planet is allowed to migrate on time-scales as short as $O(100)$ orbits, which is much shorter than the relaxation time we use. To check this conclusion we extended the low viscosity simulation to more</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Observed Band 6 continuum image (left-upper) of the HD 163 296 disc obtained from the fiducial images data on the DSHARP data release website, compared with the synthetic maps at the same wavelength obtained for our 3 models with a planet fixed at 86 au (the model's name is shown in the left-upper corner of each figure). The white × symbol represents the position of the planet. The beam size is shown as a grey ellipse in the left-lower corner of each figure.</figcaption></figure>

sized dust which is therefore more depleted in the inner disc parts in Model NM025 compared to Model NM05 (this particular point will be illustrated in the left-hand panels in Fig. 4 for Model PM025, which includes planet migration).

The upper right panel of Fig. 1 shows the synthetic image obtained with Model NM05 for $\alpha = 10^{-3}$ and $M_{\text{disc}} = 0.05 M_{\odot}$. This model provides a closer match to the observations than previous Models NM025 and NM05a4. The synthetic map of Model NM05 does not show asymmetric structures as in the lower viscosity model, and its gap width and the mean dust continuum flux are more consistent with the observation.

Concluding this section, we see that with a fixed planet model it is possible to obtain a reasonable synthetic match to the data by adjusting the disc mass and its viscosity parameter. For the problem at hand, $\alpha = 10^{-3}$ and $M_{\text{disc}} = 0.05 M_{\odot}$ provide a decent agreement with the data.

## 3.3 Planetary migration: starting point at 86 au

In Section 3.2, we discarded planet migration, assuming that the planet orbit was fixed. However, given that some of our models feature a rather large disc mass, one may worry that a substantial
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Time evolution of the planet's orbital radius in the simulations presented in Section 3.3, in which the planet starts migrating at 86 au after having remained on a fixed orbit over 100 orbital periods. Different curves correspond to different disc masses, the green curve (0.05 M☉) corresponds to the disc mass used in the best model described in Section 3.2.</figcaption></figure>

migration of the planet could occur over the duration of the simulations. In addition, the planet ability to carve an annular gap in the disc, the dust dynamics and its emission could all be affected by planetary migration as well (e.g. Meru et al. 2019; Wafflard-Fernandez & Baruteau 2020). Therefore, in this section we relax the fixed orbit assumption, letting the planet interact with the disc and hence migrate freely as directed by the gravitational torque from the disc (cf. Section 2.1.3).

For a migrating planet, its starting position may be defined as the orbital radius where the planet has reached its current mass. This radius is, however, very uncertain as the onset of runaway gas accretion and its pace depend on the physical properties of the disc in the planet's circumplanetary environment, in particular on its cooling rate, and on its heating rate via solid phase accretion (e.g. Ayliffe & Bate 2009; Szulágyi et al. 2014; Emsenhuber et al. 2021). A variety of recent observational data and simulations actually indicate that planets may grow much slower than usually believed (e.g. see the references and discussion in Nayakshin, Elbakyan & Rosotti 2022). Therefore, we have adopted two different strategies about the starting position of the planet ($R_{p,start}$): either it is released at 86 au, the same location as in Section 3.2, which is what we assume in the current section, or it is released at an arbitrarily larger separation of 130 au, which is what will be assumed from the next section onward. We also define the planet's migration time-scale, $t_m$, as

$$t_m = t_{\text{sim}} \frac{R_{p,start}}{R_{p,start} - R_{p,end}}, \qquad (7)$$

where $R_{p,end}$ denotes the planet's orbital radius reached at the end of the simulation, $t_{sim}$. Table 1 lists $t_m$ obtained in our simulations.

Since the rate of planetary migration depends sensitively on the disc's surface density, we have explored a range of disc masses from a minimum of $M_{disc} = 0.0125 M_{☉}$ to a maximum value of $M_{disc} = 0.2 M_{☉}$ in steps by a factor of 2. We display the time evolution of the planet's orbital radius in Fig. 2. The maximum duration of our simulations, $t_{sim,max}$, is approximately 630 kyr, which is about 1130 orbital periods at the planet's initial separation of 86 au. Some simulations are terminated earlier, however, if the planet reaches the inner boundary earlier than 630 kyr.

Fig. 2 shows that, not surprisingly, the planet migration rate increases with disc mass (see e.g. the review by Baruteau et al. 2014). So long as $M_{disc} \le 0.025 M_{☉}$, our planet migrates inward

smoothly in a near type I migration fashion (the planet produces just a weak depression in gas around its orbit, as can be appreciated from the gas surface densities shown in Fig. 4). For $M_{disc} \ge 0.05 M_{☉}$, however, our planet migrates inward at a highly variable rate, with several stages of fast, accelerated type III (runaway) migration, sandwiched between stages of slower migration. Intermittent runaway migration has already been observed and documented in several studies (e.g. Lin & Papaloizou 2010; McNally et al. 2019; Wafflard-Fernandez & Baruteau 2020). The intermittency of inward runaway migration is intimately related to the time evolution of the planet's vorticity-weighted coorbital mass deficit (for definition see Masset & Papaloizou 2003), which can fluctuate around the planet mass depending on the surface density profile of the background gas, and how strongly the inner wake of the planet shocks the background gas (Wafflard-Fernandez & Baruteau 2020). As can be seen with the green curve in Fig. 2, fast, intermittent inward runaway migration is obtained indeed for $M_{disc} = 0.05 M_{☉}$, which was the disc mass selected for our 'best-fit' model in Section 3.2.

As argued by Nayakshin et al. (2022), very rapid planet migration is a challenge to the widely accepted scenario in which all of the gaps/rings seen by ALMA are due to embedded planets (e.g. Long et al. 2018). In regards to HD 163 296 in particular, we find that for $M_{disc} \ge 0.05 M_{☉}$ the migration time as defined by equation (7) is shorter than ~0.4 Myr. This is very short compared with the estimated age of this system ($t_* \sim 7$ to 10 Myr). Another way to look at this is to define the gap crossing time $t_{gap}$, which is the time the planet spends while migrating over the full radial width of its annular gap. This time is shorter than $t_m$ by approximately $\Delta R/R \approx 5 \times (q/3)^{1/3}$, since the typical radial width of a planetary gap (in the gas) carved by a giant planet is about 5 Hill radii (Masset, D'Angelo & Kley 2006). Thus, to stand a decent chance of observing a giant planet passing through the 86 au gap and also opening it, HD 163 296 needs to be hatching gas giant planets beyond 86 au at a rate of approximately $t_{gap}^{-1}$. One can then estimate the number of gas giant planets that HD 163 296 needs to produce over its lifetime as

$$N_p \sim \frac{t_*}{t_{gap}} \sim 100 \qquad (8)$$

for $q = 1.3 \times 10^{-4}$, $t_* = 7$ Myr, and $t_m = 0.4$ Myr. This is probably excessive. More generally, it is unlikely that a planet is in the gap at the time of observation if it is migrating rapidly, unless there are $N_p$ planets passing through the gap in succession. On these grounds we conclude that planet migration constrains the disc mass in HD 163 296 to be no larger than ~0.025 $M_{☉}$.

## 3.4 Planetary migration: starting point at 130 au

Previous Section 3.3 shows that there ought to be an upper limit for the mass of the disc around HD 163 296 to avoid fast inward migration of the putative planet responsible for the dark ring of continuum emission near 86 au. The simulations presented in this section assume that this putative planet formed and reached its final mass further out in the disc, specifically at 130 au. The planet is introduced by increasing its mass to its final value over 10 local orbital periods during which the planet remains on a fixed orbit. Simulations then let the planet migrate in the disc, and are terminated when the planet reaches 86 au. In this section, we explore two disc masses: 0.025 $M_{☉}$ (Model PM025) and 0.05 $M_{☉}$ (Model PM05). For both disc masses, the migration from 130 au down to 86 au proceeded smoothly, without runaway, and lasted for about 0.68 Myr and 0.30 Myr for models PM025 and PM05, respectively. For completeness, we have checked by restarting both of these
<figure><img src="image_3.png" /><figcaption>Figure 3. Similar to Fig. 1's synthetic maps but for Model PM025 (left) and Model PM05 (right).</figcaption></figure>

simulations that, below 86 au, the planet qualitatively has a similar orbital evolution as in the simulations where the planet was inserted at 86 au. In particular, the planet in Model PM05 undergoes a stage of accelerated type III migration between ~65 and ~45 au, much like in Model OM05. This indicates that the exact starting position of the planet (130 au here) is not very important so long as one runs the simulations long enough so that the planet eventually reaches its current 'observed' position at 86 au.

Fig. 3 shows the synthetic map of continuum emission for both models when the planet reaches about 86 au. The same colour scale as in Fig. 1 is used, to provide a better comparison with the observed map. For Model PM025, shown on the left-hand side, we no longer see a bright ring of emission outward of the planet gap compared to the fixed planet model at same disc mass (Model NM025) displayed in the lower-right panel of Fig. 1. We only see a dim emission inward of the planet gap at a similar level of flux as in Model NM025. A similar comparison can be made between Model PM05, shown on the right-hand side of Fig. 3, and Model NM05, shown in the top-right panel in Fig. 1: the outer bright ring is much dimmer with Model PM05, and we also notice a secondary dark ring inward of the planet gap at an orbital separation ~60 au. Since the dust temperature is the same in all our synthetic maps (recall that it is simply the gas temperature in the hydrodynamical simulations), the overall flux differences between Models PM025 and PM05 therefore indicate that the dust's continuum emission is mostly optically thin for both models.

To get more insight into the differences in continuum emission between the fixed and migrating planet models, we display in Fig. 4 the azimuthally averaged surface density profiles of the gas and of each dust fluid for Models PM025 (left-hand panels) and PM05 (right-hand panels). Profiles are shown at two epochs: at 250 orbital periods at 100 au (0.17 Myr) in the upper panels, and when the planet has reached ~86 au. The panels illustrate how dust radial drift and dust trapping about the planet vary with dust size and disc mass.

For Model PM025, we see that both 1 and 3 mm dust fluids drift inward faster than the planet migrates in. This, and the smaller dust

mass, explain why the continuum emission at 1.3 mm from the inner disc is dim in this model (as already discussed in Section 3.2).

Concerning the outer bright ring, the fact that the planet migrates in Model PM025 relatively rapidly implies that it cannot trap the largest dust outside its annular gap as vigorously as in Model NM025. This explains the lack of a strong bright ring outward of the planet gap in Model PM025 compared with NM025.

Now, for model PM05, only the 3 mm dust fluid tends to drift inward faster than the planet migrates in, while the 1 mm dust drifts at a pace that is no faster than the planet's migration rate. The presence of the 1 mm dust in the inner disc regions as well as beyond the planet gap explains the different flux distributions of Models PM025 and PM05 (Meru et al. 2019; Nazari et al. 2019). Still, the lack of 3 mm dust trapped at the outer edge of the planet gap in Model PM05 explains why the outer bright ring in PM05 is dimmer compared to Model NM05, where the fixed planetary orbit helps to build up a stronger pressure maximum outside of the planet location. This can be further appreciated by inspection of the dust surface density profiles for Model NM05, which are shown in the left-hand panel of Fig. 5. We finally point out that the dark ring of emission just inward of the planet location in Model PM05 is mainly related to the local minimum in the dust surface density of the 3 mm dust fluid, which can be seen around 60 au in the bottom-right panel of Fig. 4.

## 3.5 Dust material density: compact or porous?

Previous section shows that the smaller the mass of the disc in HD 163296, the longer the migration time-scale for the putative planet at 86 au. This relieves the statistical requirement for HD 163 296 disc to hatch dozens of planets (cf. the end of Section 3.3) but then we encounter the dust drift challenge. The smaller the disc mass, the faster is dust radial drift, and hence the smaller the flux of continuum emission compared to observations. We find ourselves at an impasse. At high disc mass, the planet is lost into the star too quickly, at low disc masses the dust drifts in too rapidly; and there appears no
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Surface densities (in $g \text{ cm}^{-2}$) of gas and dust fluids in Model PM025 (left-hand panels) and Model PM05 (right-hand panels) resulting from our hydrodynamical simulations at two different times: at 250 orbits (upper panels) and when the planet reaches 86 au (lower panels). Different colours correspond to different dust sizes (see legend in each figure). The black cross symbol marks the position of the planet at each time in the simulation. The vertical dotted lines show the initial (130 au) and final (86 au) positions of the planet in the simulations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Similar to Fig. 4 but for Model NM05 (left) and Model PM025p (right). For NM05, surface density profiles are shown at the same time as in the synthetic image in the top-right panel in Fig. 1. For PM025p, the time is that when the planet reaches 86 au.</figcaption></figure>

intermediate value of disc mass that would resolve these opposite challenges.

To avoid such a fast radial drift, we explore in this section the possibility that dust particles have a smaller material density than conventionally assumed. More specifically, we now set the internal density of dust particles to $\rho_d = 0.1 \text{ g cm}^{-3}$, like in Baruteau et al. (2019), and we refer to them as porous dust particles (as opposed to compact dust particles for our default internal density). The basic

idea is that porous dust particles will have lower Stokes numbers than compact particles at a given disc mass (equation 4), and therefore a slower radial drift. Their absorption opacity will differ too, and dust radiative transfer calculations are needed to check how the overall level of flux emission is impacted.

The right-hand panel in Fig. 5 shows the gas and dust surface density profiles for simulation PM025p, for which the disc mass is $0.025 \text{ M}_\odot$ and porous dust particles are adopted. The density profiles
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Similar to Fig. 3 but for Model PM025p (left) and Model PM025d10p (right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Planet migration and dust drift time-scales versus disc mass. In particular, $t_{\text{mig1}}$ is the expected type I (equation 9) and $t_m$ is numerically determined migration time-scales (equation 7 and Table 1). The later is shorter due to the presence of runaway migration. The dust drift time-scale $t_{\text{dr}}$ are given by equation (10) for 1 mm dust particles for either compact (bottom blue line) or porous (top green) dust models. For compact dust there is no acceptable solution for the disc mass: at low disc masses the dust drifts into the star too quickly, whereas at high disc masses the planet is lost into the star by runaway migration. Porous dust model largely alleviates these issues (see text in Section 4).</figcaption></figure>

are those when the planet reaches about 86 au, and they can be directly compared to those obtained in the previous section for compact particles (bottom-left panel in Fig. 4). A comparison between both panels shows that, as expected, dust radial drift is much slower for porous particles. At the time when the planet reaches 86 au, the surface density of 1 and 3 mm dust particles is high and close to its initial value. On the other hand, the synthetic image displayed in the left-hand panel in Fig. 6, shows a very low level of continuum flux, which arises from porous dust particles having a smaller absorption opacity compared to compact particles of same size (Kataoka et al. 2014).

We thus further tested increasing the dust-to-gas mass ratio $\epsilon$. As shown in the right-hand panel in Fig. 6, increasing $\epsilon$ from 0.01 to 0.1 gives a satisfactory agreement between the properties of the synthetic map (overall level of flux, gap width etc.) and those of the ALMA observation.

## 4 CONCLUDING REMARKS

There is mounting observational evidence that annular sub-structures are very common features of discs continuum emission. Interestingly, not only are bright and dark rings frequently observed, like amongst the Taurus discs (Long et al. 2018), they also tend to come as multiple-ring systems (Huang et al. 2018). This is the case of HD 163 296 disc, for which the detection of localized velocity perturbations in the CO emission strongly hints at the presence of massive planets in (some of) its dark rings (Teague et al. 2018; Pinte et al. 2020). The wealth of observational data on the HD 163 296 disc clearly makes it an ideal target to test predictions of disc-planet interactions. Though many numerical works have intended to reproduce observed annular structures with planets (a statement that is not restricted to HD 163 296 disc), only a handful of them, including this work, have taken planetary migration into account.

Taking the HD 163 296 disc as a case study, here we have shown that planetary migration could be used to set an upper limit on the mass of the disc. For the HD 163 296 disc this upper limit is ~0.025 $M_{\odot}$, which is at the low end of values compatible with the gas or dust emission, this disc mass is also lower than previous recommended values (e.g. ~0.15 $M_{\odot}$, see Zhang et al. 2021). While a lower disc mass would result in slower migration, and therefore to a higher probability of catching a migrating planet sculpting its disc, it would also lead to a faster dust radial drift. As a result we found no value for a disc mass that is compatible with both planet migration and dust emission constraints. The situation is illustrated in Fig. 7.

The gist of the argument can be glipsed from the following semi-analytical argument. The theory of type I migration in isothermal discs leads to migration time-scales depending on disc properties as
$t_{\text{migl}} = C h^2 M_*^2 / (\Omega_K M_p \Sigma R^2)$, with dimensionless coefficient C depending on local disc properties (e.g. Tanaka, Takeuchi & Ward 2002; Paardekooper et al. 2010). The disc surface density is proportional to the total disc mass $M_d$ through equation (1), and therefore we see that this scales as $t_{\text{migl}} \propto 1/M_d$. Using the results of our lower mass OM125 and OM025 simulations that showed no runaway migration to deduce the constant C from simulations directly, we get

$$t_{\text{migl}} \sim 2 \times 10^6 \text{yr} \frac{0.025 M_{\odot}}{M_d} \frac{0.26 M_J}{M_p}. \quad (9)$$

This time-scale is short but not overly so compared with the age of HD 163296. Our hydrodynamical simulations actually place a surprisingly robust upper limit on $M_d$ that is much more stringent than can be expected based on type I migration theory (equation 9). We find indeed that at disc masses higher than ~0.025 $M_{\odot}$ runaway migration (e.g. Baruteau et al. 2014) shortens planet migration time to a fraction of a Myr. Since the gap crossing time is yet shorter than planet migration time-scale (see the end of Section 3.3), we then estimated that up to 100 planets need to be injected into the disc of this system to present a decent statistical chance of us seeing the planet.

At the same time, the dust particles drift time-scale (e.g. Armitage 2007) is also a function of disc mass, albeit one that increases with $M_d$ in the limit of small Stokes numbers:

$$t_{\text{dr}} = \frac{1}{\Omega_K} \frac{\text{St}^{-1}}{h^2 \eta} \sim 10^5 \text{ years} \frac{M_d}{0.025 M_{\odot}} \frac{1 \text{ mm}}{s} \frac{1.3 \text{ g cm}^{-3}}{\rho_d}, \quad (10)$$

where $\eta = -d\ln(P)/d\ln R \approx 3.2$ with $\langle P \rangle$ the azimuthally averaged disc pressure (this estimate is made at $R = 86$ au).

Comparing equations (9) and (10) we see that there is no value of disc mass that allows both gas and dust to be present in the disc of HD 163296 for a time comparable with the 7–10 Myr old age of this system. This is clearly a significant challenge to our understanding of how this system evolved to be in the current state.

We find that a porous rather than compact dust composition presents a potential solution to this dust drift versus planet migration tension. Simulations with strongly decreased dust density, $\rho_d = 0.1 \text{ g cm}^{-3}$, have shown that the dust can in fact stay in the disc for far longer. This is natural given equation (10). To account for the observed level of dust continuum flux in HD 163296, the disc needs to be rather dust-rich, with a dust-to-gas mass ratio about 5 to 10 times larger than the canonical value of 0.01 used in the literature. However, the total dust mass that we infer is about 0.002 $M_{\odot}$, which is actually in good agreement with that inferred from the dust continuum emission (Powell et al. 2019). Such a high dust-to-gas mass ratio in an old disc would imply that there were significantly more gas in the disc at earlier times, and that most of it, say ~90%, has since been removed from the disc. Such a process must selectively remove gas but not dust.

Gas could be lost from the disc via gas accretion onto the star. However, as we need dust particles to be small in order to minimize dust drift into the star, we would expect the dust to be carried into the star with the gas. It is hence difficult to see why gas accretion onto the star would selectively remove gas from the disc. On the other hand, there have been suggestions that dust in older protoplanetary discs may be of secondary nature (e.g. Turrini et al. 2019). In this scenario primordial dust is locked into planetesimals relatively early on. Once giant planets form, planetesimals are excited onto eccentric orbits. Like in debris discs (e.g. Wyatt 2008), planetesimal collisions then give rise to fragmentation cascade leading to dust production. Secondary dust has in fact been proposed for HD 163296 (D'Angelo & Marzari 2022). Gas outflows from the disc can

also remove gas preferentially, leaving larger dust particles behind (Hutchison et al. 2016).

As a final word, we recall that, for the sake of quantitative research, we have only studied the case of a single planet embedded in a protoplanetary disc with physical parameters inspired from the HD 163296 disc. Considering multiple planets could certainly impact their orbital evolution but it may not go against their natural tendency to undergo fast runaway migration if the disc is massive enough. As an example, we mention again the recent study of Garrido-Deutelmoser et al. (2023), who have shown that the many substructures in the dust continuum emission at 1.25 mm wavelenth, including the intriguing comma-shaped feature inside the dark ring of emission near 54 au, could be reproduced by invoking four planets in a resonant chain. However, their model did not take into account planetary migration while the disc mass was taken to be 0.15$M_{\odot}$. Should planet migration be accounted for in their simulation, planets would not end up in a resonant chain but be destabilized instead by fast encounters due to disc-induced runaway migration. This is shown in Section B in the Appendix.

## ACKNOWLEDGEMENTS

We thank Juan Garrido-Deutelmoser for sharing his problem setting and the referee for helpful comments. We thank Yi-Xian Chen, Haochang Jiang, Kaitlin Kratter, Adrien Leleu and Douglas N. C. Lin for constructive suggestions. YW further thanks Wenrui Xu, Shuo Huang, Shangjia Zhang, and Kiyoki Doi for helpful discussions. This research used the ALICE High Performance Computing Facility at the University of Leicester, and DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). Part of the numerical simulations were performed on the CALMIP Supercomputing Centre of the University of Toulouse. YW gratefully acknowledges the support by the DUSTBUSTERS RISE project (grant agreement number 823823) for his secondment in University of Arizona. SN acknowledges the funding from the UK Science and Technologies Facilities Council, grant No. ST/S000453/1.

## DATA AVAILABILITY

The FARGO3D code is publicly available from http://fargo.in2p3.fr. The RADMC-3D code is available from https://www.itau.uni-heidelberg.de/~dullemond/software/radmc-3d. The python package used for the post-processing of the simulations data and for the analysis of the radiative transfer calculations are available via its GitHub repository: fargo2radmc3d (https://github.com/charango/fargo2radmc3d). The data obtained in our simulations can be made available on reasonable request to the corresponding author.

## REFERENCES

* ALMA Partnership, 2015, ApJ, 808, L3

* Andrews S. M., Terrell M., Tripathi A., Ansdell M., Williams J. P., Wilner D. J., 2018a, ApJ, 865, 157

* Andrews S. M. et al., 2018b, ApJ, 869, L41

* Armitage P. J., 2007, ApJ, 665, 1381

* Ayliffe B. A., Bate M. R., 2009, MNRAS, 397, 657

* Baruteau C., Masset F., 2008, ApJ, 678, 483

* Baruteau C. et al., 2014, Protostars and Planets VI, University of Arizona Press, Tucson, p. 667

* Baruteau C. et al., 2019, MNRAS, 486, 304

* Baruteau C., Wafflard-Fernandez G., Le Gal R., Debras F., Carmona A., Fuente A., Rivière-Marichalar P., 2021, MNRAS, 505, 359
Benítez-Llambay P., Masset F. S., 2016, ApJS, 223, 11

Benítez-Llambay P., Krapp L., Pessah M. E., 2019, ApJS, 241, 25

Bergin E. A. et al., 2013, Nature, 493, 644

Bergin E. A., Du F., Cleeves L. I., Blake G. A., Schwarz K., Visser R., Zhang K., 2016, ApJ, 831, 101

Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148

Booth A. S., Walsh C., Ilee J. D., Notsu S., Qi C., Nomura H., Akiyama E., 2019, ApJ, 882, L31

Carmona A., 2010, Earth Moon Planets, 106, 71

Casassus S., Pérez S., 2019, ApJ, 883, L41

Crida A., 2009, ApJ, 698, 606

Crida A., Baruteau C., Kley W., Masset F., 2009, A&A, 502, 679

D'Angelo G., Marzari F., 2022, MNRAS, 509, 3181

Dong R., Li S., Chiang E., Li H., 2018, ApJ, 866, 110

Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, Astrophysics Source Code Library, recordascl:1202.015

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

Elbakyan V., Wu Y., Nayakshin S., Rosotti G., 2022, MNRAS, 515, 3113

Emsenhuber A., Mordasini C., Burn R., Alibert Y., Benz W., Asphaug E., 2021, A&A, 656, A70

Flaherty K. M. et al., 2017, ApJ, 843, 150

Franceschi R., Birnstiel T., Henning T., Pinilla P., Semenov D., Zormpas A., 2022, A&A, 657, A74

Garrido-Deutelmoser J., Petrovich C., Charalambous C., Guzmán V. V., Zhang K., 2023, ApJ, 945, L37

Huang S., Ormel C. W., 2023, MNRAS, 522, 828

Huang J. et al., 2018, ApJ, 869, L42

Hutchison M. A., Price D. J., Laibe G., Maddison S. T., 2016, MNRAS, 461, 742

Isella A. et al., 2018, ApJ, 869, L49

Jones M. G., Pringle J. E., Alexander R. D., 2012, MNRAS, 419, 925

Kama M. et al., 2020, A&A, 634, A88

Kamp I., Thi W. F., Woitke P., Rab C., Bouma S., Ménard F., 2017, A&A, 607, A41

Kanagawa K. D., Muto T., Tanaka H., Tanigawa T., Takeuchi T., Tsukagoshi T., Momose M., 2016, PASJ, 68, 43

Kataoka A., Okuzumi S., Tanaka H., Nomura H., 2014, A&A, 568, A42

Lesur G. R. J., 2021, A&A, 650, A35

Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874

Lin M.-K., Papaloizou J. C. B., 2010, MNRAS, 405, 1473

Liu S.-F., Jin S., Li S., Isella A., Li H., 2018, ApJ, 857, 87

Liu Y. et al., 2022, Science China Phys. Mech. Astron., 65, 129511

Long F. et al., 2018, ApJ, 869, 17

Masset F. S., Papaloizou J. C. B., 2003, ApJ, 588, 494

Masset F. S., D'Angelo G., Kley W., 2006, ApJ, 652, 730

McNally C. P., Nelson R. P., Paardekooper S.-J., Benítez-Llambay P., 2019, MNRAS, 484, 728

Meru F., Juhasz A., Ilee J. D., Clarke C. J., Rosotti G. P., Booth R. A., 2017, ApJ, 839, L24

Meru F., Rosotti G. P., Booth R. A., Nazari P., Clarke C. J., 2019, MNRAS, 482, 3678

Miotello A., Kamp I., Birnstiel T., Cleeves L. I., Kataoka A., 2022, preprint (arXiv:2203.09818)

Müller T. W. A., Kley W., Meru F., 2012, A&A, 541, A123

Nayakshin S., Elbakyan V., Rosotti G., 2022, MNRAS, 512, 6038

Nazari P., Booth R. A., Clarke C. J., Rosotti G. P., Tazzari M., Juhasz A., Meru F., 2019, MNRAS, 485, 5914

Paardekooper S. J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950

Penzlin A. B. T., Ataiee S., Kley W., 2019, A&A, 630, L1

Pérez S., Casassus S., Baruteau C., Dong R., Hales A., Cieza L., 2019, AJ, 158, 15

Pinte C. et al., 2020, ApJ, 890, L9

Powell D., Murray-Clay R., Schlichting H. E., 2017, ApJ, 840, 93

Powell D., Murray-Clay R., Pérez L. M., Schlichting H. E., Rosenthal M., 2019, ApJ, 878, 116

Riols A., Lesur G., Menard F., 2020, A&A, 639, A95

Schwarz K. R. et al., 2021, ApJS, 257, 20

Setterholm B. R. et al., 2018, ApJ, 869, 164

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

Suzuki T. K., Ogihara M., Morbidelli A., Crida A., Guillot T., 2016, A&A, 596, A74

Szulágyi J., Morbidelli A., Crida A., Masset F., 2014, ApJ, 782, 65

Tabone B., Rosotti G. P., Cridland A. J., Armitage P. J., Lodato G., 2023, MNRAS, 519, 2917

Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257

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

Trapman L., Miotello A., Kama M., van Dishoeck E. F., Bruderer S., 2017, A&A, 605, A69

Tsukagoshi T. et al., 2019, ApJ, 878, L8

Turrini D., Marzari F., Polychroni D., Testi L., 2019, ApJ, 877, 50

Veronesi B., Paneque-Carreño T., Lodato G., Testi L., Pérez L. M., Bertin G., Hall C., 2021, ApJ, 914, L27

Vioque M., Oudmaijer R. D., Baines D., Mendigutía I., Pérez-Martínez R., 2018, A&A, 620, A128

Wafflard-Fernandez G., Baruteau C., 2020, MNRAS, 493, 5892

Weber P., Pérez S., Benítez-Llambay P., Gressel O., Casassus S., Krapp L., 2019, ApJ, 884, 178

Woitke P. et al., 2016, A&A, 586, A103

Wu Y., Chen Y.-X., Jiang H., Dong R., Macías E., Lin M.-K., Rosotti G. P., Elbakyan V., 2023, MNRAS, 523, 2630

Wyatt M. C., 2008, ARA&A, 46, 339

Xin Z., Espaillat C. C., Rilinger A. M., Ribas A., Macias E., 2023, ApJ, 942, 4

Zhang S. et al., 2018, ApJ, 869, L47

Zhang K. et al., 2021, ApJS, 257, 5

Zhu Z. et al., 2019, ApJ, 877, L18

## APPENDIX A: HIGHER AND LOWER PLANET MASS

As described in Section 2.1, we use the results of Dong et al. (2018) – their equation (10) – to derive the mass of the planet on a fixed circular orbit that can reproduce the radial half-width of the gap around 86 au in the HD 163296 disc (this yields $M_p \approx 0.26 M_J$ for our disc’s aspect ratio). However, as a part of the parameter space study, we also carried out several simulations with $M_p$ twice higher or lower than that. The simulations are shown in Table 1 as NM05H (higher planet mass) for $M_p = 0.52 M_J$, and NM05L (lower planet mass) for $M_p = 0.13 M_J$. Like Model NM05, these additional models do not allow for planet migration, but we also ran the same two models with planet migration allowed (these are labelled as Models PM05H and PM05L).

Fig. A1 shows synthetic images of the dust continuum emission for these additional models. For model NM05H (fixed planet, higher mass), the dark ring of continuum emission is significantly wider compared to our fiducial Model NM05 (compare with the top-right panel in Fig. 1). For model NM05L (fixed planet, lower mass), the width of the dark ring is smaller than for Model NM05 and actually in pretty good agreement with the observed one. The width and flux level around the outer bright ring agree somewhat less with the observation, although this could be improved by slightly increasing the dust-to-gas mass ratio. This actually illustrates the degeneracy in the modelling of observed disc structures. Regarding the additional models with migration, the one with higher planet mass (PM05H) undergoes rather fast migration, which causes a very wide dust gap (it takes only 260 orbits for the planet to move from 130 to 86 au). And finally the model with lower mass (PM05L) shows no dust trapping at the gap’s outer edge and therefore no outer bright ring (the larger
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Synthetic maps of continuum emission for additional models with different planet masses. In the left column of panels (Models NM05H and NM05L), the planet is held on a fixed circular orbit at 86 au, and the synthetic images are shown at the same time as those in Fig. 1). In the right column of panels, the planet migrates from 130 au and images are shown when it reaches about 86 au (which is at *t* = 260 orbits for Model PM05H and *t* = 620 orbits for Model PM05L).</figcaption></figure>

dust drifts faster than the planet and mostly end up in the inner disc). This further shows how planet migration and dust drift, when both included in the simulations, are able to reduce degeneracy of the models.

## APPENDIX B: A MULTIPLE PLANET SCENARIO IN HD 163296

Although the HD 163296 disc likely features multiple planets, our model, which is inspired from this disc, only has one planet for

simplicity. To test how our conclusions would be changed by the presence of several planets, we carried out an additional simulation with multiple migrating planets. More specifically, we used the exact same disc and planet parameters as Garrido-Deutelmoser et al. (2023), where four near-resonant planets were simulated on fixed orbits. These four planets are located at 46, 54, 84.5, and 137 au with mass 0.255 $M_J$, 0.18 $M_J$, 0.4 $M_J$, and 1 $M_J$, respectively. We modified our source files to adopt the same radial profile of $\alpha$ turbulent viscosity as in equation (2) of Garrido-Deutelmoser et al. (2023). The disc mass is 0.15 $M_\odot$.
Just like Garrido-Deutelmoser et al. (2023), we first ran this simulation for 2000 orbits at 48 au without letting the planets migrate. We have checked that the gas density that we obtain agrees very well with their fig. 1. We note that this requires the gas indirect term to be discarded in the simulation, otherwise the disc gas grows an unstable *m* = 1 mode that completely destabilizes the disc after a few hundred orbits. We then restarted the simulation by allowing the planets to migrate (planets interact with the disc other and with each other).

<figure><img src="image_9.png" /><figcaption>**Figure B1.** Similar to Fig. 2 but for the four-planet model presented in Section B of the Appendix.</figcaption></figure>

The time evolution of the planets orbital radius is displayed in Fig. B1. It is clear that the planets do not maintain their initial (near-resonant) orbits, perhaps not too surprisingly given the high disc mass (much higher than the one from which we get runaway migration with a single planet). Although planet-planet interactions are strong, disc-induced runaway migration remains inevitable, more particularly for the innermost and outermost planets. We point out that the second and third planets temporarily share the same orbital radius for a few hundred orbits while being ~180 degrees far apart. To our knowledge only a few hydrodynamical simulations of disc-planets interactions have reported that two planets could reach a 1:1 mean-motion resonance through convergent migration; in particular, in the work of Penzlin, Ataiee & Kley (2019) for eccentric circumbinary planets, and that of Crida (2009), where a 1:1 resonance is obtained between Uranus and Neptune in a simulation of the four giants in the Solar system that involves a massive version of the Minimum Mass Solar Nebula.

If indeed the HD 163296 disc has several planets near mean-motion resonances, it seems unlikely that the planets acquired such near-resonant states with a high disc mass (Huang & Ormel 2023). All in all, whether a single planet or multiple planets are invoked, our main conclusions regarding the impact of disc mass on disc migration remain unaltered.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.