© C. Bergez-Casalou et al. 2022

1 Introduction

Recent observations with the Atacama Large Millimeter/Submillimeter Array (ALMA) and the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instruments show protoplanetary disks that present different kinds of substructures (rings, gaps, cavities, and asymmetries) present in the gas (e.g., Teague et al. 2018; Pinte et al. 2020) and in the dust (e.g., ALMA Partnership 2015; Avenhaus et al. 2018; Andrews et al. 2018). These substructures may have different possible origins, including: self-induced dust traps due to dust growth and dust backreaction on the gas (Gonzalez et al. 2017), dust growth in snow lines (Zhang et al. 2015), zonal flows (Flock et al. 2015), secular gravitational instabilities (Takahashi & Inutsuka 2016; Tominaga et al. 2020), sintering-induced rings (Okuzumi et al. 2016), and gap opening embedded planets (Pinilla et al. 2012).

Focusing on the features created by planets, it is hard to observe the planets directly while they are embedded in their protoplanetary disk (Sanchis et al. 2020; Kloster & Flock 2019; Asensio-Torres et al. 2021). Therefore, analyzing the dust gap size (Zhang et al. 2018) or the CO velocity perturbations (Teague et al. 2018; Pinte et al. 2020) are ways to indirectly derive the properties of potentially embedded planets. Assuming that these features are indeed caused by planets, we are able to probe forming planets that are not observable directly and that will continue to evolve by accretion and migration in disks. These objects can then be used to derive or confirm some constraints on planet formation processes.

On the other hand, our own Solar System has some characteristics representative of its birth environment. For example, meteorites are solid remains of the protoplanetary solid disk. Kruijer et al. (2017) showed that their chemical composition in the Solar System can be used to constrain the time at which Jupiter’s core formed, as it is supposed to separate the reservoirs of carbonaceous and non-carbonaceous chondrites by blocking the pebbles flowing through the disk. Using solid mass estimates from the asteroid and Kuiper belts, Lenz et al. (2020) tried to reproduce the possible gas and solid distributions of our natal protoplanetary disk.

Different models investigate how different parts of the Solar System could have formed. For example, the classical model (Wetherill 1994; Raymond et al. 2009b) attempts to reproduce the inner Solar System via impacts and the accretion of planet embryos and planetesimals; in the Nice model (Gomes et al. 2005; Nesvorný 2011; Morbidelli et al. 2018) the dissipation of the gas disk triggers a dynamical instability, spreading the solids in the system; and in the Grand Tack scenario (Walsh et al. 2011; Pierens et al. 2014), Jupiter and Saturn migrate inward during the gas phase until their capture in resonance, inducing an outward migration of the two giants. According to all these models, different noticeable features would be created in the protoplanetary disk as the models require different giant planet configurations, which influence the dust distributions in the disk. The goal of this study is to determine what kinds of features a Solar System embedded in its natal gaseous and dusty disk could produce and how they would be observed by today’s instruments. By comparing the resulting images to current observations, we are able to situate the Solar System protoplanetary disk in the wide spectrum of planet-forming disks.

In order to derive these synthetic images, we divided our study into four different steps: we start by simulating the gaseous disk containing the different giant planets using 2D isothermal hydro-simulations with the FARGO2D1D code (Crida et al. 2007). These simulations allow us to derive a detailed radial gas profile that is used as input for a dust evolution model, derived in the TWO-POP-PY code (Birnstiel et al. 2012). This radial dust model takes growth, fragmentation, and drift into account. The resulting dust size distributions are then extended in three dimensions and used in a radiative transfer code, RADMC3D (Dullemond et al. 2012), to derive the synthetic images at different wavelengths; these images are finally convolved with different beam sizes in order to represent more realistic observations.

This paper is structured as follows. First, in Sect. 2 we outline the different numerical setups of the steps listed above. We present the results in Sect. 3: the results of the dust evolution model are presented in Sect. 3.1, and then the derived images, convolved with beams, are shown in Sect. 3.2. The derived images are discussed in Sect. 4, and we summarize and conclude our results in Sect. 5.

2 Numerical setups

In this study, we want to simulate how the dust would be distributed in a protoplanetary disk where multiple fully formedgiant planets are embedded. We therefore proceed in four steps: (i) first, hydrodynamical simulations are run with FARGO2D1D¹ in order to determine the gas distribution in the disk, considering different planet configurations, described in Sect. 2.1 and exploring different disk parameters; (ii) using the time and azimuthally averaged gas distribution from the hydrodynamical simulations, we investigate how dust behaves in such disks using a dust evolution code called TWO-POP-PY², giving us the dust surface density distributions as a function of radius and grain size; (iii) we derive the synthetic images from a 3D extension of the dust distributions using the radiative transfer code RADMC3D³, which outputs can then be convolved with different beams in step (iv), giving us realistic images of the disks.

In this section, we present the different setups for each step, starting by presenting the different planet configurations explored. In Sects. 2.2 and 2.3, we present the numerical setups we used for the hydrodynamical and dust evolution simulations. The radiative transfer code as well as the beam convolutions are presented in Sect. 2.4.

2.1 Planet configurations

As we simulate planetary systems that are still embedded in their gas disks, we consider two possible formation scenarios. Assuming that the four giant planets in the Solar System are already formed (i.e., their mass are fixed to their present value), they are placed in either a Compact or a Spread configuration (see Table 1). In the Compact configuration the planets are located in a tight configuration corresponding to the initial configuration needed by the Nice instability to occur (Gomes et al. 2005); whereas in the Spread configuration, the planets semimajor axis are increased by 30% compared to their nowadays positions. This aims to take into account migration, assuming that after they formed they migrated inward from further away orbits toward their current configuration (Bitsch et al. 2015; Sotiriadis et al. 2017; Pirani et al. 2019; Öberg & Wordsworth 2019).

To investigate planet formation in a global scale, we also chose a third configuration representing an exoplanetary system. This system is studied to help the comparison between the resulting images and the observations. It is composed of three giant planets of 1 Mj located at ~5 R_H,mut (Chambers et al. 1996) from each other in the inner region of the disk. This configuration represents a possible intermediate step of N-body simulations aimed to study giant planet formation (Bitsch et al. 2020), but serves also as initial conditions for N-body simulations aimed to explain the eccentricity distribution of giant planets via scattering (Jurić & Tremaine 2008; Raymond et al. 2009a).

As our goal is to constrain giant planet formation, we also investigate the impact of the planet masses on the disks. By changing the planet masses, we probe different stages in the formation process. We study this effect for two of our configurations, the Spread and the Three Giants, with reduced planet masses corresponding to one-half or one-third of their current mass. In return, we enhance the disk gas masses by 1% and 2%, respectively, as the disks harboring the less massive planets represent an earlier evolution step. Only one set of disk parameters is used in this case (α = 10⁻⁴ and $h=0.025\times r_{\text{AU}}^{2/7}$; see next subsection).

2.2 Hydro-dynamical setup

In order to derive the gas disk profile in presence of four giant planets, we run hydrodynamical simulations with the FARGO2D1D code (Crida et al. 2007). This code is composed of a linear 2D (r, ϕ) grid, where the planets are located, surrounded by two linear azimuthally symmetric 1D grids. It allows us to simulate thewhole viscous evolution of the disk at a reasonable computational cost. The 2D grid spans from 1 AU to 52 AU and is prolonged by the inner 1Dgrid from 0.1 to 1 AU and by the outer 1D grid from 52 to 160 AU, except whenotherwise specified. The interfaces between the 2D and 1D grids are chosen to be far enough from the planets so that we can consider the disk axisymmetric.

The diskis locally isothermal and subject to an α viscosity as described in Shakura & Sunyaev (1973). We investigate the influence of viscosity by taking three different values for α = 10⁻⁴, 10⁻³, 10⁻². For the aspect ratio, we use the minimum mass solar nebula (MMSN) profile, derived by Weidenschilling (1977); Hayashi (1981): $h=0.033\times r_{\text{AU}}^{2/7}$. We also investigate the impact of the aspect ratio by taking a lower value: $h=0.025\times r_{\text{AU}}^{2/7}$. For simplicity in this paper, the first aspect ratio will be called the “MMSN-like” aspect ratio whereas the second one is referred as the “small” aspect ratio. In each case, the snow line (i.e., the location at which water vapor freezes in the disk, T = 160 K) is located at 3.3 AU in the MMSN-like case and at 0.9 AU in the small aspect ratio one. Even if the snow lines are located within the simulation domain, we neglect any physical changes that can occur around this location.

The radial extent of our gaseous initial disk (from 0.1 to 160 AU) is consistent with Lenz et al. (2020) and Kretke et al. (2012), where they derived the possible properties of the Solar System protoplanetary disk taking into account different available constraints. We also chose a gas surface density profile in agreement with Lenz et al. (2020): $\Sigma (r)={\Sigma }_0\times {(r/r_{\text{AU}})}^{-1}$, where Σ₀ = 836.1 g cm⁻² at r = 1 AU. This value was chosen so that the total mass of the disk is 0.1 M_⊙. Even if this is considered a heavy initial disk, its large radial extent allows us to neglect self-gravity. Furthermore, as FARGO2D1D features open inner and outer boundaries, the mass of the disk will decrease with viscous evolution.

The resolution is such that the innermost planet (i.e., the planet located at r_p = 5.2 AU) is resolved by five grid cells within its Hill radius. This leads, for the 2D grid, to a radial resolution of N_r,2D = 707 cells with an azimuthal resolution of N_ϕ,2D = 454 cells. This corresponds to a radial resolution of N_r,1D = 2218 when considering the whole disk (i.e., the two 1D grids combined with the 2D grid). As at low viscosity some instabilities can be triggered (Klahr & Bodenheimer 2003; Fu et al. 2014), we enhanced the resolution for α = 10⁻⁴: N_r,2D = 1414 and N_ϕ,2D = 906 leading to N_r,1D = 4436. In one particular case, the 2D–1D boundary was too close to the outer planet, creating unrealistic features. Therefore, in the Spread configuration case with low alpha (α = 10⁻⁴) and small aspect ratio, the 2D-1D outer boundary was moved from 52 to 78 AU. The radial resolution was adapted in the 2D part to match the radial resolution of the other simulations at low viscosity: here, N_r,2D = 2135.

The planets are introduced into the disk with a mass-taper function defined as $$m_{\text{taper}}={\sin }^2(t/(4n_{\text{orb}})),$$(1)

where n_orb is the number of orbits used to growfrom 0 to m_p. Here we chose n_orb = 10. The disk is integrated for 12 500 orbits at 5.2 AU (~1.5 × 10⁵ years) until it adjusts to the perturbations induced by the giant planets. The 2D density profiles after 12 500 orbits can be seen in Appendix A for each configuration and disk parameter. The disks are then evolved for another 2500 orbits at 5.2 AU (~ 3.0 × 10⁴ yr). These gas density profiles are used as an input for the dust evolution setup.

2.3 Dust evolution setup

To derive the dust distributions in the disks, we used the dust evolution code called TWO-POP-PY (Birnstiel et al. 2012, 2015). TWO-POP-PY computes the radial motion of grains as well as their growth from an initial dust and gas radial profile. These initial profiles are derived from the hydrodynamical setup presented above. The same radial resolution is used as in the hydro-simulations (N_r,1D = 2218 for the highest viscosities and N_r,1D = 4436 for α = 10⁻⁴).

The initial gas profile corresponds to the azimuthally averaged gas density and velocity profiles averaged in time over 2500 orbits at 5.2 AU (average taken between t = 12 500 and 15 000 orbits). This time average will smooth the highly perturbed disk, which is a necessary step as we do not simulate the gas and dust evolution simultaneously in 2D, as in more sophisticated simulations (Drążkowska et al. 2019). We discuss this choice in Sect. 4.4. The radial gas density profiles can be found in Appendix A.

