Issue |
A&A
Volume 647, March 2021
|
|
---|---|---|
Article Number | A174 | |
Number of page(s) | 18 | |
Section | Planets and planetary systems | |
DOI | https://doi.org/10.1051/0004-6361/202038484 | |
Published online | 29 March 2021 |
Modeling the nonaxisymmetric structure in the HD 163296 disk with planet-disk interaction
1
Institute for Theoretical Astrophysics, Zentrum für Astronomie, Heidelberg University,
Albert-Ueberle-Str. 2,
69120
Heidelberg, Germany
e-mail: p.rodenkirch@gmail.com
2
Institut für Astronomie und Astrophysik, Universität Tübingen,
Auf der Morgenstelle 10,
72076
Tübingen, Germany
3
Niels Bohr International Academy, The Niels Bohr Institute, University of Copenhagen,
Blegdamsvej 17,
2100
Copenhagen Ø, Denmark
4
Departamento de Astronomía, Universidad de Chile,
Casilla 36-D,
Santiago, Chile
5
Departamento de Física, Universidad de Santiago de Chile,
Av. Ecuador 3493,
Estación Central,
Santiago, Chile
Received:
24
May
2020
Accepted:
14
December
2020
Context. High-resolution ALMA observations such as the DSHARP campaign have revealed a variety of rich substructures in numerous protoplanetary disks. These structures consist of rings, gaps, and asymmetric features. It has been debated whether planets can be accounted for among 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 that the formation of the crescent-shaped asymmetry in the HD 163296 disk can be attributed to planet-disk interaction. The goal is to obtain constraints on planet masses, eccentricities, and disk viscosities. Furthermore, we test the reproducibility of the two prominent rings in the HD 163296 disk at 67 and 100 au.
Methods. We performed two-dimensional, multi-fluid, hydrodynamical simulations with the FARGO3D code, including three embedded planets in the setup. Dust is described via the pressureless fluid approach and distributed over eight size bins. The resulting grids were post-processed with the radiative transfer code RADMC-3D and CASA software to model the 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 in the trailing Lagrange point L5, preferably, with a mass of between 10 and 15 earth masses. The observation of such a feature constrains the level of viscosity and planetary mass. Increased values of eccentricity of the innermost Jupiter mass planet negatively impacts 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 for the existence of such a feature. Including dust feedback in 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 on the part 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 / hydrodynamics / radiative transfer
© ESO 2021
1 Introduction
The advent of high angular resolution millimeter continuum observations with the Atacama Large Millimeter and Submillimeter Array (ALMA) has made it possible to gain insights into the substructures of protoplanetary disks. 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 2015). An extensive survey with the goal to image detailed structures in 20 disks was performed with the DSHARP campaign (Andrews et al. 2018) at an 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 such as spirals and crescent-shaped features.
Currently, a matter of debate is whether these structures are signposts of embedded planets (Zhang et al. 2018). Planet-disk interaction has been a central topic in discussing 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 drifting 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 planetesimal formation (Ayliffe et al. 2012).
Dong et al. (2015) found that multiple planets hold clues to explaining 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 via 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, for example 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 what is referred to, driving angular momentum transport and thus accretion due the underlying turbulent processes. Therefore, 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, such as 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, such as 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 planets embedded in protoplanetary disks. A popular and well-studied disk is the one around the Herbig Ae star HD 163296 at a distance of 101 pc (Gaia Collaboration 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 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 the crescent-shaped asymmetry within the inner gap located at 48 au (Huang et al. 2018). The feature itself is situated at a radial distance of 55 au, that is, 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, however, a preliminary model was presented in Zhang et al. (2018), based on planet-disk interaction. In these models, asymmetries in the co-orbital region are common if the viscosity is low.
Prior to 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. Thanks to synthetic images using radiative transfer calculations, it was possible to 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.
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 yet been subject to more detailed analysis. 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, such as planet mass, turbulent viscosity, and dust size. Similarly to the study by Zhang et al. (2018), we employ two-dimensional hydrodynamical models with a fluid formulation of dust.
Section 2 introduces the physical model and code setup as well as the post processing pipeline to predict observable features. In Sect. 3, we present the main results of our study. Section 4 compares these findings with previous works and addresses limitations of the model. In Sect. 5, we summarize the main results and give our 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) and 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:
Here, Σg denotes the gas surface density, Σd,i the corresponding dust species, the gas pressure which is 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 vg and vd,i the gas and dust velocities respectively. Dust feedback is included by the term ∑iΣifi 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: (5)
with the diffusion constant Di being proportionalto the turbulent viscosity ν (Youdin & Lithwick 2007): (6)
The Stokes number Sti of dust species i is proportional to the stopping time tstop and can be written as (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 (8)
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 refer to the three planets as planet 1, planet 2, and planet 3. The same notation is used for the apparent ring structures, namely, ring 1: the observed or modeled ring at 67 au; and 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 were set, with the planet mass, Mp, adjusted 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: (9)
TG refers to the planet growth time scale ranging from 10T0 to 500T0 with T0 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, thus, we 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: (10)
with r0 = 48 au and the initial gas and dust surface densities Σg,0 and Σd,0. Similarly to the models of Liu et al. (2018), we chose a surface density slope of p = 0.8. For the dust, a sharper cutoff of s = 2 was chosen as compared to the gas cutoff of s = 1. In all simulation runs, an initial dust-to-gas ratio of is assumed.
We approximate the disk thermodynamics with a locally isothermal equation of state. The model is parametrized through locally isothermal sound speed: (11)
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 (12)
with the proton mass mp. The temperature at 48 au matches the findings of Dullemond et al. (2020).
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): (13)
where the parameters R and σ are set to 144 au and 1.25, respectively. Similarly to the values in Liu et al. (2018), the parameter R refers to the mid-point of the transition in α, whereas σ defines the slope.
Fig. 1 Sketch of the oberserved dust substructure of the HD 163296 system. The colored rings and the crescent-shaped feature mimic the dust whereas 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. |
2.3 Boundary conditions
The radial velocities in the ghost cells are set according to anti-symmetric 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 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.
2.4 Code setup and parameters
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 eight 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 cm−2 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 neglect 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 T0. Simulation runs with a nonzero growth time scale TG and the model fid_dres are run until 2000 T0. 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 Sect. 3 for more details.
2.5 Radiative transfer model and post-processing
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 diffusion model of Dubrulle et al. (1995): (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 the polar direction, the grid is expanded by 32 cells which are equally spaced up to zlim = π∕2 ± 0.3. The verticaldisk density profile is assumed to be isothermal and the conversion from the surface density to the local volume density is calculated as follows: (15)
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−1. 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 CASA-5.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 two hours whereas the more compact configuration is integrated over a reduced time with a factor of 0.22 corresponding 0.44 h.
We employed the same cleaning procedure as made available 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.
Relevant simulation runs and their respective numerical parameters.
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 T0. Lower panel:prescribed α viscosity. The blue line visualizes the radially increasing α viscosity setin the model. |
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. Figures 3 and 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 closely.
Additionally, if a crescent-shaped asymmetry is present, it is situated in the gap caused by the innermost 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 theaspect 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 whereas a second similar asymmetry at the L4 point appears for values of H∕r < 0.05. Thecrescent weakens for smaller grain sizes. Also the second prominent ring beyond planet 2 at 83 au clearly weakens for larger aspect ratios. In the case of a 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, similarly to the large aspect ratio in Fig. 3. The collection of these results also stresses the problem of a radially constant α viscosity with respect to the observed HD 163296 system. To reproduce a nonaxisymmetric feature in the vicinity of the inner planet and a smooth ring-shaped 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.
3.2 Rings
Before quantifying the nonaxisymmetric feature in the model, we want to compare properties of the ring structures with previous works and observations. The general procedure is to azimuthally average the dust density maps after 1000 orbits at 48 au and to invoke a 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.
Figure 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 work excludes contributions from the prominent nonaxisymmetric structures. Several properties become apparent. Ring 1 is wider than in the simulated profile and the peak location of ring 1 is located further outward with respect to the simulated one. Furthermore, 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 actual dust-to-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.
One conclusion drawn by Dullemond et al. (2018) was that the optical depth observed in the DSHARP survey is remarkably close to unity and that the rings were optically thin. Later Liu (2019) and 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 couldbe thus larger than expected.
For ringfitting we use a Gaussian: (16)
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.
Fig. 3 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 dustaccumulation in the co-orbital regions of the planets. The white crosses mark the position of planet 1. |
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. |
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. |
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. |
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. 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 achange 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 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 assuming 1 mm grains with the same 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.
Optical depth fitting results for ring 1.
Gaussian fit results of the simulated dust rings for all Stokes numbers, St.
3.4 Secondary planet mass
In the models p2m1 - p2m6, the mass of the secondary planet at 83 au is varied 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 2since 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 Figs. 3 and 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 to constrain physical properties of the dynamical system.
3.5.1 Structure and 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 smaller Stokes numbers.
We consider whether the dynamics of this feature change with the consideration of a dust back-reaction onto the gas. An exampleof the impact from dust feedback (model p1m1fb_dres) is given in Fig. 7 on the right hand side, with the same parameters as in 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 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 Sect. 3.2.
In Fig. 9 we display the results of model p1m6fb_dres. In contrast 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 feedback.
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.
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. |
3.5.2 Dynamical stability
In principle, the crescent-shaped features in the co-orbital region are subject to diffusive processes such as dust diffusion due to turbulent mixing or gravitational interaction with the planetary system (for example eccentric orbits). The dust trapping mechanism has to counteract these disruptive forces for the feature to be dynamically stable.
Figure 10 corroborates this line of argument. The stability of the L5 feature is sensitive to the local value of the α viscosity. Given 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 discernible 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. Therefore, in Fig. 10, we plotted results with 14 cells per scale height.
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. |
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 is preferablyconcentrated in L4. |
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 thenonaxisymmetric 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. |
3.5.3 Dust mass
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 the 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 T0 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 in the mass of planet 2 has a small influence on the 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 TG = 10, 50, 100 and 500 orbits were performed. Qualitatively, the results are mostly unaffected by the choice of TG. 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 whereas 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 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 co-orbital 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.
Fig. 11 Trapped dust masses in earth masses in the L5 feature for the largest three initial Stokes numbers. The mass values are averaged over the approximately constant mass time frame starting from 100 T0. 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. |
3.6 Synthetic images
The main question remains: whether dust trapping in the L5 point in the models presented could explain the observed feature in HD 163296. We apply the procedure described in Sect. 2.5 and thus extend the surface density maps to three dimensional grids, then we 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-resolution 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 crescent-shaped feature is the more elongated shape compared to the fiducial model. Furthermore, the models produce a radially symmetrically located feature whereas 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 Sect. 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. In 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 that 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 ≈−0.352 arcsec and δDec ≈ 0.104 arcsec. This corresponds to the coordinates RA=17h56m21.2563s, Dec= −21d57m22.3795s.
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. |
Fig. 13 Comparison of the synthetic images based on the models fid_dres (first and second columns), p1m6fb_dres (third and fourth columns) 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. |
4 Discussion
In the following, we compare our results with previous works on the HD 163296 system and equivalent simulations as well as the limits and caveats of the models presented here.
4.1 Comparison to previous works
Upon 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 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, for example an 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 1is 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 Appendix A).
The initial gas surface density at 48 au of Σg,0 = 37.4 g cm−2 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−30 g cm−2. Isella et al. (2016) used a value of Σg,0 = 10 g cm−2. 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 that include multi-species particle dust and explore 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 the relevant 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, such as 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 lifetime of the crescent-shaped asymmetry in L5 depends on diffusive processes, such as 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 observationsare 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 crescent-shaped 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. (2020).
It should be noted that dust coagulation and fragmentation is not considered here. More sophisticated models including these effects as shown in Drąż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, it is necessary to have 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. The 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, in particular, the shape and location of the crescent-shaped asymmetry, is basically unaffected within the explored parameter space of growth time scales. The synthetic images that are based on the low-mass model show narrower rings and differ significantly from the observations. Therefore, the high-mass model is favored in this study.
5 Conclusion
In this work, we present 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 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 are a qualitative match.
We show that the observation of the crescent-shaped feature puts an 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 constrains 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 timescale.
- 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 as well as, at later stages, to the fragmentation of the crescent-shaped asymmetry near the L5 point.
- 4.
Diffusiveand 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 L5feature 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 that include feedback, for example 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
P.J.R., T.R., C.P.D. and W.K. 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 program (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 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).
Appendix A Secondary planet mass
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.
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. |
Appendix B Resolution study
In the modelfid_dres, the resolution is doubled to 1120 × 1790 cells in radial and azimuthal direction respectively compared to model fid. Figure 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.
Fig. B.1 Development of the trapped dust mass in the L5 region of the simulations fid and fid_dres over time. 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. |
Since the dust trapped in this feature is sensitive to diffusive and disruptive effects, such as 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-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.
Appendix C Planet growth time scale
In Fig. C.1, simulated ALMA observations are shown for planet mass growth time scales ranging from TG = 10 T0 to 500 T0. 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.
Fig. C.1 Comparison of the synthetic images of different planet growth time scales (a): TG = 10 T0, (b): TG = 100 T0, (c): TG = 500 T0. The snapshots are taken after 2000 orbits at 48 au. |
The inner disk structure including the crescent-shaped asymmetry remains unaffected by the choice of parameters. Figure 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.
Fig. C.2 Evolution of the trapped dust mass Mna in the asymmetry around the L5 point for different planet growth time scales TG. |
Appendix D Dust temperatures
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 temperatureis approximately valid. We thus see no significant increase in temperature comparing the grain sizes to each other.
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. |
References
- ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
- Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [Google Scholar]
- Ayliffe, B. A., Laibe, G., Price, D. J., & Bate, M. R. 2012, MNRAS, 423, 1450 [NASA ADS] [CrossRef] [Google Scholar]
- Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201 [NASA ADS] [CrossRef] [Google Scholar]
- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
- Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927 [NASA ADS] [CrossRef] [Google Scholar]
- Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [NASA ADS] [CrossRef] [Google Scholar]
- Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25 [Google Scholar]
- Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [NASA ADS] [CrossRef] [Google Scholar]
- Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
- Calvet, N., D’Alessio, P., Watson, D. M., et al. 2005, ApJ, 630, L185 [NASA ADS] [CrossRef] [Google Scholar]
- Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
- Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73 [NASA ADS] [CrossRef] [Google Scholar]
- Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93 [NASA ADS] [CrossRef] [Google Scholar]
- Drążkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91 [NASA ADS] [CrossRef] [Google Scholar]
- Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
- Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool [Google Scholar]
- Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Dullemond, C., Isella, A., Andrews, S., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137 [EDP Sciences] [Google Scholar]
- Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99 [Google Scholar]
- Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150 [NASA ADS] [CrossRef] [Google Scholar]
- Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fouchet, L., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
- Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [Google Scholar]
- Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [Google Scholar]
- Hallam, P. D., & Paardekooper, S. J. 2020, MNRAS, 491, 5759 [Google Scholar]
- Hammer, M., Kratter, K. M., & Lin, M.-K. 2017, MNRAS, 466, 3533 [CrossRef] [Google Scholar]
- Hammer, M., Pinilla, P., Kratter, K. M., & Lin, M.-K. 2019, MNRAS, 482, 3609 [NASA ADS] [CrossRef] [Google Scholar]
- Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42 [Google Scholar]
- Hughes, A. M., Andrews, S. M., Espaillat, C., et al. 2009, ApJ, 698, 131 [NASA ADS] [CrossRef] [Google Scholar]
- Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
- Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
- Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49 [Google Scholar]
- Kanagawa, K. D., Muto, T., Okuzumi, S., et al. 2018, ApJ, 868, 48 [NASA ADS] [CrossRef] [Google Scholar]
- Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] [Google Scholar]
- Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] [Google Scholar]
- Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211 [Google Scholar]
- Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] [Google Scholar]
- Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
- Liu, H. B. 2019, ApJ, 877, L22 [NASA ADS] [CrossRef] [Google Scholar]
- Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87 [NASA ADS] [CrossRef] [Google Scholar]
- Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [Google Scholar]
- Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Maddison, S. T., Fouchet, L., & Gonzalez, J. F. 2007, Ap&SS, 311, 3 [NASA ADS] [CrossRef] [Google Scholar]
- Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125 [NASA ADS] [CrossRef] [Google Scholar]
- Marzari, F., & Scholl, H. 1998, Icarus, 131, 41 [NASA ADS] [CrossRef] [Google Scholar]
- Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Masset, F. S. 2002, A&A, 387, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
- McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127 [Google Scholar]
- McNally, C. P., Nelson, R. P., Paardekooper, S.-J., & Benítez-Llambay, P. 2019, MNRAS, 484, 728 [NASA ADS] [CrossRef] [Google Scholar]
- Meru, F., Rosotti, G. P., Booth, R. A., Nazari, P., & Clarke, C. J. 2019, MNRAS, 482, 3678 [NASA ADS] [CrossRef] [Google Scholar]
- Miranda, R., & Rafikov, R. R. 2019, ApJ, 878, L9 [NASA ADS] [CrossRef] [Google Scholar]
- Miranda, R., Lai, D., & Méheut, H. 2016, MNRAS, 457, 1944 [NASA ADS] [CrossRef] [Google Scholar]
- Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, A&A 642, A224 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Paardekooper, S. J., & Mellema, G. 2004, A&A, 425, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Paardekooper, S. J., & Mellema, G. 2006, A&A, 453, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15 [NASA ADS] [CrossRef] [Google Scholar]
- Picogna, G., & Kley, W. 2015, A&A, 584, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinilla, P., Flock, M., Ovelar, M. d. J., & Birnstiel, T. 2016, A&A, 596, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9 [NASA ADS] [CrossRef] [Google Scholar]
- Pohl, A., Pinilla, P., Benisty, M., et al. 2015, MNRAS, 453, 1768 [NASA ADS] [CrossRef] [Google Scholar]
- Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543 [NASA ADS] [CrossRef] [Google Scholar]
- Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2006, MNRAS, 372, L9 [NASA ADS] [CrossRef] [Google Scholar]
- Rometsch, T., Rodenkirch, P. J., Kley, W., & Dullemond, C. P. 2020, A&A, 643, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
- Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12 [NASA ADS] [CrossRef] [Google Scholar]
- Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
- van der Marel,N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
- Weber, P., Benítez-Llambay, P., Gressel, O., Krapp, L., & Pessah, M. E. 2018, ApJ, 854, 153 [Google Scholar]
- Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178 [NASA ADS] [CrossRef] [Google Scholar]
- Whipple, F. L. 1972, From Plasma to Planet, ed. A. Elvius (New York: Wiley Interscience Division), 211 [Google Scholar]
- Yoshida, F., & Nakamura, T. 2005, AJ, 130, 2900 [NASA ADS] [CrossRef] [Google Scholar]
- Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
- Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47 [NASA ADS] [CrossRef] [Google Scholar]
- Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122 [NASA ADS] [CrossRef] [Google Scholar]
- Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18 [NASA ADS] [CrossRef] [Google Scholar]
- Ziampras, A., Kley, W., & Dullemond, C. P. 2020, A&A 637, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Sketch of the oberserved dust substructure of the HD 163296 system. The colored rings and the crescent-shaped feature mimic the dust whereas 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. |
|
In the text |
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 T0. Lower panel:prescribed α viscosity. The blue line visualizes the radially increasing α viscosity setin the model. |
|
In the text |
Fig. 3 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 dustaccumulation in the co-orbital regions of the planets. The white crosses mark the position of planet 1. |
|
In the text |
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. |
|
In the text |
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. |
|
In the text |
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. |
|
In the text |
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. |
|
In the text |
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. |
|
In the text |
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 is preferablyconcentrated in L4. |
|
In the text |
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 thenonaxisymmetric 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. |
|
In the text |
Fig. 11 Trapped dust masses in earth masses in the L5 feature for the largest three initial Stokes numbers. The mass values are averaged over the approximately constant mass time frame starting from 100 T0. 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. |
|
In the text |
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. |
|
In the text |
Fig. 13 Comparison of the synthetic images based on the models fid_dres (first and second columns), p1m6fb_dres (third and fourth columns) 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°. |
|
In the text |
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. |
|
In the text |
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. |
|
In the text |
Fig. B.1 Development of the trapped dust mass in the L5 region of the simulations fid and fid_dres over time. 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. |
|
In the text |
Fig. C.1 Comparison of the synthetic images of different planet growth time scales (a): TG = 10 T0, (b): TG = 100 T0, (c): TG = 500 T0. The snapshots are taken after 2000 orbits at 48 au. |
|
In the text |
Fig. C.2 Evolution of the trapped dust mass Mna in the asymmetry around the L5 point for different planet growth time scales TG. |
|
In the text |
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. |
|
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.