A Systematic Review on Anomaly Based Intrusion Detection System

Astronomy & Astrophysics manuscript no. main December 18, 2020 ©ESO 2020 Modeling the nonaxisymmetric structure in the HD 163296 disk with planet-disk interaction P. J. Rodenkirch1 , Thomas Rometsch2 , C. P. Dullemond1 , Philipp Weber345 , and Wilhelm Kley2 1 2 3 arXiv:2012.09217v1 [astro-ph.EP] 16 Dec 2020 4 5 Institute for Theoretical Astrophysics, Zentrum für Astronomie, Heidelberg University, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany Niels Bohr International Academy, The Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen Ø, Denmark Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile Departamento de Física, Universidad de Santiago de Chile, Av. Ecuador 3493, Estación Central, Santiago, Chile December 18, 2020 ABSTRACT Context. High resolution ALMA observations like the DSHARP campaign revealed a variety of rich substructures in numerous protoplanetary disks. These structures consist of rings, gaps and asymmetric features. It is debated whether planets can be accounted for these substructures in the dust continuum. Characterizing the origin of asymmetries as seen in HD 163296 might lead to a better understanding of planet formation and the underlying physical parameters of the system. Aims. We test the possibility of the formation of the crescent-shaped asymmetry in the HD 163296 disk through planet disk interaction. The goal is to obtain constraints on planet masses and eccentricities and disk viscosities.We furthermore test the reproducibility of the two prominent rings in the HD 163296 disk at 67 au and 100 au. Methods. Two dimensional, multi-fluid, hydrodynamical simulations are performed with the FARGO3D code including three embedded planets. Dust is described with the pressureless fluid approach and is distributed over eight size bins. Resulting grids are post-processed with the radiative transfer code RADMC-3D and the CASA software to model synthetic observations. Results. We find that the crescent-shaped asymmetry can be qualitatively modeled with a Jupiter mass planet at a radial distance of 48 au. Dust is trapped preferably in the trailing Lagrange point L5 with a mass of 10 to 15 earth masses. The observation of such a feature confines the level of viscosity and planetary mass. Increased values of eccentricity of the innermost Jupiter mass planet damages the stability of the crescent-shaped feature and does not reproduce the observed radial proximity to the first prominent ring in the system. Generally, a low level of viscosity (α ≤ 2 · 10−3 ) is necessary to allow the existence of such a feature. Including dust feedback the leading point L4 can dominantly capture dust for dust grains with an initial Stokes number ≤ 3.6 · 10−2 . In the synthetic ALMA observation of the model with dust feedback two crescent-shaped features are visible. The observational results suggest a negligible effect of dust feedback since only one such feature has been detected so far. The dust-to-gas ratio may thus be overestimated in the models. Additionally, the planet mass growth time scale does not strongly affect the formation of such asymmetries in the co-orbital region. Key words. protoplanetary disks - planet-disk interactions - planets and satellites: formation - planets and satellites: rings - hydro- dynamics - radiative transfer 1. Introduction In the advent of high angular resolution millimeter continuum observations with the Atacama Large Millimeter and Submillimeter Array (ALMA) insights into substructures of protoplanetary disks have become available. The striking results of the first highly resolved observation of the protoplanetary disk HL Tau unveiled a rich variety of concentric rings in the dust continuum (ALMA Partnership et al., 2015). An extensive survey with the goal to image detailed structures in 20 disks was performed with the DSHARP campaign (Andrews et al., 2018) with unprecedented resolution. Rings and gaps seem to be ubiquitous in these disks and appear independently of the stellar luminosity (Huang et al., 2018). A subset of the observations display nonaxisymmetric structures like spirals and crescent-shaped features. Currently, it is debated whether these structures are signposts of embedded planets (Zhang et al., 2018). Planet-disk interaction has been a central topic in the dynamics of protoplanetary disks, first with analytic studies of resonances and spiral density waves (Goldreich & Tremaine, 1979, 1980) or planetary migration (Lin & Papaloizou, 1986). Jupiter mass planets can open a gap in the gas (Kley, 1999). As shown by numerical, two fluid simulations by Paardekooper & Mellema (2004) a planet of 0.1 Jupiter masses (Mjup ) is sufficient to open a gap in the dust. Even lower masses down to 0.05 Mjup can lead to gap formation if mainly mm-sized dust particles are present (Paardekooper & Mellema, 2006). Dust structures created by planet-disk interaction are generally more diverse than their counterpart in the gas (Fouchet et al., 2007; Maddison et al., 2007). For massive planets of 5 Mjup and cm-sized grains Fouchet et al. (2010) found azimuthally asymmetric dust trapping in the context of 3D SPH simulations. Embedded dust grains are prone to drift towards pressure maxima in the disk (Whipple, 1972). Thus, perturbations caused by a sufficiently massive planet can efficiently trap up to meter-sized bodies on the outer edge of the gap, form a ring structure and may aid planArticle number, page 1 of 19 A&A proofs: manuscript no. main etesimal formation (Ayliffe et al., 2012). Dong et al. (2015) found that multiple planets can explain large cavities at near-infrared and millimeter wavelengths as observed in transition disks (Calvet et al., 2005; Hughes et al., 2009). In a system with two planets dust trapped in the leading and trailing Lagrange points (L4 & L5) can be a transient feature, depending on the outer planet (Picogna & Kley, 2015). In the context of low viscosity disks multiple rings and gaps emerge with a single planet through shocks of the primary and secondary spiral arm (Zhu et al., 2014; Bae et al., 2017). Torques caused by the gravitational interaction between planets and the disk lead to migration (Kley & Nelson [2012] for a review) which in turn affects the observable dust substructures, e.g. changes in the ring intensity or asymmetric triple ring structures depending on the migration rate and direction (Meru et al., 2019; Weber et al., 2019). Migration is sensitive to the underlying disk physics and can be chaotic in very low viscosity disks (McNally et al., 2019). In general, protoplanetary disks seem to be only weakly turbulent, indicating regimes of α on the order of 10−4 to 10−3 (Flaherty et al., 2015, 2017; Dullemond et al., 2018), using the turbulent α viscosity parametrization from Shakura & Sunyaev (1973). When referring to low viscosity in disks, the magnitude of the effective α viscosity is meant which drives angular momentum transport and thus accretion due the underlying turbulent processes. Hence, a low effective viscosity can be linked to weak turbulence. Spiral waves in the gas, excited by a planet, are mostly hidden in the dust dynamics, favoring gaps and rings (Dipierro et al., 2015). The gravitational instability (Toomre, 1964) in sufficiently massive disks can however trigger spiral waves trapping large particles (Rice et al., 2004, 2006). These waves are in principle also observable in scattered light observations (Pohl et al., 2015). Nonaxisymmetric features like vortices can be created by the Rossby wave instability (Lovelace et al., 1999; Li et al., 2000) enabling dust trapping (Baruteau & Zhu, 2016). Observationally these might be visible as "blobs" or crescent-shaped features as seen in IRS 48 van der Marel et al. (2013) or HD 135344B Cazzoletti et al. (2018). Alternatively, hydrodynamic instabilities like the baroclinic instability (Klahr & Bodenheimer, 2003) or the vertical shear instability are able to form vortices (Manger & Klahr, 2018). In the presence of weak magnetic fields the magneto-rotational instability (MRI) triggers turbulence (Balbus & Hawley, 1991) and drives accretion flows. If the disk mid plane is effectively shielded from ionizing radiation the inner part of the disk becomes laminar, the so-called dead zone, with layered accretion on the surface level (Gammie, 1996). In the outer parts of the disk high energy photons may ionize the gas sufficiently to activate the MRI. The transition between the dead zone and the MRI-active region and thus the change in turbulent viscosity can also create ring structures (Flock et al., 2015). With all these possible substructure formation mechanisms at hand, it is of interest to identify markers of the presence of planets embedded in protoplanetary disks. A popular and wellstudied disk is the one around the Herbig Ae star HD 163296 at a distance of 101 pc (Gaia Collaboration et al., 2018). The appearance in the 1.25 mm continuum emission of the disk is dominated by two, already beforehand observed rings at a radial distance of 67 au and 100 au relative to the central star respectively (Isella et al., 2016, 2018). An additional faint ring has been detected at 159 au. An intriguing feature is a crescentshaped asymmetry within the inner gap located at 48 au (Huang et al., 2018). The feature itself is situated at a radial distance of Article number, page 2 of 19 55 au, thus with an offset of 7 au from the gap center (Isella et al., 2018). An image of the original observation is show in Fig. 13. The origin of such a structure is unknown and a preliminary Ring 3 - 155 au Ring 2 - 100 au Planet 3 - 137 au Ring 1 - 67 au Planet 2 - 83 au Planet 1 - 48 au Crescent-shaped asymmetry - 55 au Fig. 1: Sketch of the oberserved dust substructure of the HD 163296 system. The colored rings and the crescent-shaped feature mimic the dust while the grey dashed lines indicate the semi-major axes of the modeled planetary system. The labels introduced in this figure will be used throughout the text. model was presented in Zhang et al. (2018) involving planetdisk interaction. In these models asymmetries in the co-orbital region are common if the viscosity is low. Before the publication of the results from the DSHARP campaign, the HD 163296 disk was modeled by Liu et al. (2018). Their models incorporated 2D two-fluid hydrodynamical simulations with three planets in their respective positions matching the observed gaps. With synthetic images using radiative transfer calculations they could match the observed density profile with 0.46, 0.46, and 0.58 Jupiter masses for the three planets and a radially increasing turbulent viscosity parametrization. In the suite of simulations by Zhang et al. (2018) using hydrodynamical models with Lagrangian particles the proposed mass fits are 0.71, 2.18, and 0.14 Jupiter masses for an α-viscosity of 10−3 . In their lower viscosity models of α = 10−4 masses of 0.35, 1.07, and 0.07 Jupiter masses were fitted. Further observational constraints of the hypothetical two outer planets were provided by kinematical detections by Teague et al. (2018). Their model predicts masses of 1 and 1.3 Jupiter masses for these planets. Pinte et al. (2020) argue that velocity "kinks" observed in the CO observations with ALMA are evidence of nine planets in the DSHARP sample, including two planets in the HD 163296 disk at 86 au and 260 au. The signal-to-noise ratio is not sufficient to probe the inner gap at 48 au. In this paper we want to further explore the possibility of reproducing the observed structures by planet-disk interaction with a focus on the crescent-shaped asymmetry in the dust emission. This asymmetric feature has been present in the works discussed above but it has not been subject to more detailed analysis yet. Given the motivation of the crescent-shaped feature in the observation of HD 163296 we aim to constrain the visibility of such an agglomeration of dust caused by planet-disk interaction and its dependence on the physical parameters of the system like P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 planet mass, turbulent viscosity and dust size. In comparison to the study of Zhang et al. (2018) we employ two-dimensional hydrodynamical models with a fluid formulation of dust. Sec. 2 introduces the physical model and code setup as well as the post processing pipeline to predict observable features. In Sec. 3 we present the main results of our study. Sec. 4 compares these findings with previous works and addresses limitations of the model. In Sec. 5 we summarize the main results with concluding remarks. 2. Model All hydrodynamical models presented in this work were performed with the FARGO3D multi-fluid code (Benítez-Llambay & Masset, 2016; Benítez-Llambay et al., 2019) making use of an orbital advection algorithm (Masset, 2000). The code is based on the public version of FARGO3D with the addition of allowing a constant dust size throughout the simulation and a spatially variable viscosity. 2.1. Basic equations The FARGO3D code solves the conservation of mass (Eqs. [1] and [3]) and conservation of momentum (Eqs. [2] and [4]) for gas and dust in our model setups:   ∂t Σg + ∇ · Σg ug = 0 X h   i Σg ∂t ug + ug · ∇ ug = −∇P + ∇ · Π − Σg ∇Φ − Σi fi (1) (2) i Σd,i 
∂t Σd,i + ∇ · Σd,i ud,i + ji = 0
 ∂t ud,i + ud,i · ∇ ud,i = −Σd,i ∇Φ + Σd,i fi (3) (4) Here, Σg denotes the gas surface density, Σd,i the corresponding dust species, P = Σg c2s the gas pressure, linked to the density by a locally isothermal equation of state with the sound speed cs , Φ the gravitational potential of the star and planets, Π the viscous stress tensor, fi the interaction forces between gas and dust and ug and ud,i the gas and dust P velocities respectively. Dust feedback is included by the term i Σi fi in Eq. 2. We consider turbulent mixing and diffusion of dust grains by using the dust diffusion implementation described in Weber et al. (2019). The corresponding diffusion flux ji can be written as !   Σd,i , (5) ji = −Di Σg + Σd,i ∇ Σg + Σd,i with the diffusion constant Di being proportional to the turbulent viscosity ν (Youdin & Lithwick, 2007): Di = ν  1 + St2i 2 . 1 + St2i (6) The Stokes number Sti of dust species i is proportional to the stopping time tstop and can be written as Sti = tstop ΩK = π ai ρd , 2 Σg (7) where ΩK is the Keplerian angular frequency, ai the dust grain size and ρd the material density of the grains. The gas-dust interaction is modeled using the Epstein drag law, which is expected to be valid if the particle size is smaller than the mean-free path of surrounding gas molecules. Here, the drag force is proportional to the relative velocities between the corresponding fluids (Whipple, 1972). The drag force fi can be expressed as  ΩK  (8) fi = − ud,i − ug . Sti 2.2. Disk model In the disk model the planet-disk interaction is implemented as an additional smoothed potential term (Plummer-potential) for each planet. The smoothing length is set to 0.6H, where H(r) = cs (r)/ΩK (r) is the pressure scale height at the radial distance r. The specific factor acts as a correction for 3D effects in the 2D simulation (Müller et al., 2012). Three planets are modeled in the simulations. The two outer ones are set to the locations indicated by Teague et al. (2018). The inner planet is put at the corresponding gap location while the mass is varied in the different runs. In the fiducial models, the planet locations and masses are rp = {48, 83, 137} au, Mp0 = {1.0, 0.55, 1.0} Mjup , respectively. From this point on we will refer to the three planets as planet 1, 2 and 3. The same notation will be used for the apparent ring structures, i.e. ring 1: observed or modeled ring at 67 au; ring 2 at 100 au (see Fig. 1). We chose lower the mass of planet 2 compared to the one predicted in Teague et al. (2018) since it allows a sufficiently massive ring 2 while not significantly disturbing the crescent-shaped asymmetry by its repeated gravitational interaction. For the fiducial model the corresponding parameter variations we set the planet mass Mp to its final value Mp0 at the beginning of the simulation. The mass growth time scale however is known to have an impact on the formation of disk structures, such as vortices (Hammer et al., 2019; Hallam & Paardekooper, 2020). We therefore investigate the robustness of the results by testing different planet growth time scales: " !# t 1 1 − cos π Mp0 (9) Mp (t) = 2 TG TG refers to the planet growth time scale ranging from 10T 0 to 500T 0 with T 0 denoting the orbital period at r0 = 48 au. By using this simplified growth prescription we do not directly model the accretion of gas onto the planets and we thus do not artificially remove any mass from the simulation domain. Furthermore, for all models the displacement of the center of mass by the influence of the planets is taken into account as an additional indirect term added to the potential. The semi-major axes of the planets are kept fixed throughout the whole simulation. The initial surface density profile is assumed to be a power law with an exponential cutoff: " !−p !s # r r Σg/d = Σg/d,0 exp − , (10) r0 r0 with r0 = 48 au and the initial gas and dust surface densities Σg,0 and Σd,0 . Similar to the models of Liu et al. (2018) we choose a surface density slope of p = 0.8. For the dust a sharper cutoff of s = 2 was chosen compared to the gas cutoff P of s