The initial dust profile is derived from the gas profile by assuming a dust-to-gas ratio of 0.01. As we consider that all the planets already fully formed in the disk, we need to take into consideration that some planets already reached the pebble isolation mass and therefore are able to block the pebble flux from the outer regions of the disk (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018). The pebble isolation mass can be estimated from the aspect ratio h of the disk: $$M_{\text{iso}}\simeq 20{(\frac{h}{0.05})}^3{M}_{\text{Earth}}.$$(2)

For our different configurations and each disk scale height profile, Jupiter and Saturn are always above the pebble isolation mass. In the Three-Giants configuration, each planet is above the pebble isolation mass. As some of the giants gaps are able to block the pebble flux from the outer disk, we assume that the dust located between the inner edge of the disk and the outer gap edge of the furthest planet that has reached the pebble isolation mass had time to drift inward during planet formation, making the inner disk depleted in dust. Therefore, our initial dust profile can be written as $${\Sigma }_{\text{d}}=\{\begin{array}{ll}0\hfill & \text{if}r\le r_{\text{p},\text{trunc}}+2H_{\text{p},\text{trunc}}\hfill \\ 0.01\times {\Sigma }_{\text{g}}\hfill & \text{if}r>r_{\text{p},\text{trunc}}+2H_{\text{p},\text{trunc}}\hfill \end{array},$$(3)

where r_p,trunc is the semimajor axis of the planet considered to block the dust flux (Saturn in the Solar System cases, the third giant in the Three-Giants case) and H_p,trunc is the disk gas scale height at the location of the planet, used to estimate the position of the outer edge of the planet gap (Paardekooper & Mellema 2006). The capacity of the giant planets to block the small dust flux is dependent on the dust diffusion and on the gap depth, which in turn depends on the disk viscosity. We show in Sect. 3.1 that if the viscosity is high, the gaps are actually not strong enough and dust diffusion is important, allowing dust from the outer disk to flow in and fill the inner disk (de Juan Ovelar et al. 2016). Therefore, the inner disk will remain empty only if the viscosity is low enough to block the dust flux from the outer disk.

The model is evolved for 1 Myr, during which we consider that the gas disk does not evolve significantly (i.e., the gas profile is fixed). The impact of time evolution is discussed in Sect. 4.4. During these 1 Myr, the grains grow, fragment and drift. The maximal grain size depends on each limiting mechanism: a_max = min(a_frag, a_drift, a_growth). When they are limited by fragmentation, the maximal size the grains can reach is written as (Birnstiel et al. 2012) $$a_{\text{frag}}=\frac{2{\Sigma }_{\text{g}}{u}_{\text{f}}^2}{3\pi \alpha {\rho }_{\text{s}}{c}_{\text{s}}^2},$$(4)

where Σ_g is the gas surface density, u_f is the fragmentation velocity, α is the α-viscosity parameter, ρ_s is the internal density of the grains, and c_s is the sound speed. The internal density of the grains is taken to be ρ_s = 1.675 g cm⁻³, matching the Disk Substructures at High Angular Resolution (DSHARP) survey’s composition of grains (Birnstiel et al. 2018). Regarding the fragmentation velocities, velocities of around ~ 10 m s⁻¹ are required in simulations to explain the formation of large planetesimals (Drążkowska et al. 2016; Lenz et al. 2019). However, laboratory experiments of dust collisions under physical conditions of protoplanetary disks show that these velocities cannot be that high (~ 1 m s⁻¹) (Musiolik & Wurm 2019; Schneider et al. 2019). Pinilla et al. (2021) investigate which conditions are necessary for planet formation to happen with a low fragmentation velocity. However, it is clear from Eq. (4) that a lower fragmentation velocity willyield small grain sizes in the fragmentation limit, if the disk’s viscosity is large. In order to reach at least millimeter particles observable with ALMA, we thus use larger fragmentation velocities (1, 3, and 10 m s⁻¹) for higher alpha values(10⁻⁴, 10⁻³, and 10⁻²) in order to have a similar fragmentation limit for all simulations.

The grains are also subject to radial drift: as the gas moves on a slightly sub-Keplerian orbit, the grains feel a head-wind pushing them inward in the disk (Weidenschilling 1977; Brauer et al. 2008). The drift barrier corresponds to the limiting size where the grains drift faster than they grow: $$a_{\text{drift}}=f_{\text{d}}\frac{2{\Sigma }_{\text{d}}{v}_K^2}{\pi {\rho }_{\text{s}}\left|\gamma \right|c_{\text{s}}^2},$$(5)

where f_d = 0.55 is a correcting factor from Birnstiel et al. (2012), v_K is the Keplerian velocity and |γ| = |dln(P)∕dln(r)| is the gas pressure profile of the disk, derived directly from our hydrodynamical simulations.

Finally, the maximum size of the grains is limited by their own growth timescale: $$\begin{array}{l}{t}_{\text{g}}=\frac{{\Sigma }_{\text{g}}}{e_{\text{s}}{\Omega }_K{\Sigma }_{\text{d}}}\\ a_{\text{growth}}=a_0\times \text{exp}(t/t_{\text{g}}),\end{array}$$

where e_s = 1 is the sticking probability, a₀ = 2.5 × 10⁻⁶ cm is the initial size of all the grains in the disk and t_g is the growth timescale.

Knowing the limiting sizes, the code then computes the motion of the grains by taking into account their relative velocities, their diffusion in the disk and wether they are coupled to the gas or not. This coupling can be monitored with the Stokes number. Assuming spherical grains in an Epstein regime near the disk mid-plane, the Stokes number is defined as $$\text{St=}\frac{\pi \text{a}{\rho }_{\text{s}}}{\text{2}{\Sigma }_{\text{g}}},$$(8)

where a is the size of the grain. If St ≪ 1, then the grains are coupled to the gas and follow its motion in the disk. This has an impact on the dust velocities. TWO-POP-PY divides the different grains into two groups depending on their size (large grain and small grain populations). The velocities of the grains are calculated as an average for each group and are defined as $$\begin{array}{l}{v}_{\text{drift}}=\frac{c_{\text{s}}^2\left|\gamma \right|}{2v_K}\\ v_i=\frac{v_{\text{r},\text{gas}}}{1+{\text{St}}_{\text{i}}^2}+\frac{2}{{\text{St}}_{\text{i}}+{\text{St}}_{\text{i}}^{-1}}{v}_{\text{drift},}\end{array}$$

where i = {0, 1} for each grain population, v_drift is the drift velocity (Birnstiel et al. 2012) and v_r,gas is the radial velocity of the gas. We present the radial gas velocities from the hydro-simulations used here in Appendix A. Studies have showed that filtering of dust to the inner parts depend on the Stokes number of the dust and the gas properties (Weber et al. 2018; Haugbølle et al. 2019), meaning that the TWO-POP-PY approach gives us a first approximation only of how dust is filtered to the inner parts of the disks.

After 1 Myr of evolution, we reconstruct the full grain size distribution, determining the surface density of each grain size as a function of orbital distance (Birnstiel et al. 2015). The particle grid used by the reconstruction routine logarithmically ranges from a₀ to 6 × a_max, a_max being the maximum grain size reached after 1 Myr of evolution, resolved with 300 cells. This dust size-density distribution, presented for each disk in Sect. 3.1, is the final 1D output used to produce the images.

2.4 Synthetic image setup

The images are derived using the radiative transfer code RADMC3D (Dullemond et al. 2012). Using the 1D dust size-density distribution from TWO-POP-PY, we extrapolate the 3D distribution of the grains in the disk. We assume a volume density following $${\rho }_{\text{d}}(r,\varphi ,z,\text{St})=\frac{{\Sigma }_{\text{d}}(r,\text{St})}{\sqrt{2\pi }{H}_{\text{d}}(r,\text{St})}\times \exp \left(-\frac{z^2}{2H_{\text{d}}{(r,\text{St})}^2}\right),$$(11)

where z = rcosθ (θ being the polar angle) and H_d is the dust scale height. H_d is derived from the gas scale height, taking into account vertical settling of grains. The grains having a large Stokes number (Eq. (8)) tend to settle toward the mid-plane, while the grains with small Stokes numbers feel the vertical mixing due to turbulence in the disk and are more coupled to the gas. The dust scale height is therefore given by (Birnstiel et al. 2010; Pinilla et al. 2021) $$H_{\text{d}}=H_{\text{g}}\times \min [1,{(\frac{\alpha }{\min (\text{St},1/2)(1+S{t}^2)})}^{1/2}].$$(12)

The hydrodynamical and dust simulations are run with a very high radial resolution, especially at low viscosity. In order to add two dimensions to our disks, we need to reduce the overall resolution due to computational limitations. To do so, we interpolated the dust distributions radially, reducing the resolution by a factor of two (four at low viscosity). The new radial resolution, used to derive the images, is therefore of N_r = 1109 cells. By applying the same method to the grain size grid, we reduce the number of dust grains by 2. The images are derived considering a disk containing 150 grain size bins. This procedure resulted in no major differences compared to simulations with the full resolution. Thanks to this reduced resolution distributions, the disk can be extended over 320 cells in colatitude (N_θ = 320). Such resolution is needed to correctly derive the temperature of the dust settled to the mid-plane. Finally, we used an azimuthal resolution of N_ϕ = 4 cells, as our disks are considered axisymmetric after the hydro-simulations step.

As mentioned in the previous section, the grains are taken to be spherical grains with the DSHARP composition (Birnstiel et al. 2018). To derive the opacities of such grains, we used OpTool, developed by Dominik et al. (2021). This library derives the opacity for each grain size using the Mie calculation. OpTool computes the full scattering matrices needed by RADMC3D to include the full treatment of anisotropic polarized scattering.

After computing the temperature profiles of each grain size, RADMC3D derives the images at λ = 1.3 mm, corresponding to band 6 of ALMA. We assume face-on disks at a distance of 140 pc, which is the typical distance of observed protoplanetary disks, such as in the Taurus, Lupus, and Ophiuchus regions (Gaia Collaboration 2018). We investigate the influence of the inclination of the disk in Sect. 3.2.3 by showing how an inclination i of 30, 45 or 60 degrees influences our results. In our configuration, RADMC3D uses 10⁷ photon packages and 5 × 10⁶ scattering photon packages to derive the raw images. These images are then convolved with a Gaussian beam of FWHM = 0.04′′ × 0.04′′. In Sect. 3.2.4, we investigate different beam sizes that are coherent with different ALMA configurations. At a distance of 140 pc, the 0.04′′ × 0.04′′ beam size corresponds to a spatial resolution of 5.6 AU × 5.6 AU, which corresponds approximately to the semimajor axis of the inner giant in our different configurations.

Fig. 1 Dust densities distributions for the MMSN-like aspect ratio. Each line represents a Solar System configuration (Compact on the first row and Spread on the second row; cf. Table 1) and each column represents a set of α viscosity and fragmentation speed (increasing from left to right). The vertical dotted lines represent the positions of each planet.The vertical dashed line shows the truncation radius (Eq. (3)). The white lines represent the maximum size reached by the grains, and each line style represents a limiting mechanism: solid line for the fragmentation limit, dashed line for the drift limit, and dotted dashed line for the growth limit size. At high viscosities, dust from the outer regions flows through Jupiter and Saturn’s gap and replenishes the inner disk, which is not the case at low viscosity. At low viscosity, the perturbations induced by the planets in the gas velocity profiles produce dust over-densities (traffic jams), creating substructures not directly related to the positions of the planets.

3 Results

3.1 Dust size distributions from TWO-POP-PY

3.1.1 Solar System configurations

In this section, we present the radial dust size distributions from TWO-POP-PY after 1 Myr of evolution. In Fig. 1, we show the distributions in the different Solar System configurations (rows) at different α viscosities and fragmentation speeds (columns) for the MMSN-like aspect ratio. The distributions with a smaller aspect ratio are presented in Fig. 2. In both figures, the white lines represent the maximal size reached by the grains: the solid line shows the part of the disk where the grains are limited by fragmentation, the dashed one corresponds to the drift limit and the dotted dashed line represents the growth limit. We see that these lines are in general above the dotted horizontal line that marks the 1.3 mm size. Vertical dotted lines represent the positions of each planet, and the vertical dashed line shows the location where the dust disk is initially truncated (see Eq. (3)).

We see that depending on the viscosity, some dust could flow through Jupiter and Saturn’s gaps: as expected, at low viscosity α = 10⁻⁴, the gaps aretoo deep for the dust from the outer disk to replenish the inner region; however, at high viscosity α = 10⁻², the dust diffused through the whole disk, leaving no strong substructures. At an intermediate viscosity (α = 10⁻³), we see that depending on the planet configurations, some dust could accumulate around the giant planets locations: in the Compact configuration with an MMSN-like aspect ratio (Fig. 1, first row, second panel), dust flows from the outer disk and accumulated between Jupiter and Saturn as well as at the inner edge of Jupiter’s gap, creating two over-densities. A similar behavior happened when considering a lower aspect ratio (Fig. 2, first row, second panel). However, as a lower aspect ratio implies deeper gaps, Saturn’s gap becomes deep enough to accumulate dust at the outer edge of its gap. In this case, the inner disk is more depleted because it is harder for the dust to diffuse through.

The depletion of dust at low viscosity creates inner cavities. These cavities are observed in several disks (Espaillat et al. 2014; van der Marel et al. 2018) and are described as wide regions in disks where there is no emission observed and therefore possibly no dust present. These cavities can be explained by the presence of planets: either one giant planet is large enough to block the dust flux from the outer disk and the inner disk empties by radial drift; or multiple planets are present and large enough to create a wide common gap. In our simulations, the cavities are created by either mechanism, or a combination of them, depending on the configuration. Either way, the position of the cavity is linked to the position of Jupiter and Saturn, which are both located within 15 AU. In the next section (Sect. 3.2) we investigate if the resolution of ALMA is sufficient to see the cavities in our cases.

Comparing the Compact and Spread configurations distributions at α = 10⁻³ for each aspect ratio, we see that there is less dust flowing to the inner disk in the Spread configuration than in the Compact one. This can be explained by the quasi-common gap created by Jupiter and Saturn in the Compact configuration: as Saturn is located further in, the dust has to diffuse through one gap while it has to diffuse through two distinct gaps in the Spread configuration. Haugbølle et al. (2019) found a similar behavior: the presence of a common gap containing Jupiter and Saturn makes filtering less efficient.

Focusing on the Compact configuration with a small aspect ratio for a medium viscosity (Fig. 2, first row, second panel), the inner disk is not completely depleted in dust and limited by growth and drift. Reducing the aspect ratio influences the dust behavior in two ways: first, the planet gaps are deeper (Crida et al. 2006) and therefore filter dust more efficiently. Moreover, diffusion is reduced as the viscosity is lower for lower aspect ratios. It is therefore harder for the dust to flow through Jupiter and Saturn’s gaps. As only some part of the small dust manage to reach the inner part of the disk, the dust-to-gas ratio is very low. As a result, the growing dust particles are limited by growth in the inner region before they drift away, in contrast to the simulations where a lot of dust can reach the inner disk, resulting in grain sizes limited by fragmentation. The capacity of the planet gaps to filter out larger dust can therefore produce an inner cavity with very low surface densities in dust, without being directly linked to the positions of the planets in the disk.

In each panel of both figures, the white lines show that all the disks are roughly fragmentation-limited at radii r < 50 AU and are growth-limited at larger radii. This is due to the fact that the timescales are longer at larger radii (differential rotation). This growth limit sets the size of the disks in the millimeter dust to be ~ 50 AU. This size is dependent on the time length of the dust simulation: the millimeter dust disk starts by increasing in size until radial drift reduces the millimeter dust disk significantly. We discuss the impact of time evolution later in Sect. 4.4.

At the connection between the growth limit and fragmentation limit (~ 50 AU), in each case, we observe a depletion of small dust (a ≲ 10⁻⁴ cm). This depletion is due to a narrow region where the dust is actually drift-limited, between the growth and fragmentation limits. We clearly see this regions at high viscosity (dashed line around 30 AU). This drift-limited region makes the largest grains drift inward, before becoming fragmentation-limited. This induces a depletion in small dust as there is no mechanism to replenish this region after the growth of the small dust. This phenomenon was studied by Birnstiel et al. (2015) and creates a gap in the small dust that is not linked to planets at all, but rather to grain growth and drift.

In the Compact configuration at low viscosity and MMSN-like aspect ratio (Fig. 1, first row, first panel), we clearly see over-densities that are not linked directly to the gaps caused by planets. These over-densities, not linked to pressure bumps in the gas disk, are due to the highly perturbed gas velocities. As the gas is highly perturbed by the presence of multiple planets (see Appendix A), the small dust coupled to the gas undergoes “traffic jams”: the change of velocity creates over-densities as the dust is slowed down. In these cases, the dust is not trapped and continues to flow after staying in the traffic jam. However, this dust caught in a lower velocity region has time to grow, creating over-densities over a wide range of different sizes. These over-densities are therefore indirectly linked to the presence of the planets and can be observed depending on the resolution of the instruments (see Sects. 3.2 and 3.2.4).

Fig. 2 Same as Fig. 1 but for the smaller aspect ratio. A smaller aspect ratio induces deeper planet-induced gaps in the gas disks, creating stronger features in the dust compared to Fig. 1. In the Compact configuration, at α = 10−3, a small amount of small dust flows through the gaps of the giants. This produces a low dust-to-gas ratio in the inner region of the disk, resulting in grain sizes limited by growth and drift rather than fragmentation (see Eq. (7)).

Fig. 3 Dust densities distributions with the planet masses reduced by one-half and two-thirds in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The masses of the planets mainly change the gap shapes, allowing more or less dust to flow to the inner regions. In the Three-Giants configuration, we also see that more massive planets create more traffic jams in the disk and therefore present more substructures.

3.1.2 Influence of planet mass and Three-Giants configuration

As presented in Sect. 2.1, we also investigate the impact of different planetary masses on the dust distributions in the Spread and in the Three-Giants configurations. In the first row of Fig. 3, we show the dust distributions in the Spread configuration case, at low viscosity and small aspect ratio, where the masses are reduced by a factor of two-thirds (left panel) and one-half (middle panel). We can compare them to the total mass case, presented on the right panel. On the second row, we present the distributions of the Three-Giants case, with the different masses as mentioned above.

Changing the mass of the planets will have two large impacts. The first comes from the initial gas distribution: with increasing planetary mass, the gas is pushed away from the planet more efficiently, depleting the gas disk (Bergez-Casalou et al. 2020). Our simulations start with lower disk masses for more massive planets to depict the effect of disk evolution during planetary growth. Consequently, as the initial dust content is derived from the gas profile (dust-to-gas ratio of 0.01), the dust disk in the case of full planetary masses is less massive than in the other two cases. As consequence, the maximal grain size reduces for the simulations with increasing planetary mass. Having a more massive dust disk also has an impact on the size of the dust disk in millimeter grains, as it allows faster growth in the outer regions of the disks.

Furthermore, the less massive planets cause shallower gaps. As the gaps are less deep, more dust can flow through the different gaps. For the Spread configuration, the planetary masses results in narrower gap in the dust profile at Saturn’s location, while Jupiter is still massive enough to prevent efficient dust diffusion. Regarding the Three-Giants case, reducing the planet masses does not alter the formation of an inner cavity. Even if some dust diffuses through the gap of the third planet in the case of smallest planetary mass, the presence of the other giants is sufficient to keep an inner cavity. Therefore, as expected, bigger planets will create bigger cavities.

On the other hand, planets of different masses will perturb the gas disk differently. More massive planets will induce more perturbations in the gas and create more traffic jams (see the velocity profiles in Appendix A). This is particularly clear in the Three-Giants case, where we see on the lower panels of Fig. 3 that the outer disk presents different over-densities depending on the planet masses.

These dust distributions show that each configuration harbors different substructures, mostly at low viscosity. In the next sections we show that some of these features are observable with ALMA.

3.2 Synthetic images: RADMC3D outputs convolved with Gaussian beams

The dust distributions studied in the previous section present different features, unique to each configuration. In this subsection, we present the images and their the radial profiles derived following the setup presented in Sect. 2.4. First, we focus on the radial profiles at λ = 1.3 mm for the Solar System configurations at each aspect ratio (Sect. 3.2.1). The corresponding images can be found in Appendix C.1. Then we present the images and profiles in the Three-Giants and Spread configurations with the different planetary masses (Sect. 3.2.2). Different disk inclinations are investigated in Sect. 3.2.3 before studying the influence of the beam size on the observable features in Sect. 3.2.4.

Fig. 4 Radial profile of convolved and unconvolved images at λ = 1.3 mm with the MMSN-like aspect ratio. Each row represents a Solar System configuration, and each column represents an α viscosity. The solid lines show the radial intensity of the images normalized to the peak intensity after convolution with a beam of FWHM = 0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU. The beam is represented with a black horizontal line in the upper-right corner of each panel. The dashed lines represent the normalized intensity of the unconvolved image. Vertical lines show the positions of the planets in each configuration. The light blue area shows the region where the normalized intensity is smaller than Fmin ∕Fpeak, where Fmin = 10 μJy is the minimal flux considered to be observable. The value of Fmin∕Fpeak is different from each panel as Fpeak is unique to each image, whereas Fmin is fixed. Comparing the profiles in the case of the convolved and unconvolved images shows us how many features can be missed due to a too low resolution. Regarding the Compact configuration, all substructures are smeared out in the beam.

Fig. 5 Sameas Fig. 4 but for a smaller aspect ratio. As the dust distributions show more intense substructures, some features become observable but the majority are still smeared out in the beam.

3.2.1 Radial profiles in the Solar System configurations at λ = 1.3 mm

In order to determine which features are observable, we show in Figs. 4 and 5the radial profiles of the normalized intensities (intensity of the image normalized to the peak intensity along one radius of the disk) for images with unconvolved (dashed) beams and for images with a 0.04′′ × 0.04′′ Gaussian beam convolution (solid). The corresponding images can be found in Appendix C.1. In the images, we assumed that the minimum flux that can be received due to noise is F_min = 10μJy beam⁻¹ (Andrews et al. 2018). This minimum flux is represented in the radial profiles by the blue regions: in each profile, this minimum flux is normalized to the peak intensity and emission present in this region can be assumed to be lost in the noise of the images. The value of F_min∕F_peak is different for each panel as F_peak is unique to each image, whereas F_min is fixed.

As in the previous figures, each row represents a planet configuration and each column represents an α viscosity. It should be noted that for a better comparison with the images, we present the profiles on a linear radial scale, whereas in the previous section the dust distributions are presented along a logarithmic radial scale.

Regarding the high viscosity cases (α = 10⁻²), as expected from the dust distributions, the disks show almost no noticeable feature. In the unconvolved images, we can distinguish the gap created by Jupiter. However, the gap is too close to the star and too small to avoid being smeared out by the beam. Therefore, if the viscosity is too high, then a Solar-System-like planetary structure would be completely invisible in the dust disk. This is consistent with the work of de Juan Ovelar et al. (2016), where they show that a high viscosity disk does not present strong substructures. It is also consistent with Zhang et al. (2018) where they show in their Sect. 5.1 that a Solar System protoplanetary disk featuring our giant planets in nowadays configuration would not present strong substructures if the viscosity is too high.

The images at α = 10⁻³ show a similar pattern: the features are mostly either too small or too close to the star to be distinguishable with this resolution. Although, in the Spread configuration case, the inner disk starts to be depleted in dust (see Figs. 1 and 2), resulting in a decrease in the normalized intensity in the inner regions of the images. These small inner cavities are located at a radius within Jupiter’s orbit: the giant planets are therefore outside the cavity in this case. In this same configuration, Saturn’s gap start to be large enough to induce a small dip in the intensity profile, slightly noticeable in the convolved images.

At low viscosity (α = 10⁻⁴), the dust distributions showed strong inner cavities and several substructures. The convolution with a beam of this size kept the inner cavities in all configurations, even if they are deeper in the Spread configuration than in the Compact one. In all the configurations, the intensity profiles decrease rapidly around 50 AU, which is caused by the growth limit of the grains discussed in the previous section. As the millimeter size dust is the dust that contributes the most to the 1.3 mm emission, the growth limit sets the size of the disk in the images (as can be seen in Appendix C.1). In the Spread configuration, depending on the sensitivity of the instrument, Neptune can be located close to the edge of the disk but only the growth limit sets the location of the drop of intensity.

As discussed in the previous section, the multiple planets can create some substructures in the disk not directly linked to the positions of the planets. Even if these over-densities can be seen in the unconvolved images, the beam smeared outthe majority of them. However, in the Compact case, low viscosity, low aspect ratio (first top panel of Fig. 5), we see several dips in the intensity profile. The first one is linked to the inner cavity (r < 10 AU). The second one is located at Uranus and Neptune orbits and originates from the small gaps that the two icy giants create in the gas and dust disk. However, the two gaps are indistinguishable here due to the beam size, reducing the emission of the dust located between the two planets. Similarly, as some dust piled up at the outer edge of Neptune’s gap and due to the shape of the fragmentation limit in this case, a small part of the disk is shadowed outside of Neptune’s orbit, creating a third dip in the intensity profile.

Focusing on the Spread configuration at low viscosity (α = 10⁻⁴) at each aspect ratio, Figs. 4 and 5 show a bump in the intensity between the orbits of Saturn and Uranus. This bump originates from the pileup of dust that is blocked at Saturn’s outer gap edge. This bump creates a bright ring separating the inner giants and the icy giants, located around 15 AU. This configuration and viscosity is the only setup that presents a bright clear ring at this resolution. We discuss this peculiarity in Sect. 4.1, where we compare our disks to known observed disks with similar resolutions.

In summary, the Solar System configurations do not present a lot of substructures in general at this resolution. At low viscosity, the Compact and Spread configurations are presenting very different features, with the Spread configuration showing a clear bright ring between Saturn and Uranus, whereas the Compact configuration presents features that are not directly linked to the positions of the planets. Therefore, the detectability of substructures is highly dependent on the disk viscosity and planet configuration.

Fig. 6 Radial profile of convolved and unconvolved images at λ = 1.3 mm with the planet masses reduced by one-half and two-thirds in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The intensity profiles are presented as in Fig. 4. The convolved profiles are convolved with a beam of 0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU, represented with a black horizontal line in the upper-right corner of each panel. One can notice the similarities between each convolved profile, making it difficult to disentangle between each evolutionary stage and configuration. In the Three-Giants configuration, we highlight in orange the substructures outside of the giants’ region, originating in the traffic jams discussed in Sect. 3.1.

3.2.2 Influence of the different masses on the 1.3 mm images

As discussed previously, the planet masses have a lot of impact on the dust distributions, creating different features in the disks. We show in Fig. 6 the radial normalized intensity profiles in the simulations with different planetary masses, for the Spread and the Three-Giants configurations. Due to the different masses of the dust disk, they have different sizes: in the Spread configuration (top row), we see that the drop in intensity due to dust growth is located at different radii. This effect is even more noticeable in Fig. 7, where we show the disks as they would be observed. This effect is less present in the Three-Giants case because the planets are located in the very inner region, therefore not having a strong impact on the gas distribution in the outer regions of the disk (see Appendix A).

As noticed in the previous section, having planets of different mass influence the amount of dust that can flow through the gaps of the giants. This has a very small impact here as the beam is too large to resolve the small amount of dust present between Saturn’s gap (or the third giant’s gap in the Three-Giants case) and the inner cavity. The only case where enough dust managed to flow to the inner regions is in the reduced by two-thirds case in the Spread configuration, as Saturn is only slightly above pebble isolation mass. Here, the amount of dust is large enough to slightly enhance the intensity between Jupiter and Saturn’s orbit. This slight enhancement is particularly visible in the image in Fig. 7 (first top panel). On the other hand, as more massive planets block dust most efficiently, it accumulates more at the outer edge of Saturn’s orbit. This results in a brighter and clearer ring separating the inner giants from the icy giants.

Regarding the features created by the different traffic jams, most of them are not strong enough to be noticed in the Spread configuration. However, in the Three-Giants case, the traffic jams create rings of different intensities. Expectedly, strongest over-densities create brighter rings. As the strength of the traffic jams are dependent on the planet masses, the most massive case present the strongest features. In the end, these two consecutive rings are due to the gas radial velocity profile and are not directly linked to the presence of planets close to the bright rings.

As already mentioned, changing the planet mass allows us to probe different stages of planet formation (i.e., different times in the formation process). One can notice that the differences between the Spread configuration with half its mass is quite similar to the reduced planet masses in the Three-Giants configuration (Figs. 7 and 6, middle panel of first row compared to the first and second panels of the second row). Therefore, in order to really disentangle planet formation processes, better resolution is needed. We discuss this in Sect. 3.2.4.

Fig. 7 Images in total intensity at λ = 1.3 mm corresponding to the radial profiles presented in Fig. 6. The positions of the different planets are represented by the different white arcs. The beam size (0.04′′ × 0.04′′) is represented in the lower-left corner of each image by the white ellipse. In the Spread configuration, we see that the size of the disks depends on the masses of the planets. In the Three-Giants configuration, one can notice the substructures outside of the giants’ region corresponding to the orange regions of the profiles in Fig. 6 and originating from the traffic jams observed in Sect. 3.1; these traffic jams were created by the perturbed velocity profile of the disk (see Fig. A.7).

3.2.3 Influence of the disk inclinations

Many of theobserved disks are actually inclined compared to our line of sight (ALMA Partnership 2015; Andrews et al. 2018). In this section, we explore how the inclination of these disks can have an impact on the visible substructures discussed in the previous sections. In order to do so, as in the previous section where we explored the impact of the beam size, we derive the images in the Spread configuration at low aspect ratio and viscosity and in the Three-Giants configuration and infer three different inclinations to the disks: i = 30, 45, 60°. The images at different inclinations can be found in Appendix C.2.

In Fig. 8, we present the radial profiles of the normalized intensity with different inclinations. The radial profiles are taken to be a section of the image along the semimajor axis of the inclined image. No deprojection procedure was applied, as the disks are axisymmetric by construction. By taking the profile along the semimajor axis, we look at the section of the disk that is situated at the same distance from the observer, independently of the inclination. As in the previous profiles, we show the unconvolved profiles with dashed lines and convolved profiles with solid lines.

The profiles are very similar, presenting the same features in each case. The main difference resides in the inner cavities: a more inclined disk will hide the depth of the inner cavity as the dust present closer to the observer (lower part of our images) will hide the cavity. However, the cavities are wide enough to be visible: if the cavity is too small,the dust from the closer part of the disk will completely hide the cavity to the observer. On the other hand, larger cavities will be less impacted by the disk inclination. This can be seen with our configurations: in the Spread configuration, the inner cavity is smaller than in the Three-Giants one (see Figs. 3 and 7) and the inclination of the disk has a stronger effect on the inner cavity in the Spread configuration.

As disk millimeter emissions are geometrically flat (Birnstiel et al. 2010; Pinilla et al. 2021), we show here that the inclination of the disk does not have a strong effect on the observed profiles. This means that the inclination of the disk does not hide or create features that could originate from giant planets, important to derive constraints for planet formation.

Fig. 8 Radial profile of convolved and unconvolved images at λ = 1.3 mm with different inclinations in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). Here, the disks have the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The convolved images present a beam of 0.04′′ × 0.04′′. The profiles are derived from images that have inclinations ranging from i = 0° to i = 60°. The light blue area shows the region where the normalized intensity is smaller than Fmin ∕Fpeak (see Fig. 4). The only highly impacted region is the inner cavity: more inclined disks hide the inner regions more efficiently, influencing the depth of the observed cavities.

3.2.4 Influence of the beam size

ALMA can reach different resolutions depending on the observed wavelength and configuration. As we derived the images at λ = 1.3 mm, we are interested in the Band 6 observations. With the different configurations available, the most common resolutions reached are therefore equivalent to beams of three different size: the most resolved one has a beam of 0.02′′ × 0.02′′ (as in Benisty et al. 2021), the most common one has a high resolution with a beam of 0.04′′ × 0.04′′ (as in Andrews et al. 2018) (used in the previous sections) and the last configuration gives a beam of 0.1′′ × 0.1′′ (as in Long et al. 2018; Kurtovic et al. 2021). In Fig. 9, we present the radial profiles of the normalized intensity in the Spread and Three-Giants configurations, with a small aspect ratio and low viscosity. The corresponding images are shown in Fig. 10. As we present the images in intensity per beam, the sensitivity are different for each resolution: the colorbars range from F_min to different maxima depending on the resolution.

In the Spread configuration, we noticed in the previous section that the dust situated between theinner giants and the icy giants creates a bright ring when the resolution is 0.04′′ × 0.04′′. When the resolution is lower (right panels), the emission is spread over Jupiter and Saturn’s orbit, completely hiding Saturn’s gap. However, the inner truncated disk is still noticeable, with a small decrease in the intensity. On the other hand, when we compare the high resolution case (at 0.04′′) with the highest resolution (at 0.02′′), we see that the ring is clearly located at the outer edge of Saturn’s gap. The resolution is even high enough to start distinguishing Uranus’s gap and the small over-density of dust located between Jupiter and Saturn. On the other hand, with a minimum flux situated at 10 μJy beam⁻¹, Neptune is completely missed and lost in the noise.