= 1. In all simulation runs an initial dust-to-gas ratio of i Σd,i /Σg = 0.01 is assumed. We approximate the disk thermodynamics with a locally isothermal equation of state. The model is parametrized through locally isothermal sound speed !1 H(r0 ) r 4 cs = ΩK r. (11) r0 r0 Article number, page 3 of 19 A&A proofs: manuscript no. main 1.91 10-1 0.51 0.99 Grain Size Stokes Numbers 100 0.26 10-2 0.14 0.07 10-3 0.04 α ( r) 10-3 100 Radius [au] 200 Fig. 2: Upper panel: azimuthally averaged values of the Stokes number St of model fid. The dashed lines indicate the initial values and the solid lines represent the state after 500 T 0 . Lower panel: prescribed α viscosity. The blue line visualizes the radially increasing α viscosity set in the model. This corresponds to a flared disk with flaring index 0.25. The aspect ratio H/r at r0 is set to a value of 0.05. Assuming a mean molecular weight of µ = 2.353 the mid plane temperature profile can be written in the following way r !− 1 µmp r 2 cs ≈ 25 T mid (r) = K, (12) kB r0 with the proton mass mp . The temperature at 48 au matches the findings of Dullemond et al. (2019). We assume a radially smoothly increasing turbulent viscosity profile, motivated by Eq. 4 in the work of Liu et al. (2018) and similar dead zone parametrizations presented in Pinilla et al. (2016) and Miranda et al. (2016): !" !#) ( αmax r−R 1 1− 1 − tanh − , (13) α(r) = αmin 1 − 2 αmin σr0 where the parameters R and σ are set to 144 au and 1.25 respectively. Similar to the values in Liu et al. (2018) the parameter R refers to the mid-point of the transition in α whereas σ defines the slope. 2.3. Boundary conditions The radial velocities in the ghost cells are set according to antisymmetric boundary conditions. On a staggered mesh this corresponds to vr (rghost ) = −vr (ract ), where rghost is the ghost cell and ract the equivalent mirrored cell on the active hydro mesh near the boundary. The value at the staggered boundary itself is set to zero. The azimuthal velocities are set to the initial Keplerian profile in the ghost zone. The surface density is extrapolated according to the density gradient exponent p. Additionally, wave damping is applied within 15% of the respective boundary radius. In this region the density is exponentially relaxed towards the initial values within 0.3 local orbital time scales. The procedure follows the wave damping boundary conditions used in de Val-Borro et al. (2006). We tested the robustness of the wave damping with respect to the chosen inner boundary condition. No significant wave reflections are detected for both a symmetric and anti-symmetric Article number, page 4 of 19 2.4. Code setup and parameters 0.02 10-2 10-4 20 inner boundary condition. The relative difference between these two approaches is on the order of 10−5 compared to the reflected perturbations of ≈ 10−2 without wave damping. Various ranges of individual disk parameters were considered for constraining the impact onto dust features in the simulations. Table 1 gives an overview over the performed simulations and their parameter choices. For most of the runs a resolution of 560 radial and 895 azimuthal cells was chosen where the grid is logarithmically spaced in radial direction. Compared to the local disk scale height a ratio of roughly 7 cells per scale height is achieved in each direction. The models with the suffix _dres were run with a resolution of 14 cells per scale height. The low resolution runs are subject to a larger numerical diffusion compared to the high resolution models. Although these runs can be used for mass estimates of the crescent-shaped feature, the results should be taken with caution concerning the dust substructure lifetime. Whenever this difference becomes significant we default to the high resolution runs. Further details are given in appendix B. For all simulations the mass of the central star is set to 1.9 M⊙ The parameters αmin and αmax are set to 10−5 and 5 · 10−3 which results in a value of α(r0 ) ≈ 2 · 10−4 at the location of r0 = 48 au. The radial profile of α(r) is shown in the lower panel of Fig. 2. The dust is sampled by 8 separate fluids with Stokes numbers spaced logarithmically ranging from 1.3 · 10−3 to 1.3 · 10−1 at r0 = 48 au. In the code, the equivalent grain size at r0 is applied to the whole domain and is kept constant throughout the whole simulation. With an initial value of Σg,0 = 37.4 g/cm2 the minimum and maximum grain sizes are amin = 0.19 mm and amax = 19 mm. No dust size evolution is modeled here. It should be noted that these state the initial values of St and changes with time depending on the gas surface density, as shown in the upper panel of Fig. 2. The most prominent change of an increase of about one order of magnitude in St occurs at the gap carved by planet 1 after 500 orbits. The majority of models neglects dust feedback to the motion of the gas. Each dust species can thus be scaled in density individually without violating the validity of the dynamical features. Per default, the simulations are executed until 1000 T 0 . Simulation runs with a nonzero growth time scale T G and the model fid_dres are run until 2000 T 0 . After ≈ 900 orbits the gas and dust structure converges. The crescent-shaped asymmetry builds up to a stable niveau after less than 100 orbits. We refer to Sec. 3 for more details. 2.5. Radiative transfer model and post processing In order to compare the results of the hydrodynamical simulations with the observational data synthetic images are produced with RADMC-3D (Dullemond et al., 2012) and the CASA package (McMullin et al., 2007). The results of the hydrodynamical simulations have to be extended to three dimensional dust density models which then serve as input for the radiative transfer calculations with RADMC-3D. Using the given Stokes numbers of the hydro model, the respective dust sizes are computed via Eq. 7. The number density size distribution of the dust grains follows the MRN distribution n(a) ∝ a−3.5 (Mathis et al., 1977) where a is the grain size. Dust settling towards the mid-plane is considered following the P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 Table 1: Relevant simulation runs and their respective numerical parameters. Simulation fid fid_dres hr4 hr45 hr55 hr6 taper10 taper50 taper100 taper500 p1m1 p1m2 p1m3 p1m4 p1m5 p1m1fb p1m2fb p1m3fb p1m4fb p1m5fb p1m6fb p1m6fb_dres p2m1 p2m2 p2m3 p2m4 p2m5 p2m6 ecc1 ecc2 ecc3 ecc4 ecc5 cut1 cut2 cut3 cut4 cut5 alpha1 alpha1_dres alpha2 alpha3 alpha3_dres alpha4 alpha4_dres alpha5 alpha6 alpha6_dres H/r Mpl1 [Mjup ] 0.05 0.05 0.04 0.045 0.055 0.06 1 1 Mpl2 [Mjup ] e 0.55 0.55 Nr 560 1120 Nφ rcut [au] α TG −4 895 1790 2 · 10 2 · 10−4 0 10 50 100 500 0.20 0.36 0.52 0.68 0.84 0.20 0.36 0.52 0.68 0.84 1120 1790 0.30 0.42 0.54 0.66 0.78 0.90 0.02 0.04 0.06 0.08 0.10 150 175 200 225 250 1120 1790 1120 1790 1120 1790 1120 1790 1 · 10−5 1 · 10−5 1 · 10−4 5 · 10−4 5 · 10−4 1 · 10−3 1 · 10−3 2 · 10−3 3 · 10−3 3 · 10−3 Notes. Blank spaces assume an identical parameter as the fiducial model fid. The simulation labels starting with hr denote models with a variation in aspect ratio. The models starting with p include a change in the mass of planet 1 or 2 depending on the following digit. A suffix of fb describes models with dust feedback activated. All models use a radially varying alpha viscosity model with α = 2 · 10−4 except the simulations denoted by alpha1 to alpha6 where a radially constant value of α was chosen. diffusion model of Dubrulle et al. (1995) r α Hd = H α + St (14) Here, we assume a Schmidt number on the order of unity. The grid resolution for the radiative transfer is identical to the hydrodynamical mesh. In polar direction the grid is expanded by 32 cells which are equally spaced up to zlim = π/2 ± 0.3. The vertical disk density profile is assumed to be isothermal and the conversion from the surface density to the local volume density is calculated as follows: ρcell = √ ! π zlim · Hd · erf −1 √ 2πHd 2Hd 2 Σ  erf  z √+ 2Hd  − erf z+ − z−  z √− 2Hd  . (15) Article number, page 5 of 19 A&A proofs: manuscript no. main Table 2: Optical depth fitting results for ring 1. 2 2 Constant Σg,0 (r = 48 au) [ cmg 2 ] Σg,0 (r = 1 au) [ cmg 2 ] amin [mm] amax [mm] κmin [ cmg ] κmax [ cmg ] Grain size d/g ratio 37.4 1.31 827.7 28.91 0.191 0.007 19.1 0.67 1.32 0.40 0.11 3.30 Notes. Results for both the high and low mass model as derived from the optical depth fitting of ring 1 in Sec. 3.2. Both the initial gas surface densities Σg,0 at 1 and 48 au are listed. Depending on the gas density a different dust grain size distribution with the minimum and maximum size of amin and amax was chosen. Their respective opacity values are listed under κmin and κmax . Table 3: Gaussian fit results of the simulated dust rings for all Stokes numbers, St. a [mm] St wring1 [au] wring2 [au] rring1 [au] rring2 [au] Mring1 [Mearth ] Mring2 [Mearth ] 0.2 0.4 0.7 1.4 2.6 5.1 9.9 19.1 sum (high mass model) sum (low mass model) 1.3 · 10−3 2.6 · 10−3 5.0 · 10−3 9.6 · 10−3 1.9 · 10−2 3.6 · 10−2 6.9 · 10−2 1.3 · 10−1 4.09 ± 0.09 3.13 ± 0.08 2.27 ± 0.07 1.68 ± 0.05 1.28 ± 0.04 1.00 ± 0.03 0.80 ± 0.03 0.64 ± 0.02 10.5 ± 0.17 7.10 ± 0.09 4.80 ± 0.04 3.34 ± 0.01 2.37 ± 0.01 1.76 ± 0.01 1.33 ± 0.01 1.01 ± 0.01 62.26 ± 0.09 62.05 ± 0.08 61.90 ± 0.07 61.86 ± 0.05 61.89 ± 0.04 61.92 ± 0.03 61.93 ± 0.03 61.94 ± 0.02 100.18 ± 0.16 99.90 ± 0.09 99.25 ± 0.04 98.87 ± 0.01 98.69 ± 0.01 98.67 ± 0.01 98.67 ± 0.01 98.66 ± 0.01 2.1 2.9 4.0 5.5 7.9 11.5 17.3 28.0 79.3 11.4 1.9 2.9 4.7 8.0 12.3 17.0 22.2 26.9 95.9 13.8 Notes. Ring widths are denoted by wring1,2 while the ring position is labeled rring1,2 . The dust size individual mass trapped in the rings and the total mass are shown in the last two columns. The error function term is a correction for the limited domain extend in the vertical direction that would otherwise lead to an underestimation of the total dust mass. Similarly, the second correction term accounts for the finite vertical resolution, especially important for thin dust layers with strong settling towards the mid-plane. The coordinates z+ and z− denote the cell interface locations in polar direction along the numerical grid. For each grain size bin and a wavelength of 1.3 mm the corresponding dust opacities were taken from the dsharp_opac package which provides the opacities presented in Birnstiel et al. (2018). These opacities are based on a mixture of water ice, silicate, troilite and refractory organic material. The grains are assumed to be spherically shaped and to have no porosity. In the RADMC-3D model, the central star is assumed to have a mass of 1.9 solar masses with an effective temperature of 9333 K which results in a luminosity of 2.62 · 1034 erg/s. The system is assumed to be at a distance of 101 pc. For the dust temperature calculation a number of nphot = 108 photon packages and for the image reconstruction nphot_scat = 107 photon packages were used. The thermal Monte-Carlo method is based on the recipe of Bjorkman & Wood (2001). We include isotropic scattering in the radiative transfer calculation. A comparison between the prescribed gas temperatures and the computed dust temperatures is given in Appendix D. For simulating the detectability of the various features present in the model we use the task simalma from the CASA5.6.1 software. A combination of the antenna configurations alma.cycle4.8 and alma.cycle4.5 was chosen. The simulated observation time for configuration 8 is 2 hours while the more compact configuration is integrated over a reduced time with a factor of 0.22 corresponding 0.44 hours. We employed the same cleaning procedure as made available Article number, page 6 of 19 in Isella et al. (2018) to reduce the artificial features from the incomplete uv-coverage and to allow a comparison to the observation. The procedure involves the CASA task tclean with a robust parameter of -0.5 and manual masking of the disk geometry. Consistent with the radiative transfer model, the observed wavelength is simulated to be at 1.3 mm, corresponding to ALMA band 6. 3. Results In the following parts the outcome of the simulation runs listed in Table 1 will be presented and analyzed. First, the variety of substructures emerging from the interaction of gas, dust and the three planets will be described. Afterwards parameter dependence and observability will be addressed. 3.1. Dust substructure overview A variety of dust substructures emerges from the planet-disk system during its dynamical evolution. Fig. 3 and Fig. 4 show the dust surface density structure for a selected parameter space of aspect ratio and the turbulent viscosity, characterized by α. In the following analysis of the crescent-shaped asymmetries simulation snapshots after 500 orbits at 48 au are compared with each other since their evolution is comparable for all resolutions. Most prominently multiple rings form in most cases. As expected, fluids with larger Stokes numbers St and thus larger grain sizes exhibit thinner rings and more concentrated substructures. Especially, nonaxisymmetric features can be seen mostly for a ≥ 2.6 mm or St > 10−2 . For smaller dust sizes the dust is better coupled to the gas and resembles its structure more P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 closely. Additionally, if a crescent-shaped asymmetry is present, it is situated in the gap caused by the most inner planet at 48 au for the majority of the parameter space. Further asymmetries appear if the α-viscosity is radially constant with values of α ≤ 1 · 10−3 or for low values of the aspect ratio h. Dust is preferably trapped in the Lagrange point L5. In Fig. 3 rings and asymmetries become weaker with increasing aspect ratio. A crescent-shaped feature at the L5 position is present for all aspect ratios while a second similar asymmetry at the L4 point appears for values of H / r < 0.05. The crescent weakens for smaller grain sizes. Also the second prominent ring beyond planet 2 at 83 au weakens clearly for larger aspect ratios. In the combination of the largest grain sizes and H/r the ring completely vanishes. To highlight the importance of the turbulent α viscosity parameter a subset of results with radially constant values of α are shown in Fig. 4. Not surprisingly, larger values of alpha generally lead to a more diffuse and symmetric distribution of dust. Below α = 5 · 10−4 no concentric rings form due to vortices in the gas. Crescent-shaped features in both Lagrange points of the innermost planet are visible in the very low viscosity case of α = 10−5 . A sufficiently large viscosity on the other hand also leads to the disappearance of the second ring in the limit of larger grains and Stokes numbers, similar to the large aspect ratio in Fig. 3. The collection of these results also stresses the issue with a radially constant α viscosity with respect to the observed HD 163296 system. In order to reproduce a nonaxisymmetric feature in the vicinity of the inner planet and a smooth ringshaped outer structure, a radially increasing value of α would be the natural choice. This is also consistent with an embedded dead zone at the inner part of the disk and a more active outer disk region with a higher degree of ionization (Miranda et al., 2016; Pinilla et al., 2016). We thus chose a radially increasing α viscosity for all remaining simulation runs. Models with a variation in the dust cutoff radius only show little changes in the resulting dust structures. Only this subset of the possible parameter space already exposes the degeneracy of the emerging substructures with respect to the chosen disk models. work excludes contributions from the prominent nonaxisymmetric structures. Several properties become apparent: – Ring 1 is wider than in the simulated profile. – The peak location of ring 1 is located further outward with respect to the simulated one. – The peak value of the optical depth of ring 2 from smaller grains and Stokes numbers is lower in the simulations compared to the estimated value from the observation. A partial explanation of these differences could be that the dustto-gas ratio could become larger than in the models so that dust feedback shifts the ring further outward and spreads the ring as shown in Weber et al. (2018) and Kanagawa et al. (2018). Here, the dust-to-gas ratio is not sufficient to cause a significant effect in the models including dust feedback. A conclusion of Dullemond et al. (2018) was that the optical depth observed in the DSHARP survey were remarkably close to unity and that the rings were optically thin. Later Zhu et al. (2019) argued that dust scattering could account for this phenomenon and that the actual optical depth could be larger. In the case of HD 163296 the mass hidden in ring 1 could be thus larger than expected. For ring fitting we use a Gaussian: ! (r − r0 )2 Σfit (r) = A exp − (16) 2w with the peak value A, the ring location r0 and the ring width w. In Fig. 6 the fitted ring widths from model fid are plotted and compared to the observed values in Dullemond et al. (2018). Grain sizes of the high mass model are used. The width of ring 2 is matched close to a grain size of 1 mm. On the other hand, the width of ring 1 is not reached with the parameters chosen in our models. It should be noted that the gas ring width is about 8 au, just slightly larger than the observed value of about 7 au. Smaller Stokes numbers could in principle reproduce these findings. The equivalent model p1m6fb_dres including dust feedback shows no significant differences. 3.2. Rings Before quantifying the nonaxisymmetric feature in the model, we want to compare properties of the ring structures with previous works and the observations. The general procedure is to azimuthally average the dust density maps after 1000 orbits at 48 au and to invoke Gaussian fitting of the dust rings, comparable to Dullemond et al. (2018). Since the ring structure is approximately converged after 900 orbits the following procedure is based on the snapshots at 1000 orbits. Fig. 5 shows the results of the high resolution fiducial model fid_dres as an example of the radial surface density structure. The increased resolution was chosen since the lower resolution models may overestimate the ring width and the trapped dust mass due to numerical diffusion (appendix B). In the upper panel the dust species is rescaled to gas density peak of ring 1. The dust rings are clearly thinner than the gaseous envelope and the ring width decreases with increasing dust grain sizes due to the stronger drift. With the corresponding opacities κi a rough estimate of the resulting optical depth can be computed by τsim = κi Σi , where i denotes the dust species index. The results of this estimate are displayed in the lower panel of Fig. 5. All optical depths are rescaled to the values at the position of ring 1 from the profile τobs derived in Huang et al. (2018). The profile provided in their 3.3. Surface density estimation With ring 1 being the most prominent substructure in the observed system, its estimated lower bound for the surface density would be a reasonable choice for rescaling the simulated dust density maps. Furthermore, the crescent-shaped feature of interest is located closely to ring 1. In order to estimate a minimum mass of trapped dust in this feature an appropriate normalization of the density with respect to ring 1 would be a natural choice. There are two possible methods in achieving a simple normalization. First, we rescale the sum of all azimuthally averaged dust densities so that the combined optical depth at the location of ring 1 equals the observed value. To maintain the validity of the dynamics of the system, the Stokes number corresponding to the grain size of a fluid has to be unmodified by this process. The immediate consequence is then a change in the dust-to-gas ratio if the densities are rescaled, since a change in the gas surface density with a constant grain size would modify the Stokes number (see Eq. 7). With a change in the dust-to-gas ratio the dust dynamics only remain comparable if no dust feedback is considered. The second choice would be to maintain an initial dust-to-gas ratio of 0.01 and to change the dust grain size and the gas surface density. A modification of the grain size affects in turn the Article number, page 7 of 19 A&A proofs: manuscript no. main St = 1.3 10 1 h = 0.045 h = 0.050 h = 0.055 h = 0.060 100 0 100 0.04 0.03 0.02 0.01 0.15 100 0.10 0 0.05 100 0.5 0.4 0.3 0.2 0.1 100 0 100 100 0 100 100 0 100 100 0 au 100 100 0 100 100 0 Surface density [g/cm2] St = 1.9 10 2 St = 5.0 10 3 h = 0.040 100 Fig. 3: Shown are dust surface density maps for a subset of three fluids with varying values of the aspect ratio h = H/r at 500 orbits at 48 au. Crescent-shaped asymmetries are visible for all aspect ratios. Large values of h weaken dust accumulation in the co-orbital regions of the planets. The white crosses mark the position of planet 1. dust opacities and thus the optical depth. Consequently, the process of generating τsim , rescaling it to τobs and inferring the corrected dust surface densities has to be iterated until convergence is achieved. Keeping the dust-to-gas ratio constant is important for the model runs with dust feedback enabled. Table 2 provides relevant results from these two approaches which will be denoted by the high and low mass model in the following parts. With the unchanged grain sizes the dust-to-gas ratio diminishes to ≈ 2.4 · 10−3 . For the iterative approach with the dust-to-gas ratio unchanged, the grain size distribution shifts towards smaller grains with a maximum size amax ≈ 0.67 mm and a minimum size amin = 0.007 mm. Results of the Gaussian fits onto the dust rings of the fiducial model are listed in Table 3 for all simulated dust species. The inferred ring widths are decreasing with increasing grain sizes and Stokes numbers. The peak maxima shift towards the star for larger grain sizes since the dust drift becomes more dominant in this regime. The total dust mass of all species is computed for both the low and high mass model. Dullemond et al. (2018) found masses of 56 Mearth and 43.6 Mearth for ring 1 and ring 2 respectively while assuming 1 mm grains with equal opacities values as used in the model presented here. The results of both the low and high mass model encompass the values of Dullemond et al. (2018). The high mass model will be the preferred choice in the following diagrams and analysis since the opacity is dominated by smaller grain sizes compared to the low mass model. The choice is motivated by the broad ring structures apparent in the observations. Article number, page 8 of 19 3.4. Secondary planet mass In the models p2m1 - p2m6 the mass of the secondary planet at 83 au is varied in order to verify its impact on the ring structure. The results of Teague et al. (2018) indicate a planet mass of 1 Mjup within an error margin of 50%. In our model runs we chose a mass of 0.55 Mjup for planet 2 since a larger mass causes a stronger dissipation of the crescent-shaped asymmetry. Lower masses significantly decreased the dust content in ring 2 and thus the fiducial planet mass value was chosen as the sweet spot between a strong ring contrast and an maximized asymmetric dust accumulation. Further details are given in appendix A. 3.5. Asymmetries Of particular interest is the crescent-shaped asymmetry in the vicinity of ring 1 in HD 163296. Such a feature arises naturally in planet-disk interaction models including dust in the form of dust trapping in a Lagrange point of the gap carving planet. In this case planet 1 is responsible for dust trapping in the trailing L5 point which is also visible in Fig. 3 and Fig. 4 for a significant subset of the parameter space. An equivalent result was presented in Isella et al. (2018). In the following subsections we aim to perform a more extensive analysis of this feature in order to constrain physical properties of the dynamical system. P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 5 = 1.0 10 = 5.0 10 4 4 = 1.0 10 3 = 2.0 10 3 = 3.0 10 3 0.04 0.03 0.02 0.01 St = 5.0 10 3 100 0 100 0.15 St = 1.9 10 2 100 0.10 0 0.05 100 St = 1.3 10 1 100 Surface density [g/cm2] = 1.0 10 0.4 0 0.2 100 100 0 100 100 0 100 100 0 100 100 0 100 au 100 0 100 100 0 100 g cm 2 ] Fig. 4: Dust surface density maps for different values of a radially constant α viscosity after 500 orbits at 48 au. For α ≤ 10−4 strong asymmetries are present. ring 2 weakens or vanishes for large viscosities. The white crosses mark the position of planet 1. Optical Depth 2 6 · 10 −3 9 6 · 10 −3 −2 St = 3 6 · 10 −1 St = 1 3 · 10 2 St = . St = . 10 . 10 1 . 8 gas 10 0 10 0 Width [au] Surface Density [ / 10 1 3 · 10 −3 2 6 · 10 −3 −3 St = 9 6 · 10 −2 St = 3 6 · 10 −1 St = 1 3 · 10 St = . St = . Ring 1 (Dullemond+2018) Ring 2 (Dullemond+2018) Ring 1 (Fiducial model) Ring 2 (Fiducial model) 6 4 . 10 2 . -1 . τ 10 obs -2 25 50 75 100 125 150 175 200 Radius [au] Fig. 5: Upper panel: azimuthally averaged dust and gas surface densities of model fid_dres after 1000 orbits at 48 au. Dust density profiles are normalized to the gas peak value at the location of ring 1. Lower panel: azimuthally averaged optical depths of model mid after 1000 orbits at 48 au compared to the observed optical depth τobs . Simulated optical depths are normalized to τobs at the location of ring 1. 3.5.1. Structure & dust feedback In Fig. 7 in the left panels dust surface density maps are shown in polar coordinates for four different dust fluids of model fid_dres. The region is focused around the co-orbital region of planet 1. Clearly, dust is concentrated in the trailing Lagrange point L5 of the Jupiter mass planet at 48 au. Several trends become apparent: Not surprisingly, dust grains are trapped more efficiently for larger grain sizes and Stokes numbers due to the stronger drift. Furthermore, the shape is more elongated for 103 Grain size [µm] 104 Fig. 6: Ring widths of model fid_dres after 1000 orbits at 48 au. Simulated values are compared to the inferred ring widths of Dullemond et al. (2018), displayed as the horizontal dashed lines. smaller Stokes numbers. Does the dynamics of this feature change with the consideration of a dust back-reaction onto the gas? An example of the impact from dust feedback (model p1m1fb_dres) is given in Fig. 7 on the right hand side with the same parameters as model fid_dres. The crescent-shaped feature exhibits different structures compared to model fid_dres at the location of L5. The dust back-reaction onto the gas triggers and instability leading to fragmentation of the dust feature. Unlike the fiducial model, the additional crescent in the leading Lagrange point L4 of planet 1 is more pronounced than the L5 feature for smaller grain sizes and Stokes numbers. In the explored parameter space dust is significantly trapped for small aspect ratios or very low values of α < 10−4 without the modeling of dust feedback. With the fiducial set of parameters the L4 feature slowly dissipates after more Article number, page 9 of 19 A&A proofs: manuscript no. main 80 St = 5.0 10 60 40 80 3 2.0 1.5 1.0 0.5 St = 9.6 10 60 40 80 2 2.0 1.5 1.0 0.5 0.0 St = 3.6 10 60 40 1 80 St = 1.3 10 60 40 2 /3 /3 0 [rad] /3 2 /3 2 /3 /3 0 [rad] /3 2 /3 Surface density (Normalized) 3 2.0 1.5 1.0 0.5 2.0 1.5 1.0 0.5 0.0 Fig. 7: Dust surface density maps in polar coordinates for model fid_dres without dust feedback (left hand side) and p1m6fb_dres with dust feedback after 500 orbits at 48 au for four different initial Stokes numbers and dust sizes. A dust agglomeration around the location of the trailing Lagrange point L5 is present for all Stokes numbers. Unlike the nonfeedback case the right panels show that dust is trapped more efficiently in the leading Lagrange point L4 for smaller Stokes numbers. For St . 3.6 · 10−2 the L4 feature is more pronounced than the dust over density around L5. Dust feedback leads to an instability at the L5 point and fragments the asymmetric feature. Dust densities are normalized to the peak value of ring 1. than a 1000 orbits. A closer look on the radial and azimuthal extent of these dust shapes is provided in Fig. 8 for all simulated dust fluids. The surface densities of the crescent-shaped asymmetry are normalized to its maximum value. The radial and azimuthal width increases with decreasing values of the Stokes number. In the lower panel the density peak at the trailing L5 dominates the leading peak at L4 for all dust species. For smaller Stokes numbers and grain sizes the density maximum moves away from the planet location similar to the peak shift of the concentric dust rings described in Sec. 3.2. In Fig. 9 we display the results of model p1m6fb_dres. Contrary to the nonfeedback case the density peak at L4 becomes significant for St ≤ 3.6 · 10−2 . In the vicinity of the L5 point two density peaks appear due to the fragmentation by dust feedArticle number, page 10 of 19 back. The azimuthal gas density profile reveals the momentary location of the spiral wakes caused by the planets as well as the gas accumulation around planet 1 itself at φ = 0. 3.5.2. Dynamical stability In principle the crescent-shaped features in the co-orbital region are subject to diffusive processes like dust diffusion due to turbulent mixing or gravitational interaction with the planetary system (e.g. eccentric orbits). The dust trapping mechanism has to counteract these disruptive forces for the feature to be dynamically stable. Fig. 10 corroborates this line of argument. The stability of the L5 feature is sensitive to the local value of the α viscosity. Given P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 . 3 6 · 10 −2 . 52 . 1 3 · 10 −3 . π −2 3 π/ − 3 π/ φ 0 [rad] 3 π/ 2 3 π π/ Mna Fig. 8: Radial and azimuthal cut through the maximum of the crescent-shaped asymmetry around L5 of the fiducial model fid_dres after 500 orbits at 48 au. The color map indicates the different initial Stokes numbers of the dust fluids. The dashed lines represent the gas density. The surface densities are normalized to their respective maximum value in the co-orbital region of planet 1. The light blue and light green vertical lines indicate the Lagrange points L5 and L4 respectively. 1.0 gas density 3.6 · 10 −2 1.9 · 10 −2 α = 1 · 10 −3 9.6 · 10 −3 5.0 · 10 −3 2.6 · 10 −3 1.3 · 10 −3 α = 3 · 10 −3 200 400 600 Orbits (at 48 au) 800 1000 . 6 9 · 10 −2 . Initial Stokes Number / Σ Σ max 1.3 · 10 −1 6.9 · 10 −2 1 3 · 10 −1 0.5 3 6 · 10 −2 . 0.0 44 46 1.0 48 r [au] 50 1 9 · 10 −2 52 . 9 6 · 10 −3 . gas density 5 0 · 10 −3 Σ Σ max . / = 5 · 10 −4 ] 0.0 − Mearth 0.5 . α [ 2 6 · 10 −3 Mearth / Σ Σ max . Mna . 5 0 · 10 −3 ] 9 6 · 10 −3 gas density [ 50 Mna 1.0 48 r [au] = 1 · 10 −5 ] 46 Mearth 44 1 9 · 10 −2 α [ 0.0 Mearth 6 9 · 10 −2 5 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0 4 3 2 1 0 [ 0.5 Mna . Initial Stokes Number / Σ Σ max 1 3 · 10 −1 ] gas density Initial Stokes Number 1.0 2 6 · 10 −3 0.5 . Fig. 10: Development of the trapped dust mass in the L5 region for different α viscosities over time with a resolution of 14 cells per scale height. Dust masses Mna in the nonaxisymmetric feature denoted are integrated down to a cut-off value of 10−2 Σmax and normalized to the high mass model stated in Table 2. 1 3 · 10 −3 . 0.0 − π −2 3 π/ − 3 π/ φ 0 [rad] 3 π/ 2 3 π/ 3.5.3. Dust mass π Fig. 9: Radial and azimuthal cut through the maximum density value in the co-orbital region of the model p1m6fb_dres with dust feedback after 500 orbits at 48 au. The light blue and light green vertical lines indicate the Lagrange points L5 and L4 respectively. Depending on the location of the feature the density cuts are normalized to the peak value in the region of L4 or L5. Smaller dust with St . 3.6 · 10−2 preferably concentrates in L4. a value of α = 10−5 the feature remains stable throughout almost the entirety of the simulation whereas α = 10−3 shortens the existence down to about 300 dynamical time scales. For larger viscosities no discernable feature develops and the co-orbital region simply empties its dust content from the initial condition. Studying the results of the models p1m1 to p1m5 we find that below about 0.4 to 0.5 Jupiter masses for planet 1 no stable feature forms since the gravitational interaction is not sufficient to enforce dust trapping in the Lagrange points. Qualitatively the feature life time is thus very sensitive to the given physical parameters. It should be noted that the absolute value in dynamical time scales could be underestimated due to the numerical diffusion present in the lower resolution runs. This additional diffusive effect prevents a stable trapping region around L4 and L5. Further details are provided in appendix B. In Fig. 10 we therefore plotted results with 14 cells per scale height. A substantial amount of dust can be trapped in the feature at L5. The total mass ranges from roughly 1 Mearth for the low mass model to 10 − 15 Mearth for the high mass model. An overview of relevant parameters and their impact on the trapped dust mass is given in Fig. 11. Generally, if a sufficiently stable crescent-shaped asymmetry develops, the order of magnitude of the trapped dust mass is comparable for all parameters. Here, we only consider the three largest dust species since smaller grains are prone to be weakly trapped. Only the low resolution simulations are compared to each other in this parameter study. The mass of the crescent-shaped feature was averaged over 200 orbits starting from 100 T 0 when convergence is reached. Looking at the aspect ratio dependence, we find that an increase in H/r also leads to a higher dust mass of the L5 feature. A more massive planet also causes an increase in the trapped dust mass. The lighter the planet, the less stable the agglomeration of dust becomes. As mentioned above, for less than about 0.4 to 0.5 Mjup no stable feature forms at L5. Considering the influence of the α viscosity parameter a local value of α ≥ 10−3 causes a significant loss of dust mass. For α ≥ 2 · 10−3 no feature develops. Simulation results with a radially constant value of α are used here Finally, the introduction of an eccentric planetary orbit leads to an almost linear decrease of the trapped dust mass with respect to the eccentricity value. For values ≥ 0.06 no feature forms for dust Stokes numbers of 3.6 · 10−2 and below. An increase of the mass of planet 2 has a small influence on the Article number, page 11 of 19 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.040 0.045 0.050 0.055 0.060 Aspect ratio 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.00 0.02 0.04 0.06 0.08 Eccentricity 0.6 0.8 1.0 Planet 1 mass [Mjup] 10 5 10 4 10 3 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.30 0.42 0.54 0.66 0.78 0.90 0.5 Planet 2 mass [Mjup] 101 102 Growth time scale 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 1.3 10 1 6.9 10 2 3.6 10 2 Initial Stokes number Trapped dust mass [Mearth] Trapped dust mass [Mearth] A&A proofs: manuscript no. main Fig. 11: Trapped dust masses in earth masses in the L5 feature for the largest three inital Stokes numbers. The mass values are averaged over the approximately constant mass time frame starting from 100 T 0 . Dust masses are normalized according to the high mass model stated in Table 2. Missing data points are equivalent to a nonexistence of a stable crescent-shaped asymmetry in the co-orbital region around 48 au. trapped dust mass. The continuous gravitational interaction perturbs the crescent-shaped asymmetry decreases the amount of mass trapped in the feature. The difference between 0.3Mjup and 0.9Mjup however accounts to roughly 5% of trapped dust mass. 3.5.4. Growth time scale The growth time scale of the planet mass can have a significant impact on the formation of vortices (Hammer et al., 2017, 2019; Hallam & Paardekooper, 2020). Therefore, simulation runs including longer growth time scales with values of T G = 10, 50, 100 and 500 orbits were performed. Qualitatively, the results are mostly unaffected by the choice of T G . In Fig. 11 the masses only deviate about 10 % from the fiducial model. With the longest growth time scale of 500 orbits smaller grains are trapped more efficiently while the mass contained in the largest grains decreases slightly. Deviations from the fiducial model for these long growth time scales could also be caused by a loss of dust content and local redistribution of grain sizes due to dust drift in the simulation domain. More details are given in Appendix C. The formation of vortices is sensitive to the planet growth time scale since the vortex smooths out the gap edge and reduces the Article number, page 12 of 19 steepness of the corresponding edge slope which in turn weakens the Rossby wave instability (Hammer et al., 2017). This effect is important for longer planet growth time scales and leads to weaker, elongated vortices. In the context of dust trapping in the Lagrange points however, the process takes place in the coorbital region and the amount of mass concentrating in the asymmetries is determined by the initial dust content available within this region (Montesinos et al., 2020). Since the dust trapping here is related to the horseshoe motion in the co-orbital region, it is a different mechanism and the evolution of the gap edge does not seem to have a major influence on the dynamical origin of the asymmetric features. 3.6. Synthetic images The main question remains if dust trapping in the L5 point in the models presented could explain the observed feature in HD 163296. We apply the procedure described in Sec. 2.5 and thus extend the surface density maps to three dimensional grids, perform dust radiative transfer calculations RADMC-3D and simulate the observation with ALMA by using the CASA package accordingly. Snapshots of the density maps for these synthetic images are shown in Fig. 12. Both the high reso- P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 St = 1.3 10 1 feedback (500 T0) fiducial (1000 T0) feedback (1000 T0) fiducial (2000 T0) 100 0 100 1.0 0.8 0.6 0.4 0.2 1.0 0.8 0.6 0.4 0.2 100 0 100 1.0 0.8 0.6 0.4 0.2 100 0 100 100 0 100 100 0 100 100 0 au 100 100 0 100 100 0 Surface density (Normalized) St = 1.9 10 2 St = 5.0 10 3 fiducial (500 T0) 100 Fig. 12: Comparison of the dust surface density map for three different dust sizes of the models fid_dres (fiducial) and p1m6fb_dres (feedback) at different simulation times. The surface densities are normalized to the peak value at ring 1. lution fiducial model fid_dres and the dust feedback model p1m6fb_dres are used. The resulting dust density grid crucially depends on the vertical dust scale height Hd and thus the dust settling prescription (see Eq. 14). With the assumption of constant grain sizes throughout the disk the local Stokes number is used for the calculation of Hd . The grain sizes distribution is the same as the one in the simulation runs. Depending on the low or high mass model, the grain sizes and opacities and densities were adjusted accordingly. We chose to leave α as a free parameter for the vertical dust settling recipe which will be denoted as αz . In Fig. 13 synthetic images from various snapshots of the simulation are compared to the observation. The images qualitatively reproduce the observed features. The key difference with respect to the crescentshaped feature is the more elongated shape compared to the fiducial model. Furthermore, the models produce a radially symmetrically located feature while the observed one is situated closer to ring 1. After 500 and 1000 orbits at 48 au the feature at L4 is still visible due to the slow dissipation at the L4 point. In the image computed from the snapshot at 2000 orbits of the simulation fid_dres, only the L5 feature appears, as also seen in Fig. 12. As already discussed in Sec. 3.5.1 and Fig. 7 a significant amount of dust agglomerates in the L4 region for smaller grain sizes and Stokes numbers. This effect is clearly visible in the synthetic image of p1m6fb_dres at 500 orbits with dust feedback enabled in Fig. 13. No such feature is present in the observation. The fragmentation of the crescent-shaped asymmetry around L5 becomes more apparent in the later stages of the simulation. After 1000 orbits dust is concentrated in clumps of small azimuthal extent. These features are substantially different from the observation. Since the high mass model assumes a lower dust-to-gas ratio, the synthetic images created from the high mass model and the respective grain size distribution only serve as a comparison of the visibility of such features in the co-orbital region. Additionally, synthetic observations based on the low mass model are shown in Fig. 13. As expected, the larger values of the opacity for the grains with the maximum Stokes numbers compared to the high mass model leads to much narrower rings and substructures. The high mass model is thus more suitable for the HD 163296 system. Focusing on ring 2, the intensity is slightly lower compared to the observations. This result is consistent with Fig. 5 where the optical depth of ring 2 does not reach the derived values of Huang et al. (2018) for smaller grains. Looking at the inner part of the disk within the gap at 48 au, the simulated images show a secondary gap caused by planet 1 which is not present or visible in the observed structure. Furthermore, the influence of vertical mixing of dust grains can be investigated with the three synthetic images in Fig. 14 of the nonfeedback model. The model run with an increased dust scale height (αz = 10−4 ) displays a more diffuse intensity map and a slight decrease in intensity perpendicular to the axis of inclination. A difference in ring thickness depending on the azimuthal location is not visible in the observed system. For the weakest vertical mixing (αz = 10−6 ) the dust substructures appear completely flat. Ring 2 is much fainter than the observed intensity and ring 1 appears significantly thinner. The model with αz = 10−5 comes closer to the observed ring thickness while having a mostly azimuthally constant ring structure. Assuming dust trapping in the L5 point of the observation we can propose potential coordinates for a yet undetected planet. Comparing the results of model fid_dres with the ALMA image, the planet offset relative to the disk center is δRA ≈ Article number, page 13 of 19 A&A proofs: manuscript no. main Fig. 13: Comparison of the synthetic images based on the models fid_dres (first and second column), p1m6fb_dres (third and fourth column) with the observation taken from Andrews et al. (2018); Isella et al. (2018). The ellipsis at the lower right of each panel visualizes the synthesized beam. The beam size is 49.7 × 41.4 mas for the synthetic images. The rms noise reaches ≈ 50 µas. The synthetic images are projected with an inclination of 46.7° and a position angle of 133.33°. Fig. 14: Comparison of the synthetic images based on the model fid_dres after 2000 orbits at 48 au with respect to the vertical dust mixing parametrized by αz . −0.352 arcsec and δDEC ≈ 0.104 arcsec. This corresponds to the coordinates RA=17h56m21.2563s, DEC=-21d57m22.3795s. 4. Discussion In the following parts we compare our results with previous works on the HD 163296 system and equivalent simulations as well as limits and caveats of the models presented here. 4.1. Comparison to previous works Studying the observed gap widths Isella et al. (2016) postulated a range of 0.5 Mjup to 2 Mjup for planet 1 at 48 au, 0.05 Mjup to 0.3 Mjup for planet 2 at 83 au and 0.15 Mjup to 0.5 Mjup for planet 3 at 137 au. With more detailed hydrodynamical models Article number, page 14 of 19 by Liu et al. (2018) using a multi-fluid dust approach the planet masses were constrained to 0.46 Mjup , 0.46 Mjup and 0.58 Mjup for the three planets respectively. At this point no asymmetries were observationally resolved. Among the publication of the DSHARP survey Zhang et al. (2018) performed an extensive parameter study with hydrodynamical planet-disk interaction simulations using a Lagrangian particle dust formalism. Their results indicate planet masses of 0.35 Mjup , 1.07 Mjup and 0.07 Mjup if a radially constant α viscosity of 10−4 is assumed. The predicted masses of 1 Mjup and 1.3 Mjup by Teague et al. (2018) for the two outer planets exceed the hydrodynamical results. However, with the uncertainties of about 50% the planet masses used in the models can be consistent with the kinematical detections. Our models indicate that for fiducial model parameters, e.g. an P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 aspect ratio of 0.05 and a radially increasing α viscosity similar to Liu et al. (2018), a minimum mass of ≈ 0.5 Mjup for planet 1 is necessary to produce a stable dust trap in the trailing Lagrange point L5. For higher masses, the amount of dust trapped in the crescent-shaped asymmetry can be slightly decreased (see A). The initial gas surface density at 48 au of Σg,0 = 37.4g/cm2 for the high mass model assuming a local dust-to-gas mass ratio of ≈ 2.4 · 10−3 is close to the findings of Zhang et al. (2018) with Σg,0 = 3 − 30g/cm2 . Isella et al. (2016) used a value of Σg,0 = 10g/cm2 . Given the proximity of the crescent-shaped asymmetry and ring 1 in the observations, it is a natural choice to normalize the dust density to the values derived from the optical depth comparison. Marzari & Scholl (1998) found that if planetesimals are small enough to be affected by the gas drag, the stability of the L4 point is reduced and the density distribution of L4 and L5 becomes asymmetric. A similar effect was observed in the gas by Masset (2002) if viscosity is included. In this case, compared to the gas drag affecting the dust, the viscous gas drag acts as the effect causing the asymmetric gas distribution. Similar results in the context of hydrodynamical simulations including gas and dust were found by Lyra et al. (2009). Recently, Montesinos et al. (2020) presented hydrodynamical simulations including multispecies particle dust exploring the stability of the L4 and L5 in the presence of a massive planet with at least one Jupiter mass. Their findings basically agree with the results presented in this paper without the effect of dust feedback. They state, that L5 captures a larger amount of dust compared to the L4 point. They argue that colder disks allow for more efficient dust trapping in these Lagrange points, lower viscosity leads to a more symmetric distribution of dust in L4 and L5 and dust entering the co-orbital region from the outer part of the disks seems to not significantly contribute to the mass of the clumps in L4 and L5. Interestingly, the Trojans populating L4 and L5 around Jupiter seem to be more numerous around the L4 point (Yoshida & Nakamura, 2005). In our models this effect only appears if dust feedback plays a significant role in this region. 4.2. Model assumptions Dust opacities are highly sensitive to its material composition and spatial structure. Estimating the surface densities from the optical depth is thus subject to a significant uncertainty. This is amplified by the choice of the dust size distribution and dust size limits. However, the features of interest, i.e. the crescent-shaped feature and ring 1 match the observations reasonably well with the high mass model, setting amin = 0.19 mm and amax = 19 mm with the MRN size distribution. The life time of the crescent-shaped asymmetry in L5 depends on diffusive processes like the turbulent viscosity and mixing of dust grains. In the resolution study (see appendix B) the resulting life time seems not to be limited within the simulated time frame. In lower resolution studies investigated in this paper the numeric diffusion artificially truncates the feature’s life time. Even by employing highly resolved simulations the age estimation of the crescent-shaped feature and thus approximately the planet itself would be difficult with the degenerate parameter space. Eccentricity and viscosity both shorten the time scale of dispersal significantly. More detailed studies and observations are necessary to constrain dynamical age of the substructures which need to be performed at higher spacial resolution. In the observation presented in Isella et al. (2018) the crescentshaped asymmetry is located at r = 55 au instead of r = 48 au. No combination of parameters in our models are able to reproduce this effect. An eccentric planet would be an intuitive choice but only leads to a disruption of the crescent-shaped feature. Dust feedback can lead to an unstable feature, ultimately leading to small clumps. In the earlier stages, dust feedback promotes dust trapping in the L4 point. No such effect is seen in the observations. Another explanation of the positioning of the crescent-shaped asymmetry could be planet migration. Depending on the migration direction and speed, the locations of the rings and features in the co-rotation region can be asymmetrically shifted in radial direction (Meru et al., 2019; Pérez et al., 2019; Weber et al., 2019). Additionally, sudden migration jumps in a system of multiple planets can temporarily create trailing asymmetries with respect to the migrating planet as shown in Rometsch et al. (submitted). It should be noted that dust coagulation and fragmentation is not considered here. More sophisticated models including these effects as shown in Drażkowska ˛ et al. (2019) could be used in this case but are computationally demanding. Ring 2 is slightly fainter in our models compared to the observations. The amount of dust that can be trapped in ring 2 depends on planet 3 since it truncates the dust flow from the outer part of the disk. One hypothesis might be that planet 3 formed later than planet 2 and thus allowed a larger amount of dust to be accumulated in the second ring. As shown in Fig. 4 it is furthermore possible to confine the range of permissible values of α by the disappearance of ring 2 for large viscosities (α > 2 · 10−3 ) and low viscosities (α < 5 · 10−4 ) due to vortex activity. The synthetic images in Fig. 13 display an additional gap in the inner dust disk. This secondary is caused by the interaction with the spiral wakes originating from planet 1. The effect is mostly visible for large Stokes numbers and dust sizes as well as high planet masses. However, a close to Jupiter mass planet is necessary to trap the needed amount of dust in the L5 point to be comparable to the observations. Results of Miranda & Rafikov (2019) indicate that radiative effects are important, even at large distances of the central star, since locally isothermal models over-pronounce the effect of the spiral wakes and secondary gaps. The same effect was shown to be important for the inner gas disk of HD 163296 in the work of Ziampras et al. (2020). It can be expected that the additional secondary ring in the inner disk disappears when radiative effects are taken into account. Nevertheless, the inner dust disk is not the main aspect of our work and the locally isothermal approach can be considered to be sufficient for modeling the crescent-shaped feature. The planet growth time scale has only a minor impact on the overall dust substructure emerging in the simulations. Differences are likely caused by the change in dust content and local dust size distribution due to dust drift. Longer growth time scales lead to a lower intensity of ring 2 due to the lack of material that has already drifted inwards before being trapped by the outer planets. The dynamical structure, especially the shape and location of the crescent-shaped asymmetry, is basically unaffected within the explored parameter space of growth time scales. The synthetic images based on the low mass model show narrower rings and differ significantly from the observations. The high mass model is thus favored in this study. 5. Conclusion We presented a parameter study of the crescent-shaped feature of the protoplanetary disk around HD 163296 using multi-fluid hydrodynamical simulations with the FARGO3D code. The model includes eight dust fluids with initial Stokes numbers ranging Article number, page 15 of 19 A&A proofs: manuscript no. main from St = 1.3 · 10−3 to St = 1.3 · 10−1 and grain sizes of amin = 0.19 mm and amax = 19 mm for the high mass model. Additionally, synthetic ALMA observations based on radiative transfer models of the hydrodynamical outputs are presented. Comparing the model with the observation, the results match qualitatively. In this work we showed that the observation of the crescentshaped feature puts important constraint on the disk and planet parameters – always under the assumption that the feature is truly caused by dust accumulation in the planet’s trailing Lagrange point L5. Most importantly, it confines the level of viscosity and planetary mass. The main findings can be summarized as follows: 1. The observed crescent-shaped asymmetry in the observation (Isella et al., 2018) can be reproduced with a Jupiter mass planet in the respective gap location at 48 au. Dust is effectively trapped in the trailing Lagrange point L5. In the case of negligible dust feedback the L4 point is not sufficiently populated to be observable. The peak of the asymmetric dust density distribution shifts towards the planet location for larger Stokes numbers and grain sizes. 2. Rescaling the dust densities to the observed optical depth of ring 1 at 67 au dust masses of 10 to 15 earth masses can be trapped in a crescent shaped feature located at the L5 point. The trapped dust mass is relatively insensitive to the choice of viscosity, aspect ratio, planet mass and eccentricity as well as the planet growth time scale. 3. Including the dust back reaction onto the gas can lead to dust trapping preferably at the leading Lagrange point L4 for initial Stoke numbers of St ≤ 3.6 · 10−2 and at later stages to fragmentation of the crescent-shaped asymmetry near the L5 point. 4. Diffusive and disruptive effects counter the stability of the dust trap in L5. Values of α ≥ 2 · 10−3 prevent the formation of an asymmetric and stable feature. Introducing eccentricity leads to the same result. The shifted location of the observed crescent-shaped feature at 55 au is not justified by an eccentric planet carving the corresponding gap in the given parameter space. 5. If the L5 feature is caused by an embedded planet, the models allow an estimation of the azimuthal planet position in the gap. The planet offset relative to the disk center is δRA ≈ −0.352 arcsec and δDEC ≈ 0.104 arcsec which corresponds to the coordinates RA=17h56m21.2563s, DEC=-21d57m22.3795s. We can thus conclude that a combination of ≈ 1Mjup and ≈ 0.5 Mjup for the inner planets in combination with a MRN dust size distribution with amin = 0.19 mm and amax = 19 mm as well as a local value of α = 2 · 10−4 can reproduce the observed crescent-shaped asymmetry and ring structures sufficiently well. The dust-to-gas ratio in the models may be overestimated since none of the features emerging in the simulations including feedback, e.g. two crescent-shaped asymmetries and fragmentation, are present in the observation. Additional high resolution studies are necessary to constrain the parameter space further, also in regard to the long-term stability of the feature. Acknowledgements. Authors Rodenkirch, Rometsch, Dullemond and Kley acknowledge funding from the DFG research group FOR 2634 ”Planet Formation Witnesses and Probes: Transition Disks” under grant DU 414/23-1 and KL 650/29-1, 650/30-1. The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 638596; P.W.). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the Article number, page 16 of 19 state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 37/935-1 FUGG. Plots in this paper were made with the Python library matplotlib (Hunter, 2007). References ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3 Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 Ayliffe, B. A., Laibe, G., Price, D. J., & Bate, M. R. 2012, MNRAS, 423, 1450 Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201 Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 Baruteau, C. & Zhu, Z. 2016, MNRAS, 458, 3927 Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25 Benítez-Llambay, P. & Masset, F. S. 2016, ApJS, 223, 11 Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 Bjorkman, J. E. & Wood, K. 2001, ApJ, 554, 615 Calvet, N., D’Alessio, P., Watson, D. M., et al. 2005, ApJ, 630, L185 Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161 de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73 Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93 Drażkowska, ˛ J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 Dullemond, C., Isella, A., Andrews, S., Skobleva, I., & Dzyurkevich, N. 2019, arXiv e-prints, arXiv:1911.12434 Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multipurpose radiative transfer tool Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150 Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99 Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68 Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16 Fouchet, L., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037 Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1 Gammie, C. F. 1996, ApJ, 457, 355 Goldreich, P. & Tremaine, S. 1979, ApJ, 233, 857 Goldreich, P. & Tremaine, S. 1980, ApJ, 241, 425 Hallam, P. D. & Paardekooper, S. J. 2020, MNRAS, 491, 5759 Hammer, M., Kratter, K. M., & Lin, M.-K. 2017, MNRAS, 466, 3533 Hammer, M., Pinilla, P., Kratter, K. M., & Lin, M.-K. 2019, MNRAS, 482, 3609 Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42 Hughes, A. M., Andrews, S. M., Espaillat, C., et al. 2009, ApJ, 698, 131 Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90 Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49 Kanagawa, K. D., Muto, T., Okuzumi, S., et al. 2018, ApJ, 868, 48 Klahr, H. H. & Bodenheimer, P. 2003, ApJ, 582, 869 Kley, W. 1999, MNRAS, 303, 696 Kley, W. & Nelson, R. P. 2012, ARA&A, 50, 211 Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 Lin, D. N. C. & Papaloizou, J. 1986, ApJ, 309, 846 Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125 Maddison, S. T., Fouchet, L., & Gonzalez, J. F. 2007, Ap&SS, 311, 3 Manger, N. & Klahr, H. 2018, MNRAS, 480, 2125 Marzari, F. & Scholl, H. 1998, Icarus, 131, 41 Masset, F. 2000, A&AS, 141, 165 Masset, F. S. 2002, A&A, 387, 605 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, Astronomical Society of the Pacific Conference Series, Vol. 376, CASA Architecture and Applications, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127 McNally, C. P., Nelson, R. P., Paardekooper, S.-J., & Benítez-Llambay, P. 2019, MNRAS, 484, 728 Meru, F., Rosotti, G. P., Booth, R. A., Nazari, P., & Clarke, C. J. 2019, MNRAS, 482, 3678 Miranda, R., Lai, D., & Méheut, H. 2016, MNRAS, 457, 1944 Miranda, R. & Rafikov, R. R. 2019, ApJ, 878, L9 Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, arXiv eprints, arXiv:2009.10768 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 Paardekooper, S. J. & Mellema, G. 2004, A&A, 425, L9 Paardekooper, S. J. & Mellema, G. 2006, A&A, 453, 1129 Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15 Picogna, G. & Kley, W. 2015, A&A, 584, A110 P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 Pinilla, P., Flock, M., Ovelar, M. d. J., & Birnstiel, T. 2016, A&A, 596, A81 Pinte, C., Price, D. J., Ménard, F., et al. 2020, The Astrophysical Journal Letters, 890, L9 Pohl, A., Pinilla, P., Benisty, M., et al. 2015, MNRAS, 453, 1768 Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543 Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2006, MNRAS, 372, L9 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33 Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12 Toomre, A. 1964, ApJ, 139, 1217 van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 Weber, P., Benítez-Llambay, P., Gressel, O., Krapp, L., & Pessah, M. E. 2018, ApJ, 854, 153 Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 Yoshida, F. & Nakamura, T. 2005, AJ, 130, 2900 Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588 Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47 Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122 Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18 Ziampras, A., Kley, W., & Dullemond, C. P. 2020, arXiv e-prints, arXiv:2003.02298 Article number, page 17 of 19 A&A proofs: manuscript no. main St = 1.3 10 1 M = 0.42 M = 0.54 M = 0.66 M = 0.78 M = 0.90 0.04 0.03 0.02 0.01 100 0 100 0.15 100 0.10 0 0.05 100 0.5 0.4 0.3 0.2 0.1 100 0 100 100 0 100 100 0 100 100 0 100 100 0 au 100 100 0 100 100 0 Surface density [g/cm2] St = 1.9 10 2 St = 5.0 10 3 M = 0.30 100 Fig. A.1: Dust surface density maps for a subset of three fluids with varying values of the planet 2 mass in Mjup at 500 orbits at r0 . The white crosses mark the position of planet 1. fiducial model ] 6.9 · 10 −2 Article number, page 18 of 19 1.9 · 10 −2 fiducial model - 2x resolution Mearth ] Mna 3.6 · 10 −2 9.6 · 10 −3 5.0 · 10 −3 2.6 · 10 −3 [ In the model fid_dres the resolution is doubled to 1120 × 1790 cells in radial and azimuthal direction respectively compared to model fid. Fig. B.1 indicates that an increased resolution leads to a more stable asymmetry in the Lagrange point L5. No decline in mass can be observed in the simulated time frame (1000 orbits at 48 au). The feature in the low resolution model fid however depletes rapidly after 600 orbits. Since the dust trapped in this feature is sensitive to diffusive and disruptive effects, like e.g. viscosity, eccentricity and the passing of the outer planet, the accelerated dispersal may be attributed to numerical diffusion. Quantifying the absolute value of the numerical diffusion is complex, however the order of magnitude can be estimated by a simple comparison of the high resolution run alpha3_dres with α = 5 · 10−4 and the fiducial model. As shown in Fig. 10 the feature lifetime is approximately comparable to the one of the lower resolution run fid with a local α = 2·10−4 . Therefore, the effect of the numerical diffusion in the low resolution model should be approximately equivalent to α ≈ 5 · 10−4 . Since FARGO3D is second-order accurate in space and the error of the dust module has been found to be proportional to a power law with an exponent of -2.2 as a function of the number of grid cells (Benítez- Mna Appendix B: Resolution study 1.3 · 10 −1 Initial Stokes Number 5 4 3 2 1 0 4 3 2 1 00 [ In Fig. A.1 a parameter study of the planet 2 mass influence is shown, involving the models p2m1 to p2m6. The nonaxisymmetric feature at the L5 point is sensitive to the mass of planet 2. In general, the passing of the planet acts as a perturber, inhibiting an effective dust trap in L5. The effect is visible in Fig. A.1 for smaller grain sizes. Lower planet masses weaken the dust trapping in ring 2. We identify the balance between effective trapping in ring 2 and optimal dust trapping in the crescent-shaped feature to be on the order of ≈ 0.5Mjup , thus the choice of the fiducial model parameter. Llambay et al., 2019), the numerical diffusion is expected to be equivalent to α ≈ 1·10−4 in the model fid_dres. The resolution is thus sufficient to describe the effect of prescribed local viscosity of α = 2 · 10−4 . Nevertheless, the absolute values of the amount of trapped dust mass in the stable phase is not significantly affected by the low resolution effect and lower resolution models are thus acceptable to quantify these values. Mearth Appendix A: Secondary planet mass 1.3 · 10 −3 200 400 Orbits 600 800 1000 Fig. B.1: Development of the trapped dust mass in the L5 region of the simulations fid and fid_dres over time. In the lower panel the grid resolution is doubled in both the radial and azimuthal direction. The dispersal of the features are prolonged in the high resolution setup. Appendix C: Planet growth time scale In Fig. C.1 simulated ALMA observations are shown for planet mass growth time scales ranging from T G = 10 T 0 to 500 T 0 . Ring 2 becomes less massive for longer planet growth time scales since dust drift depletes the outer regions before the planets reach a sufficiently high mass for efficient trapping. P. J. Rodenkirch et al.: nonaxisymmetric structures in HD 163296 Fig. C.1: Comparison of the synthetic images of different planet growth time scales (a): T G = 10 T 0 , (b): T G = 100 T 0 , (c): T G = 500 T 0 . The snapshots are taken after 2000 orbits at 48 au. Appendix D: Dust temperatures TG = 10 T0 2 2 0 4 Mna [Mearth] TG = 100 T0 2 0 4 1.3 6.9 3.6 1.9 9.6 5.0 2.6 1.3 10 10 10 10 10 10 10 10 1 2 2 2 3 3 3 Mna [Mearth] TG = 500 T0 3.6 · 10 −2 1.9 · 10 −2 9.6 · 10 −3 5.0 · 10 −3 2.6 · 10 −3 101 250 500 750 1000 1250 1500 1750 Orbits (at 48 au) Fig. C.2: Evolution of the trapped dust mass Mna in the asymmetry around the L5 point for different planet growth time scales TG. 1.3 · 10 −1 6.9 · 10 −2 102 101 2 0 0 Gas 3 Initial Stokes Number Mna [Mearth] TG = 50 T0 Initial Stokes Number 0 4 In Fig. D.1 dust temperatures from the radiative transfer calculation of the fiducial model and the prescribed gas temperatures are shown. While gas temperature gradient is smaller than the one for the dust, the temperatures and the slope of both match well at the location of planet 1, the primary region of interest where the crescent-shaped asymmetry forms. For the large grain sizes studied in this paper, a gray body approximation for the temperature is approximately valid. We thus see no significant increase in temperature comparing the grain sizes to each other. Temperature [K] Mna [Mearth] 4 r [au] 102 1.3 · 10 −3 Fig. D.1: Comparison of prescribed gas temperature in the hydrodynamical simulation and the dust temperatures of the thermal Monte-Carlo calculation at the disk mid plane for all grain sizes. The inner disk structure including the crescent-shaped asymmetry remains unaffected by the choice of parameters. Fig. C.2 reveals the temporal evolution of the dust content in the asymmetry around the L5 point for all four growth time scales. For all runs grains smaller than ≈ 1 mm become depleted after more than 1000 orbits of evolution. After initial jumps in dust mass all simulations reach a stable stationary state considering millimeter grain sizes and above. Article number, page 19 of 19