Regarding the Three-Giants configuration, we see a similar behavior: with a beam of 0.1′′ × 0.1′′ (right panels), the substructures created by the perturbations by the giants are completely washed out and the only feature remaining is the inner cavity. In both configurations, one can notice that the cavity is shifted compared to the planets orbits: the giants are located in the decrease in intensity, not at the minimum. As we noticed in the previous section, a resolution of 0.04′′ × 0.04′′ is sufficient to start to distinguish the over-densities of dust located at the outer edge of the further giant gap. The highest resolution is needed to really resolve the main over-densities that are due to the perturbations of the gas velocity by the multiple giants. With this very high resolution, two rings are observed, corresponding to the two brightest set of over-densities seen in the dust distributions (bottom left panel of Fig. 3).

Changing the resolution has an impact on the observable features, but has also an impact on the observed size of the disk. As the beam spread the intensity in the disk, this enhances the value of F_peak, reducing the value of F_min∕F_peak, as we assume a fixed F_min value. Therefore, as the intensity decreases with the radius, the observed size of the disk will depend on the resolution used, presenting a larger disk at lower resolutions. As the disk size can change depending on the image resolution, alternative methods such as uv-modeling (Hendler et al. 2020) should be considered to retrieve such a quantity.

4 Discussion

4.1 Comparisons to known observed disks

The DSHARP survey (Andrews et al. 2018) studied several bright massive disks around stars located in the vicinity of the Sun with a beam size of~0.035′′. These disks present several features, such as gaps, rings, spirals and asymmetries. The most axisymmetric disks show several configurations of gaps and rings (Huang et al. 2018). Some of the disks present bright rings located at large radii (r > 50 AU), such as AS209 and HD163296. Our synthetic images never showed features at such large radii, even in the Spread configuration where Neptune is located at 39 AU. This can be explained by the sizes of the disks: our gas disk is small (160 AU in radius) whereas HD163296 is thought to be wider than 500 AU in radius (Isella et al. 2007; Muro-Arena et al. 2018) with a potential outer planetlocated at 260 AU (Pinte et al. 2018).

One of the explanations for the existence of structures at large radii is the presence of planets carving gaps and creating rings. Lodato et al. (2019) show that with giant planets migrating fast enough, it is possible to produce gaps and rings at large radii and still reproduce the distribution of eccentric giant planets observed in radial velocity. However, such migration speeds require a too high viscosity (Ndugu et al. 2019) compared to the viscosity needed to allow the formation of planets that could explain the observed substructures. On the other hand, it is also very unclear how planets can form that far in disks with the core accretion model (Morbidelli 2020). Moreover, our images with the Three-Giants configuration do not present features in the outer disk because the planets are located in the inner disk. It is possible that this configuration will lead afterward to some scattering events that will produce systems with giant eccentric planets (Bitsch et al. 2020). We show here that planets forming in the inner disk do not result in features (rings or gaps) in the outer disk region as observed in the DSHARP survey. It is clear these two different giant planet formation channels result in different observable disk structures.

Another important point is linked to the substructures induced by the planets but not directly linked to their orbit and gap. For example, Zhang et al. (2018) analyze the rings and gaps structures present in the DSHARP disks and derive which possible planet mass could produce such substructures. Even if they take into account the fact that some planets can create multiple gaps at low viscosity (Dong et al. 2018; Bae & Zhu 2018), as in AS 209, we found that the gas radial velocity structure can also create rings, blurring even more the link between the number of planets and the number of gaps present in the disk. We discuss in Sect. 4.3 how this problem could be addressed.

However, the disk surveys are biased toward the brightest disks. The differences in the images between the massive DSHARP disks and our study confirm that planet formation happening in the inner regions of the disk results in different features in observed disks. However, these surveys contain some bright disks that are similar in size and show comparable features as the disks studied here. In the Ophiuchus DIsc Survey Employing ALMA (ODISEA) (Cieza et al. 2019; Williams et al. 2019), DoAr44 presents a bright inner ring and a dimmer one exterior to it, resembling the image of the Spread configuration disk at low viscosity and low aspect ratio with a resolution of 0.02′′ (see Fig. 10, first top panel). Similarly, Facchini et al. (2020) observed two disks, LkCa15 and J1610 showing features similar to our Three-Giants configuration disk with a resolution of 0.04′′ (see Fig. 10, middle panel of second row). Observations of the V4046 Sgr circumbinary disk by Martinez-Brunner et al. (2022) present features that are very similar to the Spread Solar System observed with a similar resolution (in Fig. 10, top left panel), around a binary, unlike our configuration. This images prove that ALMA is capable of reaching such high resolution. Observations of such disks can therefore give us some insights on how planet formation can occur in the inner regions of the disks, compared to the DSHARP observations giving insights on how it occurs in the outer regions of the disks. Moreover, the constrains derived from the local study of the Solar System could help understand planet formation in disks such as DoAr44 that present similar features as our Solar System disks.

Fig. 9 Radial profile of convolved and unconvolved images at λ = 1.3 mm with different beam sizes in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The different beams investigated range from 0.02′′ × 0.02′′ = 2.8 AU × 2.8 AU (left) to 0.1′′ × 0.1′′ = 14 AU × 14 AU (right). We present the fiducial resolution (0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU) in the middle panel for comparison. Each beam is represented by a horizontal line in the upper-right corner of each panel. Orange regions represent the substructures produced by traffic jams, clearly visible in the Three-Giants configuration. The different profiles show that the highest resolution is really needed to start to correctly represent the features of the dust disk.

Fig. 10 Images in total intensity at λ = 1.3 mm corresponding to the radial profiles presented in Fig. 6. The positions of the different planets are represented by the different white arcs. The different beams investigated range from 0.02′′ × 0.02′′ (left) to 0.1′′ × 0.1′′ (right) and are represented in the lower-left corner of each image by the white ellipses. We show the fiducial resolution (0.04′′ × 0.04′′) in the middle panel for comparison. The colorbars are adjusted for each resolution as they have different sensitivities. These images show the importance of resolution: the substructures start to be well represented at very high resolution.

4.2 Comparing to exoplanet populations

In Sect. 3.2.4, we show that with the highest resolution, it is possible to observe features originating from the ice giants if they are in the very outer regions of the disk and far away from the inner gas giants. Microlensing surveys, such as that presented in Suzuki et al. (2016), claim that the most common type of planets observed are ice giants located at a few AU from their star. Our study shows that it is possible to observe, with the highest resolution that ALMA can reach, features caused by such planets if they are at a few tens of AU. These observations could therefore help to constrain the formation pathways of the ice giants found in microlensing surveys, under the assumption that these ice giants do not turn into gas giants. We should note here that microlensing surveys mostly observe dwarf stars, which should have less massive disks in the first place, making observations unfortunately very difficult.

Constraining planet formation during the disk phase is important to improve our understanding of different formation scenarios. Indeed, the gas disk phase contains information about the initial conditions of planet formation and the initial structure that could lead to dynamical instabilities after the disk phase. The final structure of the planetary system depends highly on the processes happening during the gas disk phase.

Currently, the formation of giant planets is still unclear. In one hand, observations of large disks, such as in the DSHARP survey, motivate the idea that giant planets must form in the outer part of the disks and then migrate inward, explaining the presence of bright outer rings and the planet distributions observed by different surveys (Lodato et al. 2019). However, this scenario requires a rather high viscosity in order to have an efficient migration of the giant planets. Ndugu et al. (2019) shows that if the viscosity of the disk is lower, as disk observations seem to suggest (Dullemond et al. 2018; Flaherty et al. 2018), then these giant planets do not have time to migrate to semimajor axis corresponding to distances within the reach of radial velocity surveys (e.g., Fulton et al. 2021).

Another possible giant planet formation scenario is to have giant planets forming in the inner regions of the disks, where the orbital timescales favor planet formation and where a slower migration of the planets can still explain the observed giant distributions (Bitsch et al. 2020). Our study shows that giant planets forming in the inner part of the disks do not produce bright features as observed in the DSHARP survey.

Higher resolution observations of disks can therefore help us distinguish between the formation of planets in the outer disk or in the inner regions of protoplanetary disks. This can give constraints on the initial conditions needed for planet formation to occur and improve the link between the different observed planet populations and the theoretical models studying different planet formation scenarios.

4.3 Features created by traffic jams

In Sect. 3.1, we present the dust distributions in each configuration. Some of the distributions show multiple narrow dust over-densities, especially at low viscosity. The configuration showing the clearest over-densities is the Three-Giants configuration (Fig. 3). These dust rings are created by traffic jams, as shown in Appendix B. As these traffic jams originate from the highly perturbed gas radial velocities and not from a pressure bump present in the gas, the dust is not trapped and will flow to the inner parts of the disk.

The presence of these traffic jams has several impacts on our understanding of planet formation. First, as they create features observable by ALMA, it blurs further the link between the number of planets present in disks and the number of gaps and rings created. Considering that a single planet can create multiple gaps in low viscosity disks (Dong et al. 2018; Bae & Zhu 2018), having features created by velocity perturbations on top of the one created by pressure perturbations complicates our estimations of planet masses needed to create observed features.

However, in order to trigger the formation of these traffic jams, the disk needs to be highly perturbed in velocity. In our simulations, it requires the presence of multiple giant planets. In systems where only one or two giant planets are embedded, the velocity perturbations do not create strong traffic jams (Pinilla et al. 2015). The presence of traffic jams is therefore linked to the presence of the ice giants in our simulations, meaning that their impact is non-negligible on the dust substructures. Depending on the masses of these planets, many of the disks observed might therefore present features originating in some traffic jam effect rather than from dust trapped in pressure bumps. This effect has been encountered in the past (e.g., Rosotti et al. 2016) and some studies show that it is possible to disentangle between an over-density of dust caused by a pressure trap or by a traffic jam (Pinilla et al. 2017a,b; Dullemond et al. 2018). Observations at multiple wavelengths is a possible way to distinguish between each mechanism and can therefore help to unveil the number of planets contained in the observed disks.

Another method that can be used to determine if an observed ring is coming from a pressure trap or a traffic jam is to study the CO velocity perturbations (Teague et al. 2018; Pinte et al. 2018, 2020). In these studies, the presence of a pressure gradient in the disk can be linked to a change of rotational velocity. Traffic jams, originating from a perturbation of the radial gas velocity, would not influence the rotational velocity profile of the disk. The presence of this “kink” in the CO rotational velocity perturbations could be used to distinguish between a traffic jam or a pressure bump (Izquierdo et al. 2021).

Fig. 11 Images at λ = 1.3 mm of the different configurations, at low viscosity and low aspect ratio, at different times: 0.5 Myr, 1 Myr, 2 Myr, and 3 Myr, from left to right. The white lines represent the positions of the different planets in each configuration. We assume that the sensitivity of the instrument is limited to fluxes larger than 10 μJy beam−1. As time evolves, the sizes of the disks shrink due to inward drift, leaving stable rings in the inner disks after 1 Myr.

4.4 Impact of time evolution

With this project, we made the choice to implement dust growth with TWO-POP-PY at the expense of a 2D or 3D evolution model, allowing us to study the time evolution of dust growth. As Drążkowska et al. (2019) show in their 2D study, including dust coagulation has a non-negligible impact on the dust distribution. In our case, we implemented dust growth by assuming that the gas is fixed during dust growth (here evolved for 1 Myr). This assumption makes a good approximation at low viscosity for two reasons: the first one is that the viscous timescale is way larger than 1 Myr, meaning that the disk would remain almost static during this time; on the other hand, migration is slow at this viscosity (Baruteau et al. 2014), meaning that the dynamics of the planets would not strongly influence the gas disk structure. Even if instabilities can be triggered at low viscosities creating azimuthal asymmetries (as in Zhang et al. 2018, Sect. 5.1 with a low viscosity Solar System disk), we assume that these asymmetries vanished by the time the planets are fully formed, as supported also by other hydrodynamical simulations (e.g., Hammer et al. 2017; Bergez-Casalou et al. 2020) and shown in Appendix A.

However, at high viscosity, gas evolution and planet–disk interactions over 1 Myr start to be non-negligible: the gas is accreted toward the star and planets migrate faster. Migration of planets can alter the dust distributions in the disk (e.g., Meru et al. 2019; Weber et al. 2019), but, as mentioned in Sect. 3.1, the substructures in the gas (e.g., pressure bumps) are not strong enough to block the inward diffusing dust, preventing the creation of notable features in the disk. Therefore, our setup is a good approximation to estimate how dust is distributed and further studies would be needed to detail the impact of gas evolution.

As mentioned, in this paper, we chose to evolve the dust for 1 Myr during which the gas profile is considered constant. However, as Eqs. (5) and (7) show, the millimeter dust disk size will expand with time as dust will grow further and further out in the disk, until drift will deplete the outer regions of the disk and reduce the millimeter disk size. In Figs. 11 and 12, we present the images and their radial intensity profiles at different times: after 0.5 Myr, 1 Myr (as in Sect. 3.1), 2 Myr, and 3 Myr. After 0.5 Myr, the dust disk is large as drift only starts to deplete the outer region of the disk, showing more substructures than at later times. In the Spread configuration, the gaps created by Uranus and Neptune are slightly distinguishable before being washed away by drift. In the Three-Giants case, the velocity perturbations induced by the planets create numerous gaps and rings in the outer regions, before being washed at later times by drift as for the other configurations. The presence of numerous substructures in young disks matches the observations of the disk around IRS 63, supposedly younger than 0.5 Myr (Segura-Cox et al. 2020).

After 2 Myr, the majority of the dust located in the outer disk had time to grow and drift to the inner regions. The only dust remaining is the dust trapped in the pressure bump located at Saturn’s gap or at the outer giant gap. Long et al. (2020) also investigated the impact of time evolution of the dust size of the disk: they show that without any dust traps, the millimeter dust size of the disk increases until drift reduces the disk; on the other, in presence of a pressure bump created by a planet and acting like a dust trap, the size of the disk is first dominated by growth and then by the position of the dust trap.

Furthermore, at low viscosity, some asymmetries in the gas can arise due to instabilities (e.g., Rossby wave instability Lovelace et al. 1999; Li et al. 2001): to study these asymmetries, a 2D (at least) analysis has to be done. Future simulations, with sufficient computational power, should then consider both dust growth in multidimensional grids (Drążkowska et al. 2019). In our study, we avoided these asymmetries by simulating the gas disk for sufficient time, letting the possible instabilities dissipate in the disk (Hammer et al. 2017).

Time evolution also has an impact on the characteristics of the planets. Indeed, we assume that the planets already have their final mass and do not migrate. Bergez-Casalou et al. (2020) showed that gas accretion can have an impact on the gap depth, specially at low viscosities. This would influence the pressure profile of the disks and therefore the dust distributions and the created substructures. Computing the evolution of gas, planet mass and semimajor axis with dust growth and evolution simultaneously would be ideal but too computationally expensive. In order to investigate what the effect of the planet mass can have on the results, we showed in Sect. 3.2.2 how the images could be different at different stages of planet evolution. The similarities between the images showed us that it is more important to improve the resolution of our observations, as a resolution of 0.02′′ or better is needed to disentangle between different planetary systems.

Fig. 12 Radial intensity profiles after 0.5 Myr, 1 Myr, 2 Myr, and 3 Myr of evolution in each configuration with low viscosity and a low aspect ratio. The total size of the disk depends on the time evolution, which determines the dust distribution.

4.5 Disk dust mass

One of the challenges of planet formation is to overcome the too large radial drift speed of dust. Cosmochemical studies show that planetesimal formation happened at all times during the Solar System disk age (e.g., Connelly et al. 2012). Therefore, if the dust drifts too rapidly during the disk lifetime, then there is not enough material to to ensure a continuous planetesimal formation as in the Solar System. One way to prevent fast inward drift that empties the disk is to trap the dust (Pinilla et al. 2012), as showed in Sect. 3.1. We investigate here how much dust is actually remaining in our disks as a function of disk parameters and planet configurations. In Fig. 13, we show thetotal dust mass contained in each of our disks at two different times: initially and after 1 Myr of evolution.

As expected from Sect. 3.1, the high viscosity disks (α = 10⁻²) do not present strong dust traps and the high diffusion allows the dust to drift toward the star. Therefore, after 1 Myr,they lost a significant amount of dust. On the other hand, all the disks with lower viscosities (α = 10⁻³ and 10⁻⁴) trap the dust efficiently. Even if the dust is distributed differently as a function of time (see Figs. 11 and 12 from the previous section), the amount of trapped dust available to form the small bodies of our Solar System stays constant (from the time the giant planets form and assuming that the inner disk is already depleted in dust).

In all the disks with efficient dust trapping, the total amount of dust is above 100 M_⊕. This can be compared to the solid masses needed to form the small bodies of the Solar System. The solid budget needed to trigger the Nice instability is of at least 20 M_⊕ (Nesvorný & Morbidelli 2012; Nesvorný et al. 2013). As the mass contained in the asteroid belt is negligible (5 × 10⁻⁴ M_⊕, Kresak 1977), our dust disks have enough material to form the small bodies of the Solar System even if the formation of planetesimal from pebbles is not 100% efficient. Moreover, having the majority of the dust located in the outer disk (i.e., outside Saturn’s location) is in agreement with Izidoro et al. (2021) where the authors show that the inner Solar System, after formation of some planetesimals, should be depleted in pebbles to explain the formation of the current terrestrial planets. This also clearly indicates that the viscosity of the gas in the protoplanetary disk must have been low, because at high viscosity even Jupiter is not able to block inward drifting pebbles, which is required by cosmo-chemical studies (e.g., Kruijer et al. 2017; Weber et al. 2018).

On the other hand, dust masses derived by observations are on the order of a few Earth masses to a few tens of Earth masses (Andrews et al. 2013; Ansdell et al. 2016). In order to understand the origin of this discrepancy between the dust masses in our simulations and the one derived by observations, we derive the dust masses from our disk images using the same methods as in the observations and compare them to the actual dust mass of our simulations. Assuming an optically thin disk, the observed dust mass M_d,obs is given by (Hildebrand 1983) $$M_{\text{d},\text{obs}}=\frac{F_{\nu }{d}^2}{{\kappa }_{\nu }^{\text{abs}}{B}_{\nu }(T_{\text{d}})},$$(13)

where F_ν is the total integrated flux density of the disk, d the disk’s distance, ${\kappa }_{\nu }^{\text{abs}}$ the absorption opacity at the observed wavelength (here λ = 1.3 mm) and B_ν(T_d) the Planck function at temperature T_d.

In general, observers use the assumptions made by Andrews et al. (2013) regarding the opacity and dust temperature. The opacity is assumed to be ${\kappa }_{\nu }^{\text{abs}}=2.3$ cm² g⁻¹, following Beckwith et al. (1990). The averaged temperature is taken as T_d = 20 K. This value was derived from the observations of Taurus disks by Andrews & Williams (2005), using Eq. (13): the authors assumed a simple disk model in order to determine the disk dust mass and derived the disk average temperature assuming the same opacity as mentioned above. To be consistent with our simulations, we derive the observed masses using the absorption opacity and average temperature from the RADMC3D outputs and compare the resulting masses to the ones obtained by following Andrews et al. (2013). Assuming that mostly millimeter grains contribute to the opacity at this wavelength (Dullemond et al. 2018), we use the absorption opacity derived from OpTool following the assumptions made in Sect. 2.3, resulting in ${\kappa }_{\nu }^{\text{abs}}=2.04{\text{cm}}^2{\text{g}}^{-1}$ for a grain size of a = 0.1 cm at λ = 1.3 mm. The dust disk temperatures and masses can be found in Table 2.

Equation (13) relies on the assumption that the disk is optically thin. However, from the dust evolution models in Sect. 3.1, we know that the dust is trapped in relatively dense rings that can become optically thick. As some significant amount of mass can be hidden in optically thick regions, we compare the observed dust mass M_d,obs to the dust mass contained in the optically thin regions. In order to do so, assuming ${\tau }_{\nu }(r)={\Sigma }_{\text{d}}(r)\times {\kappa }_{\nu }^{\text{abs}}$, the comparison is made with the amount of dust contained in the regions where τ < 1.

By comparing the total amount of dust to the mass contained in the optically thin regions (Cols. 3 and 5 or circles and diamonds in Fig. 13), we show that the majority of the mass is hidden in the optically thick regions. This results in masses $M_{\text{d},{\tau }_{\nu }<1}$ that can be morethan ten times lower than the actual total mass. Therefore, the presence of optically thick rings in disk images can hide a large amount of dust.

We show in Fig. 14 the difference between the masses derived from Eq. (13), $M_{\text{d},\text{obs}}^{\lambda =1.3\text{mm}}$ and the masses contained in the optically thin regions of our disks, $M_{\text{d},{\tau }_{\nu }<1}^{\text{tot}}$. The diamonds and hexagons represent the different disk parameters where the colors match the planet configuration (Table 1). Markers surrounded by a black outline are the masses assuming the Andrews et al. (2013) setup. We clearly see that the masses are under-estimated when using the actual dust temperature average, whereas a colder temperature (as in Andrews et al. 2013) tends to over-estimate it.

The underestimation of the observed dust mass derived with Eq. (13) with self-consistent opacities and disk temperatures originates from the calculation of $M_{\text{d},{\tau }_{\nu }<1}^{\text{tot}}$. As this mass comes from our simulations, it contains all grain sizes present in the disk. However, as $M_{\text{d},\text{obs}}^{\lambda =1.3\text{mm}}$ is based on F_ν (Eq. (13)), it represents only the grains emitting consequently at the observed wavelength. In our case, the image at λ = 1.3 mm is dominated by millimeter grains, meaning that the mass contained in the smaller grains is not probed here, leading to aless massive dust disk. As this mass underestimation is therefore expected, this shows that the T_d = 20 K assumption is leading toan unrealistic overestimation of the optically thin dust masses from observations. From our study, an average temperature of 45 K (see Table 2) might give more reasonable results; however, further complete studies can help in improving the estimation of T_d.

In conclusion, observations might completely under-estimate the total amount of dust mass contained in disks due to optically thick regions, as it was shown in other studies (Dullemond et al. 2018). Improving our understanding on opacities and disks temperatures are crucial to unveil the mystery around the amount of material available for planet formation.

Fig. 13 Masses of the dust disks for each configuration and each disk parameter derived from the dust evolution model. The squares correspond to the initial amount of dust mass, and the circles show how much of this mass is left after 1 Myr of dust evolution. The differences in initial dust masses for all disks originate from the assumption that the dust-to-gas ratio is initially 0.01 and the inner disk is depleted in dust by the giant planets (see Sect. 3.1). In most cases, all the dust is trapped in the disks, except at high viscosity because the planet gaps are not strong enough to prevent the dust from flowing toward the star. To compare to observed dust masses as in Fig. 14, we also derive the dust contained in the optically thin regions (represented by the markers corresponding to Fig. 14). We seehere that the majority of the mass is hidden in the optically thick regions.

Fig. 14 Observed dust masses for each configuration (colors) and each disk parameter (markers). A black outline surrounds the masses derived using the Andrews et al. (2013) assumptions for the disk’s opacity and temperature. $M_{\text{d},{\tau }_{\nu }<1}^{\text{tot}}$ corresponds to Cols. 5 and 7 (see also Fig. 13) and $M_{\text{d},\text{obs}}^{\lambda =1.3\text{mm}}$ to Cols. 6 and 8 in Table 2. The solid gray line represents the masses for which the observations match the masses from the simulations. The top (bottom) dashed line shows the disks for which the observations overestimate (underestimate) the actual mass contained in the optically thin part by 25%.

5 Conclusions

In this paper, we derived images at λ = 1.3 mm of different planetary system configurations representing the potential Solar System protoplanetary disk. We also derived the images of a giant system composed of three planets of 1 Jupiter mass each, representing a potential initial state for scattering events to happen and produce eccentric planets that match the radial velocity observations (Jurić & Tremaine 2008; Raymond et al. 2009a; Sotiriadis et al. 2017; Bitsch et al. 2020). Using 2D hydrodynamical simulations we determined the gas disk profile in the presence of four (or three) giant planets. This profile was then used as an input for a dust evolution model. After 1 Myr of dust evolution, the resulting dust distributions were used to compute synthetic images of these different disks. Our main conclusions are:

• The dust distributions show that the perturbations created by multiple planets in one disk can lead to substructures that are not directly linked to the positions of the planets. These features are created by traffic jams in the disk, revealing the importance of the gas radial motion in the case of multiple giant planets. Considering that a single planet can also create multiple gaps and rings in a low viscosity disk by perturbing the gas surface density (Dong et al. 2018; Bae & Zhu 2018), this complicates the relation between the number of features created by single or multiple planets and the actual number of perturbers. Our study thus highlights that not all individual gaps and rings are caused by individual planets perturbing the pressure profile of the disks, complicating the link between protoplanetary disk observations and exoplanets;

• By comparing the synthetic images obtained to known observed disks, we showed that the disk phase can be used to derive robust constraints on planet formation scenarios. The presence of bright substructures located at large radii in the DSHARP survey can be explained by the large size, mass, and brightness of these disks. Here, we show that planet formation occurring in smaller disks can easily be missed at low resolutions in the observations (i.e., with a beam larger than 0.04′′ × 0.04′′). One way to improve our understanding of planet formation is thus to observe small protoplanetary disks (i.e., of a few tens of AU) at high resolution to probe the formation environments of different planetary populations;

• The Three-Giants configuration, representing a future system that could experience scattering events after the disk phase, only presents substructures within 40 AU. While Lodato et al. (2019) speculate that the bright rings observed by DSHARP can be explained by the presence of fast migrating giant planets matching the radial distribution of eccentric planets observed by radial velocity, Ndugu et al. (2019) show that this requires a migration at high viscosity, which is contrary to the recent derivation of disk viscosity. Our study here shows that a giant planet system that is susceptible to scatter later during its formation would not produce bright rings in the outer regions during its gas disk phase;

• At highviscosity, too much dust diffuses through the gaps generated by Jupiter and Saturn, inconsistent with terrestrial planet formation (e.g., Izidoro et al. 2021) and cosmo-chemical evidence (e.g., Kruijer et al. 2017). At low viscosity, dust can be retained in a pressure trap exterior to the giant planets, generating large optically thick dust pileups. Self-consistently constraining the dust mass of the disks observationally revealed that the observationally inferred dust mass can be a factor of ten below the real dust mass in optically thick rings in our simulations. Moreover, improving dust temperature estimates can highly improve the estimation of dust from the observations.

This study shows the importance of resolution in observations for our understanding of planet formation. For example, in the Compact configuration (Figs. 4, 5, C.1, and C.2), the features created by the four giant planets were smeared out by the beam of the instrument, making it impossible to determine how many planets are located in this disk. If future surveys focus on very high resolution observations of smaller protoplanetary disks, then it will be possible to distinguish the conditions needed for giant planets to form in the outer or inner regions of the disk. As discussed in Sect. 4.1, an interferometer such as ALMA already has the power to produce images with a high enough resolution. Such observations should be combined with further studies that model the disk structures in the presence of multiple planets. Finally, in order to improve our understanding of the origin of the dust substructures (traffic jams or pressure bumps as discussed in Sect. 4.3), multiwavelength imaging willhelp us determine how many planets are trapped in disks, as well as help us determine the optical properties of the dust. This last point is important for deriving how much mass is available in disks for planet formation.

Acknowledgements

C. Bergez-Casalou and B. Bitsch thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. C. Bergez-Casalou thank K. Dullemond for his help in managing RADMC3D and T. Birnstiel,C. Lenz and M. Garate for their help in managing TWO-POP-PY and F. Cantalloube and J. Milli for their interesting discussions on observations at different wavelengths. P. Pinilla and N.T. Kurtovic are also thankful for the support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. The authors thank the referee for their interesting remarks and questions that helped improve thispaper.

Appendix A Gas hydrodynamical profiles

As discussed in Sect. 2.2, we present in this appendix the outputs of the 2D hydrodynamical setups. In Figs. A.1, A.2, and A.4 we show the perturbed surface densities of the 2D grids of the disks, for the different configurations and disk parameters. Each row represents a configuration, and each column represents a different α viscosity, ranging from 10⁻⁴ to 10⁻² from left to right. The two first figures represent the two different aspect ratios investigated: an MMSN-like aspect ratio in is shown in Fig. A.1 and a smaller aspect ratio, as described in Sect. 2.2, in Fig. A.2). In Fig. A.4 we present the perturbed surface densities for the different masses investigated in the Spread and Three-Giants configurations. In each of the panels of these three figures we see that the gas disk is axisymmetric after t = 12 500 orbits. The vortices triggered by some instabilities or planet growth that could form at low viscosity at the edges of the giant gaps have time to vanish (Hammer et al. 2017; Bergez-Casalou et al. 2020), meaning that we can take the azimuthal average needed as inputs for TWO-POP-PY.

In Figs. A.3 and A.5 we present the azimuthal and time average gas profiles used as inputs for the dust evolution model (Sect. 3.1). The profiles are time-averaged over 2 500 orbits. For each configuration and each viscosity, the profiles at each aspect ratio are plotted in the same panel: the MMSN-like aspect ratio is presented in solid line while the smaller aspect ratio is shown in dashed lines. Each planet’s orbit is represented by a vertical dotted gray line.

These 2D surface densities show the importance of the viscosity. In the Solar System configurations (Figs. A.1 and A.2), Jupiter and Saturn create a common gap at low viscosity whereas only Jupiter is able to start to form a gap at high viscosity. Depending on the planet configuration, at α = 10⁻³, the two inner giants create different features: when they create a common gap in the Compact configuration, some gas is accumulated in between Jupiter and Saturn in the Spread configuration.

In the Compact configuration, at low viscosity and for both aspect ratios, Uranus and Neptune are massive enough to start creating a pileup of gas outside of Neptune’s orbit. It is particularly visible with the small aspect ratio and in the 1D profiles (see Fig. A.3). When we compare these two panels to the two corresponding panels for the Spread configuration, we notice that Uranus and Neptune barely have an effect on the gas disk.

Regarding the planets of different mass and the Three-Giants configuration (Figs. A.4 and A.5), we see that the Three-Giants configuration always creates a deep common gap as the planets are close to one another. However, in the Spread configuration the amount of gas present between Jupiter and Saturn clearly create two different gaps when the planets have reduced masses.

We show in Sect. 3.1 that the velocities of the planets can create traffic jams that produce noticeable substructures. This is due to the fact that the gas disk is highly perturbed by multiple giant planets. In Figs. A.6 and A.7 we present the radial azimuthally and time-averaged profiles used in our simulations. We see that even after averaging the profiles for 2 500 orbits, the disk remain highly perturbed by the planets. In Fig. A.7 we clearly showthat these perturbations are due to the planets as we see that they are stronger for more massive planets.

Fig. A.1 Perturbed gas surface densities (Σ∕Σ0) for an MMSN-like aspect ratio at t = 12 500 orbits of the inner planet, in the Compact (first row) andSpread (second row) configurations. The positions of the planets are marked by dots in each panel (blue corresponds to Jupiter, orange to Saturn, green to Uranus, and red to Neptune). The disks can be considered axisymmetric, which is important for the dust evolution model that takes the 1D radial gas profile as an input.

Fig. A.2 Same as A.2 but for a smaller aspect ratio.

Fig. A.3 Time- and azimuthal-averaged gas surface density profiles for each aspect ratio. Vertical dotted lines represent the positions of each planet. Jupiter and Saturn are the only ones creating substructures in the disks, except in the Compact configuration with a low α and small aspect ratio, where Uranus and Neptune create a small gap and an over-density outside of Neptune’s gap.

Fig. A.4 Perturbed gas surface densities (Σ∕Σ0) at t = 12 500 orbits of the inner planet, in the Spread (first row) and Compact (second row) configurations. The disks have a small aspect ratio and low viscosity (α = 10−4). The masses of the planets are reduced by a factor of two-thirds (left panels) and one-half (middle panels). They can be compared to the total mass configuration in the right panels. As in Fig. A.1, the positions of the planets are marked by dots in each panel and the disks can be considered axisymmetric.

Fig. A.5 Time- and azimuthal-averaged gas surface density profiles for the different configurations. Vertical dotted lines represent the positions of each planet. As expected, more massive planets create deeper gaps. In the Three-Giants case, the giants are close enough to one another to always create a common gap. The masses of the planets then dictate how deep the gap is and how much gas is present in the gap between them. The two orange vertical lines show the positions of the rings seen in the synthetic millimeter images (Fig. 6). We see that they do not correspond to strong features in the gas disk and are located far from the giants’ orbits.

Fig. A.6 Time- and azimuthal-averaged gas radial velocity profiles for the different configurations. Vertical dotted lines represent the positions of each planet. Multiple planets highly perturb the gas velocities, having an important impact on the dust distributions (see Sect. 3.1).

Fig. A.7 Same as Fig. A.6 but for the Spread and Three-Giants configurations and different planet masses. The two orange vertical lines show the positions of the rings observed in Fig. 6.

Appendix B Impact of the radial gas velocity on the dust distributions

In Sect. 3.1 some over-densities are observed at positions that are not directly related to the orbits of the planets or to any perturbations in the gas surface density. These over-densities originate in the radial velocity profile of the gas disk, highly perturbed by the presence of multiple giant planets. These perturbations create traffic jams, where the dust can accumulate without being trapped. In order to determine if these traffic jams are indeed producing such over-densities, we study the dust evolution distribution also with a gas radial velocity forced to be zero.

We take the example of the Three-Giants configuration as it produces the most perturbed disk. We present in Fig. B.1 the dust distributions in the case where the same gas surface density profile is given to the model but the radial velocity profiles are either averaged as in this paper (left panels) or set to zero (right panels). In the first row, we show the integrated dust surface density over all the grain sizes: these profiles allow us to see that the dust is distributed differently in both cases. When the radial gas profile is set to zero, the dust mostly accumulate in the pressure bumps present in the disk. Even if the gas surface density profile present a very slight bump located at 26 AU, creating a small over-density in the dust at this location, it is too small tocreate a noticeable feature in the observations. However, when the gas radial profile is taken into account, thedust gets stuck in these different traffic jams, explaining this spiky behavior. When compared to the positions of the rings observed specially in Fig. 10, represented in orange in this figure, we see that the second ring located at 26 AU is clearly originating in the strong traffic jam located at the same semimajor axis.

Fig. B.1 Dust distributions in the Three-Giants configuration case: in the left panels, the velocity of the gas used as an input for the dust evolution model is the averaged velocity as presented in Fig. A.7; in the right panels, the gas velocity is taken to be null. The first row presents the integrated dust distribution along all the dust sizes, representingthe total dust distribution in the disk, whereas the second row presents the classic dust distributions as in Sect. 3.1. The presence of spikes in the left panel shows that the radial gas velocities are indeed responsible for the dust accumulations in additional rings exterior to the positions of the planets.

As these traffic jams can create noticeable substructures, we conclude that the gas radial velocity profile has a non-negligible impact on the dust distributions when multiple planets are present in the disk. This is important for the derivation of synthetic images but also for dust evolution models.

Appendix C Complementary images

C.1 Solar System images

We show in this appendix the images of the Solar System images corresponding to the normalized intensity profiles presented in Sect. 3.2.

C.2 Images of inclined disks

Here, we present the images of the disks with different inclinations in the Spread and Three-Giants configurations. These images correspond to the radial profiles presented in Fig. 8.

Fig. C.1 Images at λ = 1.3mm for each Solar System configuration for an MMSN-like aspect ratio. These are the images that correspond to the intensity profiles presented in Fig. 4. The beam is 0.04"×0.04" and is represented as the white circle in the lower-left corner of each panel. The white lines represent the distances of the differentplanets.

Fig. C.2 Same as Fig. C.1 but for a smaller aspect ratio.

Fig. C.3 Images at λ =1.3mm in the Spread (first row) and Three-Giants (second row) configuration, at low viscosity and low aspect ratio, for different inclinations. The inclination is increasing from left to right, going from a face-on disk (i = 0°) to a highly inclined disk (i = 60°). The white lines represent the positions of the different planets in each configuration.

References

All Tables

All Figures

Fig. 1 Dust densities distributions for the MMSN-like aspect ratio. Each line represents a Solar System configuration (Compact on the first row and Spread on the second row; cf. Table 1) and each column represents a set of α viscosity and fragmentation speed (increasing from left to right). The vertical dotted lines represent the positions of each planet.The vertical dashed line shows the truncation radius (Eq. (3)). The white lines represent the maximum size reached by the grains, and each line style represents a limiting mechanism: solid line for the fragmentation limit, dashed line for the drift limit, and dotted dashed line for the growth limit size. At high viscosities, dust from the outer regions flows through Jupiter and Saturn’s gap and replenishes the inner disk, which is not the case at low viscosity. At low viscosity, the perturbations induced by the planets in the gas velocity profiles produce dust over-densities (traffic jams), creating substructures not directly related to the positions of the planets.

In the text

Fig. 2 Same as Fig. 1 but for the smaller aspect ratio. A smaller aspect ratio induces deeper planet-induced gaps in the gas disks, creating stronger features in the dust compared to Fig. 1. In the Compact configuration, at α = 10−3, a small amount of small dust flows through the gaps of the giants. This produces a low dust-to-gas ratio in the inner region of the disk, resulting in grain sizes limited by growth and drift rather than fragmentation (see Eq. (7)).

In the text

Fig. 3 Dust densities distributions with the planet masses reduced by one-half and two-thirds in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The masses of the planets mainly change the gap shapes, allowing more or less dust to flow to the inner regions. In the Three-Giants configuration, we also see that more massive planets create more traffic jams in the disk and therefore present more substructures.

In the text

Fig. 4 Radial profile of convolved and unconvolved images at λ = 1.3 mm with the MMSN-like aspect ratio. Each row represents a Solar System configuration, and each column represents an α viscosity. The solid lines show the radial intensity of the images normalized to the peak intensity after convolution with a beam of FWHM = 0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU. The beam is represented with a black horizontal line in the upper-right corner of each panel. The dashed lines represent the normalized intensity of the unconvolved image. Vertical lines show the positions of the planets in each configuration. The light blue area shows the region where the normalized intensity is smaller than Fmin ∕Fpeak, where Fmin = 10 μJy is the minimal flux considered to be observable. The value of Fmin∕Fpeak is different from each panel as Fpeak is unique to each image, whereas Fmin is fixed. Comparing the profiles in the case of the convolved and unconvolved images shows us how many features can be missed due to a too low resolution. Regarding the Compact configuration, all substructures are smeared out in the beam.

In the text

	Fig. 5 Sameas Fig. 4 but for a smaller aspect ratio. As the dust distributions show more intense substructures, some features become observable but the majority are still smeared out in the beam.
In the text

Fig. 6 Radial profile of convolved and unconvolved images at λ = 1.3 mm with the planet masses reduced by one-half and two-thirds in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The intensity profiles are presented as in Fig. 4. The convolved profiles are convolved with a beam of 0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU, represented with a black horizontal line in the upper-right corner of each panel. One can notice the similarities between each convolved profile, making it difficult to disentangle between each evolutionary stage and configuration. In the Three-Giants configuration, we highlight in orange the substructures outside of the giants’ region, originating in the traffic jams discussed in Sect. 3.1.

In the text

Fig. 7 Images in total intensity at λ = 1.3 mm corresponding to the radial profiles presented in Fig. 6. The positions of the different planets are represented by the different white arcs. The beam size (0.04′′ × 0.04′′) is represented in the lower-left corner of each image by the white ellipse. In the Spread configuration, we see that the size of the disks depends on the masses of the planets. In the Three-Giants configuration, one can notice the substructures outside of the giants’ region corresponding to the orange regions of the profiles in Fig. 6 and originating from the traffic jams observed in Sect. 3.1; these traffic jams were created by the perturbed velocity profile of the disk (see Fig. A.7).

In the text

Fig. 8 Radial profile of convolved and unconvolved images at λ = 1.3 mm with different inclinations in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). Here, the disks have the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The convolved images present a beam of 0.04′′ × 0.04′′. The profiles are derived from images that have inclinations ranging from i = 0° to i = 60°. The light blue area shows the region where the normalized intensity is smaller than Fmin ∕Fpeak (see Fig. 4). The only highly impacted region is the inner cavity: more inclined disks hide the inner regions more efficiently, influencing the depth of the observed cavities.

In the text

Fig. 9 Radial profile of convolved and unconvolved images at λ = 1.3 mm with different beam sizes in the Spread configuration (top row) and in the Three-Giants configuration (bottom row). In these simulations, the disk has the smaller aspect ratio and α = 10−4 with vfrag = 1 m s−1. The different beams investigated range from 0.02′′ × 0.02′′ = 2.8 AU × 2.8 AU (left) to 0.1′′ × 0.1′′ = 14 AU × 14 AU (right). We present the fiducial resolution (0.04′′ × 0.04′′ = 5.6 AU × 5.6 AU) in the middle panel for comparison. Each beam is represented by a horizontal line in the upper-right corner of each panel. Orange regions represent the substructures produced by traffic jams, clearly visible in the Three-Giants configuration. The different profiles show that the highest resolution is really needed to start to correctly represent the features of the dust disk.

In the text

Fig. 10 Images in total intensity at λ = 1.3 mm corresponding to the radial profiles presented in Fig. 6. The positions of the different planets are represented by the different white arcs. The different beams investigated range from 0.02′′ × 0.02′′ (left) to 0.1′′ × 0.1′′ (right) and are represented in the lower-left corner of each image by the white ellipses. We show the fiducial resolution (0.04′′ × 0.04′′) in the middle panel for comparison. The colorbars are adjusted for each resolution as they have different sensitivities. These images show the importance of resolution: the substructures start to be well represented at very high resolution.

In the text

Fig. 11 Images at λ = 1.3 mm of the different configurations, at low viscosity and low aspect ratio, at different times: 0.5 Myr, 1 Myr, 2 Myr, and 3 Myr, from left to right. The white lines represent the positions of the different planets in each configuration. We assume that the sensitivity of the instrument is limited to fluxes larger than 10 μJy beam−1. As time evolves, the sizes of the disks shrink due to inward drift, leaving stable rings in the inner disks after 1 Myr.

In the text

	Fig. 12 Radial intensity profiles after 0.5 Myr, 1 Myr, 2 Myr, and 3 Myr of evolution in each configuration with low viscosity and a low aspect ratio. The total size of the disk depends on the time evolution, which determines the dust distribution.
In the text

Fig. 13 Masses of the dust disks for each configuration and each disk parameter derived from the dust evolution model. The squares correspond to the initial amount of dust mass, and the circles show how much of this mass is left after 1 Myr of dust evolution. The differences in initial dust masses for all disks originate from the assumption that the dust-to-gas ratio is initially 0.01 and the inner disk is depleted in dust by the giant planets (see Sect. 3.1). In most cases, all the dust is trapped in the disks, except at high viscosity because the planet gaps are not strong enough to prevent the dust from flowing toward the star. To compare to observed dust masses as in Fig. 14, we also derive the dust contained in the optically thin regions (represented by the markers corresponding to Fig. 14). We seehere that the majority of the mass is hidden in the optically thick regions.

In the text

Fig. 14 Observed dust masses for each configuration (colors) and each disk parameter (markers). A black outline surrounds the masses derived using the Andrews et al. (2013) assumptions for the disk’s opacity and temperature. $M_{\text{d},{\tau }_{\nu }<1}^{\text{tot}}$ corresponds to Cols. 5 and 7 (see also Fig. 13) and $M_{\text{d},\text{obs}}^{\lambda =1.3\text{mm}}$ to Cols. 6 and 8 in Table 2. The solid gray line represents the masses for which the observations match the masses from the simulations. The top (bottom) dashed line shows the disks for which the observations overestimate (underestimate) the actual mass contained in the optically thin part by 25%.

In the text

Fig. A.1 Perturbed gas surface densities (Σ∕Σ0) for an MMSN-like aspect ratio at t = 12 500 orbits of the inner planet, in the Compact (first row) andSpread (second row) configurations. The positions of the planets are marked by dots in each panel (blue corresponds to Jupiter, orange to Saturn, green to Uranus, and red to Neptune). The disks can be considered axisymmetric, which is important for the dust evolution model that takes the 1D radial gas profile as an input.

In the text

	Fig. A.2 Same as A.2 but for a smaller aspect ratio.
In the text

	Fig. A.3 Time- and azimuthal-averaged gas surface density profiles for each aspect ratio. Vertical dotted lines represent the positions of each planet. Jupiter and Saturn are the only ones creating substructures in the disks, except in the Compact configuration with a low α and small aspect ratio, where Uranus and Neptune create a small gap and an over-density outside of Neptune’s gap.
In the text

Fig. A.4 Perturbed gas surface densities (Σ∕Σ0) at t = 12 500 orbits of the inner planet, in the Spread (first row) and Compact (second row) configurations. The disks have a small aspect ratio and low viscosity (α = 10−4). The masses of the planets are reduced by a factor of two-thirds (left panels) and one-half (middle panels). They can be compared to the total mass configuration in the right panels. As in Fig. A.1, the positions of the planets are marked by dots in each panel and the disks can be considered axisymmetric.

In the text

Fig. A.5 Time- and azimuthal-averaged gas surface density profiles for the different configurations. Vertical dotted lines represent the positions of each planet. As expected, more massive planets create deeper gaps. In the Three-Giants case, the giants are close enough to one another to always create a common gap. The masses of the planets then dictate how deep the gap is and how much gas is present in the gap between them. The two orange vertical lines show the positions of the rings seen in the synthetic millimeter images (Fig. 6). We see that they do not correspond to strong features in the gas disk and are located far from the giants’ orbits.

In the text

	Fig. A.6 Time- and azimuthal-averaged gas radial velocity profiles for the different configurations. Vertical dotted lines represent the positions of each planet. Multiple planets highly perturb the gas velocities, having an important impact on the dust distributions (see Sect. 3.1).
In the text

	Fig. A.7 Same as Fig. A.6 but for the Spread and Three-Giants configurations and different planet masses. The two orange vertical lines show the positions of the rings observed in Fig. 6.
In the text

Fig. B.1 Dust distributions in the Three-Giants configuration case: in the left panels, the velocity of the gas used as an input for the dust evolution model is the averaged velocity as presented in Fig. A.7; in the right panels, the gas velocity is taken to be null. The first row presents the integrated dust distribution along all the dust sizes, representingthe total dust distribution in the disk, whereas the second row presents the classic dust distributions as in Sect. 3.1. The presence of spikes in the left panel shows that the radial gas velocities are indeed responsible for the dust accumulations in additional rings exterior to the positions of the planets.

In the text

	Fig. C.1 Images at λ = 1.3mm for each Solar System configuration for an MMSN-like aspect ratio. These are the images that correspond to the intensity profiles presented in Fig. 4. The beam is 0.04"×0.04" and is represented as the white circle in the lower-left corner of each panel. The white lines represent the distances of the differentplanets.
In the text

	Fig. C.2 Same as Fig. C.1 but for a smaller aspect ratio.
In the text

	Fig. C.3 Images at λ =1.3mm in the Spread (first row) and Three-Giants (second row) configuration, at low viscosity and low aspect ratio, for different inclinations. The inclination is increasing from left to right, going from a face-on disk (i = 0°) to a highly inclined disk (i = 60°). The white lines represent the positions of the different planets in each configuration.
In the text