© J. Müller-Horn et al. 2022

1 Introduction

Detailed observations by the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed a high substructure level in the protoplanetary disks around young stars (ALMA Partnership et al. 2015). These substructures frequently present themselves in the form of nearly axisymmetric gaps and rings in the dust continuum emission (Walsh et al. 2014; Andrews et al. 2016; Fedele et al. 2017). A common interpretation of these structures is that they originate from the interaction between the disk and massive planets embedded therein (Long et al. 2018; Huang et al. 2018; Pinilla et al. 2012). A sufficiently massive planet increases the gas-velocity exterior to its orbit to super-Keplerian speeds and in doing so generates a pressure bump in the disk. The pressure bump stops the inward drift of pebbles and the particles accumulate outside the planet’s orbit (Paardekooper & Mellema 2006), thus providing an explanation for the observed ring- and gaplike structures. These models were shown to be capable of reproducing the observed disk emission using hydrodynamical and radiative transfer simulations (Dipierro et al. 2015; Clarke et al. 2018; Nazari et al. 2019).

If the observed rings are in fact induced by gap-opening planets, then we might expect these embedded planets to eventually form the massive planet population detected by the extensive exoplanet surveys of the last years (Wittenmyer et al. 2020; Fulton et al. 2021). However, radial distances of known exoplan-ets lie in the range 0.1−10 AU as opposed to typical radii of the rings in disks that are of on the order of 10−100 AU (Huang et al. 2018). With few exceptions from direct imaging, this region is not yet accessible with current exoplanet detection techniques. Hence, the question arises of whether the yet unseen planets are in fact the progenitors of the discovered exoplanet population. This would require sufficient orbital evolution and inward migration to bring them closer to the star and into the observable regime. From the perspective of planet formation theory, it is therefore interesting to study how these presumed embedded planets would evolve and to compare the emerging population with the one of known exoplanets.

The question of how the suspected embedded planets will evolve has been addressed by Lodato et al. (2019) who studied a large sample of 48 gaps in 32 protoplanetary disks and approximated the masses of the suspected planets from the gap morphology. After evolving the planets for 3−5 Myr, they found a good agreement between the final configuration of the simulated planets and the distribution of cold Jupiters. The authors note, however, that the evolutionary model they used is very simplified and that several aspects could be improved. Caveats include the fact that the disk structure remained fixed throughout the integration and that the same generic lifetimes, disk profiles and viscosities were assumed for all disks in the sample. Furthermore, only the one-planet case has been investigated neglecting possible planet-planet interaction in multiple planet systems. The same is true for a similar study carried out by Ndugu et al. (2019), who pointed out that the disk viscosity is indeed key in shaping the final planet population and that the low viscosities needed to explain the core growth at these large distances are not consistent with inward migration from several 10 AUs to sub 10 AU orbits.

In this paper, we choose an opposite approach where instead of simulating a large sample of disks we investigate the evolution of planets in only one disk, but in turn pay closer attention to the properties of the disk structure as well as considering multiplanet systems. Furthermore, we take care to motivate our initial conditions exclusively from observational constraints.

Our focus lies on the planet candidates of Herbig Ae star HD 163296 whose existence has been surmised both from hydro-dynamical simulations and gas kinematics in the disk. With a radius of more than 500 AU inferred from CO line emission (Isella et al. 2007), the protoplanetary disk around HD 163296 is one of the largest known disks with resolved substructure. It is located at a distance of 101 pc (Gaia Collaboration 2018) and has an estimated system age of approximately 5 Myr (Montesinos et al. 2009). Due to the likely presence of forming planets the disk has been the subject of extensive research in recent years. Observations revealed a high degree of substructure in the spatial distribution of millimeter-sized particles and cold molecular gas. These substructures include the presence of 4 gaps at radii of ~ 10, 48, 86, and 145 AU and a crescent-shaped asymmetry at 55 AU which have been inferred from millimeter continuum and CO molecular line emission (Isella et al. 2016, 2018; Huang et al. 2018). Hydrodynamical simulations point toward a likely planetary origin of these features (Zhang et al. 2018; Liu et al. 2018; Rodenkirch et al. 2021) and recent detections of meridional flows and deviations from Keplerian rotation (Teague et al. 2018; Pinte et al. 2018; Teague et al. 2019) provide further evidence of active planet-disk interaction.

Assuming the gaps in HD 163296 are indeed of planetary origin, the system represents a direct example of the initial stages of planet formation. For the purpose of this paper, we adopt this hypothesis to investigate the future evolution of the planetary system. We would like to point out, however, that this interpretation is not undisputed and several other mechanisms have been suggested to account for the presence of these features, among other things particle trapping at the outer edge of dead zones (Ruge et al. 2016), pebble growth around condensation fronts (Zhang et al. 2015), self-induced dust pile-ups (Gonzalez et al. 2017), as well as zonal flows (Flock et al. 2015).

The paper is organised as follows; in Sect. 2 we describe our disk model, the applied accretion and migration schemes and the setup of simulations. In Sect. 3, we show the results for two exemplary individual systems and present the ensemble statistics of the simulated planet population. We summarize our findings and discuss possible implications in Sects. 4 and 5.

2 Methods

Our simulations were performed using the planet formation code FLINTSTONE. The code models the evolution of planets in a gas-dominated protoplanetary disk thereby taking planetary migration, accretion, the gas damping of eccentricities and inclinations, as well as the evolution of the disk structure into account. The orbital integration was carried out by the MERCURY hybrid symplectic integrator which resolves close encounters (Chambers 1999). In our simulations, we treated collisions between planets as inelastic mergers. A similar version of the code was used by Bitsch et al. (2019) and Izidoro et al. (2021) so we refer to them for a more detailed description and will limit ourselves here to a summary of the main concepts and to highlight any differences.

2.1 Disk Model

The structure of the protoplanetary disk determines the growth and evolution of the planets it contains and is characterized by physical properties such as its geometry, temperature, and density distribution. In contrast to Ndugu et al. (2019) who used disk profiles derived from hydrodynamical simulations and Lodato et al. (2019) who used a generic fixed disk profile for all disks in their sample, here we focused only on one disk, namely the one around HD 163296, and aimed to model its structure based on the available observational data. For this purpose, we followed Kama et al. (2020) who constrained the disk gas mass of HD 163296 using hydrogen deuteride rotational line emission. Their estimate of the total gas mass amounts to M_gas = 0.067 M_⊙. If not indicated otherwise, we used the disk parameters obtained by Kama et al. (2020) which are summarized in Table 1.

Following the viscous accretion formalism (Lynden-Bell & Pringle 1974; Hartmann et al. 1998), the initial distribution of gas in the disk is described by the vertically integrated surface density profile Σ_gas,0 as $${\text{Σ}}_{\text{gas,0}}\left(r\right)={\text{Σ}}_{\text{c}}{\left(\frac{r}{r_{\text{c}}}\right)}^{-\gamma }\exp \left[{\left(-\frac{r}{r_{\text{c}}}\right)}^{2-\gamma }\right]$$(1)

where ${\Sigma }_{\text{c}}=\left(2-\gamma \right)\frac{M_{\text{gas}}}{2\pi r_{\text{c}}^2}$ is the surface density at the characteristic radius r_c = 125 AU and γ is the power-law index. At the characteristic radius r_c = 125 AU this corresponds to an initial surface density value of Σ_gas,0(r = r_c) = 1.29 g cm⁻². Assuming a locally isothermal state the temperature profile can be expressed as $$T\left(r\right)=\frac{\mu m_{\text{p}}}{k_{\text{B}}}\frac{G{M}_{*}}{r}h{\left(r\right)}^2$$(2)

with the proton mass m_p, the gravitational constant G, the mass of the central star M_*, = 1.9 M_⊙ (Rodenkirch et al. 2021), and a mean molecular weight μ which we set to 2.3. In Eq. (2), we introduced the aspect ratio $h\left(r\right)=\frac{H}{r}$ with H being the disk scale height. For the aspect ratio we used a power-law parameterization $$h\left(r\right)=\frac{c_{\text{s}}}{v_{\text{K}}}=h_{\text{c}}{\left(\frac{r}{r_{\text{c}}}\right)}^{\Psi }.$$(3)

In this equation c_s is the sound speed, v_K the Keplerian orbital velocity, h_c the opening angle at the characteristic radius r_c and Ψ is the flaring index of the disk. With Ψ = 0.05, as found by Kama et al. (2020), the disk is slightly less flared than expected from theory; Chiang & Goldreich (1999) derived a flaring index of Ψ_theo = 2/7 for flared disks in radiative equilibrium. A lower flaring index results in a larger torque exerted onto the planets which in turn leads to a shorter migration timescale (see Sect. 2.3). We would like to mention here that a recently published study by Calahan et al. (2021) reinvestigated the thermal structure of HD 163296’s disk and finds updated parameters for its profile. In our notation their findings correspond to parameter values of h_c = 0.087, r_c = 165 AU, M_gas = 0.14 M_⊙ and a flaring index Ψ = 0.08. These results are in good agreement with the values derived by Kama et al. (2020) with the exception of a higher total gas mass, M_gas. We do, however, not expect this difference to influence our results significantly, see Sect. 4 for further details.

In all simulations, the disk extended from r_in = 0.7 AU to r_out = 300 AU where the inner edge of the disk was chosen in accordance with Flock et al. (2019) to account for an inner cavity. As was done in Izidoro et al. (2021), we parameterized the drop in surface density in the inner region of the disk by multiplying the gas density with and additional factor ℜ $$\Re =\text{tahh}\left(\frac{r-r_{\text{in}}}{C{r}_{\text{in}}}\right)$$(4)

with C = 0.05 to allow for a smooth transition. The outer boundary represents the applied ejection distance and was chosen for computational reasons.

Following Shakura & Sunyaev (1973), we used an α disk model to parameterize the turbulent viscosity v of the disk. In this widely used prescription, the parameter α is linked to the effective viscosity of the disk via v = αc_sH. There exist observational constraints on α in the literature but they vary widely in order of magnitude. The distribution of α values inferred from disk observations by ALMA spans a range from 10⁻⁴ to 4 × 10⁻² (Rafikov 2017; Dullemond et al. 2018). For the particular case of the HD 163296 disk, Flaherty et al. (2017) place an upper limit on the turbulent viscosity corresponding to α ≤ 3 × 10⁻³. Recently, Rodenkirch et al. (2021) derived a refined upper bound of α ≤ 2 × 10⁻³ which they find is needed to account for an asymmetric feature in the inner disk. Due to the prevailing uncertainty regarding the precise value, we left α as a free parameter of the study, varying it between 10⁻⁴ and 10⁻³.

The disk structure is, however, not static but evolves on a viscous timescale mainly due to accretion onto the central star. In our simulations, the surface density decayed exponentially on a τ_diss,visc = 3 Myr timescale to mimic gas dissipation (Mamajek 2009). In order to save computation time, the disk structure was recalculated only every 500 yr rather than for every time step of the numerical integrator. Since the disk evolution takes place on Myr timescales this measure does not affect the validity of our results. The total gas lifetime t_disk was treated as a free parameter in the simulations and samples were randomly drawn from an exponential distribution as suggested by Mamajek (2009). We set the lower limit of the gas lifetime to 1 Myr in addition to the current age of the system of t₀ = 5 Myr and cut off the high-end tail of the distribution above t₀ + 7.5 Myr to limit the required computational resources. This setup corresponds to total disk lifetimes ranging from 6 to 12.5 Myr in accordance with the findings of Michel et al. (2021). Since only on the order of 10% of the samples lie above 7.5 Myr we do not expect this limit to have a significant impact on the results. The distribution of gas lifetimes in our simulations is displayed in Fig. 1. Longer disk lifetimes imply more time for the planets to grow and migrate resulting in a more massive and closer-in final population.

During the final stages of the disk lifetime the density becomes low enough for photoevaporation to set in and effectively clear the remainders of the disk on a shorter timescale (Alexander et al. 2014; Owen et al. 2013). In our simulations, this dissipation process was modeled by decreasing the surface density exponentially with a timescale of τ_diss,ph.ev. = 35 Kyr over a period of 300 Kyr from t_diss,ph.ev. = t_disk − 300 Kyr until the end of gas disk life time, as in other N-body simulations (Bitsch et al. 2019; Izidoro et al. 2021). The time-dependent decrease in surface density due to dissipation in our model was, hence, calculated as $$\begin{array}{l}{\text{Σ}}_{\text{gas}}\left(r,t\right)={\text{Σ}}_{\text{gas,0}}\left(r\right)\exp \left[-\left(\frac{t}{{\tau }_{\text{diss,visc}}}\right)\right].\hfill \\ \{\begin{array}{ll}1,\hfill & t<t_{\text{diss,ph.ev.}}\hfill \\ \exp \left[-\left(\frac{t-t_{\text{diss,ph.ev.}}}{{\tau }_{\text{diss,ph.ev.}}}\right)\right],\hfill & t>t_{\text{diss,ph.ev.}}\hfill \end{array}\hfill \end{array}$$(5)

where the time t is set to zero at the beginning of the simulation.

Fig. 1 Distribution of gas disk lifetimes tdisk used in the simulations. The 100 samples (blue histogram) drawn randomly from the underlying exponential distribution (red line) are displayed as well as the cumulative frequency (black outline). The distribution is limited to lifetimes between 1 ≤ tdisk ≤ 7.5 Myr additional to the current age of the system of approximately 5 Myr.

2.2 Gas Accretion

The growth and migration of planets occurs for the duration of the disk lifetime. Planetary cores grow predominantly through the accretion of pebbles which experience an inward drift due to the aerodynamic drag (Brauer et al. 2008). For the more massive planets gas accretion is the dominant growth mechanism. The transition occurs when the growing planet carves a partial gap in the density of the protoplanetary disk. The radial pressure gradient exterior to the orbit of the planet gets reversed and the inward drift of pebbles is halted (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018). Following Bitsch et al. (2018), we express the pebble isolation mass as $$M_{\text{iso}}=25f_{\text{fit}}{M}_{\text{Earth}}$$(6)

where $$f_{\text{fit}}={\left(\frac{h}{0.05}\right)}^3\left(0.34{\left(\frac{\log 0.001}{\log \alpha }\right)}^4+0.66\right)\left(1-\frac{\frac{\partial \ln P}{\partial \ln r}+2.5}{6}\right).$$

It depends on the local properties of the disk, namely the aspect ratio h, the pressure gradient $\frac{\partial \ln P}{\partial \ln r}$, and the α viscosity parameter.

The observed ringlike gaps in the protoplanetary disk of HD 163296 indicate that the inferred planets have already reached pebble isolation mass. Therefore, we assumed that their growth is driven by gas accretion. For the process of rapid gas accretion we followed Machida et al. (2010) who performed three-dimensional hydrodynamical simulations to investigate the mass accretion rate of Jovian planets. Their findings indicate that the effective accretion rate is given by Ṁ_gas = min{Ṁ_gas,low, Ṁ_gas,high} where $${\dot{M}}_{\text{gas,low}}=0.83{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{gas}}{H}^2{\left(\frac{r_H}{H}\right)}^{\frac{9}{2}}$$(7)

and $${\dot{M}}_{\text{gas,high}}=0.14{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{gas}}{H}^2$$(8)

Here r_H denotes the Hill radius of the planet and Ω_K is the Keplerian frequency. However, we do note that there is great uncertainty in the literature about gas accretion rates and suggested values range over several orders of magnitude (Ayliffe & Bate 2009; Machida et al. 2010; D’Angelo & Bodenheimer 2013; Schulik et al. 2019). Ultimately, accretion is limited by what the disk can provide. Lubow & D’Angelo (2006) found that the mass flow rate across the gap opened by a planet lies typically between 10 and 25% of the mass accretion rate outside the orbit of the planet. As was done by Bitsch et al. (2019), we chose to restrict the maximum accretion rate for gap-opening planets to an intermediate value of 80% of the disk’s accretion rate on the star, given by $$\dot{M}=3\pi \alpha h^2{r}^2{\text{Ω}}_{\text{K}}{\text{Σ}}_{\text{gas}}.$$(9)

For a large fraction of planets in our simulations this limit is actually lower than the rates of Eqs. (7) and (8).

For very massive, gap-opening planets the reduction in surface density at the planet’s position may lead to a decrease of the accretion rate onto the planet (Tanigawa & Tanaka 2016) which was not taken into account in our simulations. We consider this choice to be justified since for planets less massive than 10 M_Jup, which make up the vast majority of simulated planets, accretion is not heavily impacted by the formation of gaps (Tanigawa & Tanaka 2016) and gas can traverse the gap even for massive planets (Lubow & D’Angelo 2006). Hydrodynamic simulations indicated, moreover, that massive planets can excite radial motion in the surrounding gas leading to an increase in eccentricity for the outer gap edge (Kley & Dirksen 2006; Li et al. 2021; Tanaka et al. 2022). The onset of gap eccentricity can significantly enhance the feeding zone of the planet and in turn cause larger mass accretion rates. Kley & Dirksen (2006) suggest that this effect becomes considerable for planet to star mass ratios of q > 0.005, which corresponds to planet masses m_pl ⪆ 10 M_Jup in the case of HD 163296. Such massive planets are formed only in rare cases in our simulations, and when they are, it is toward the very end of the disk’s lifetime or through collisions. For this reason we have not considered the effect of gap eccentricity on the mass accretion rates in our simulations but we would like to point out the uncertainty of the mass accretion rates for the most massive planets.

In cases where planets cross the inner edge of the disk at 0.7 AU and the surface density vanishes, accretion rates were set to zero and prevent further growth. Even though in principle gas accretion could continue in the inner region because the stellar accretion rate remains constant, it is unclear whether gas accretion onto planets would proceed normally (Cimerman et al. 2017; Lega et al. 2019; Lambrechts et al. 2019). In light of this uncertainty and the fact that accretion at the inner edge does not play a significant role for the majority of our simulations, we set accretion to zero in this region for definiteness.

2.3 Orbital Migration

Growing planets interact with the protoplanetary disk, thereby exchanging energy and angular momentum. The torque induced on planets by the disk leads to planetary migration and the timescale thereof determines the final configurations of the planetary system. Lower-mass embedded planets are subject to migration in the type-I regime for which we followed the prescription by Paardekooper et al. (2011) who determined a torque formula for nonisothermal type-I migration which has been tested against 3D hydrodynamic simulations (Lega et al. 2015; Bitsch & Kley 2011): $${\text{Γ}}_{\text{tot}}={\text{Δ}}_{\text{C}}{\text{Γ}}_{\text{C}}+{\text{Δ}}_L{\text{Γ}}_{\text{L}}.$$(10)

The total torque is computed as the sum of the Lindblad torque Γ_L which originates from density waves launched at Lindblad resonances and the entropy-related corotation torque Γ_C. The rescaling factors Δ_L and Δ_C were introduced by Izidoro et al. (2017) to account for torque reduction due to the orbital eccentricity and inclination. The torque depends strongly on the local properties and gradients of the disk and scales with the square of the planetary mass. For the specific equations used for implementation we refer to Izidoro et al. (2021; see their Appendix A for the detailed equations). The timescale of type-I migration is given by ${\tau }_{\text{mig,I}}=-\frac{L}{{\Gamma }_{\text{tot}}}$ with L being the angular momentum of the planet.

More massive, gap-opening planets enter the slower type-II migration regime. Classically, a planet carves a gap in the protoplanetary disk if (Crida et al. 2006) $$=\frac{3}{4}\frac{H}{r_{\text{H}}}+\frac{50}{qℛ}\le 1,$$(11)

where q denotes the star to planet mass ratio and ℛ = r²Ω_K/ν is the Reynolds number. Recently, a new approximation for gap opening, especially for the transition between the two regimes, has been put forward by Kanagawa et al. (2018). They relate the timescale of type-II migration τ_mig,II to type-I migration via $${\tau }_{\text{mig,II}}=\frac{{\Sigma }_{\text{un.p}}}{{\Sigma }_{\min }}{\tau }_{\text{mig,I}}.$$(12)

The scaling factor in this equation is given by the ratio between the unperturbed surface Σ_un,p and the surface density at the bottom of the gap Σ_min and can be expressed as (Duffell & MacFadyen 2013; Kanagawa et al. 2015) $$\frac{{\Sigma }_{\text{un.p}}}{{\Sigma }_{\min }}=1+0.04K,K={\left(\frac{m_{\text{pl}}}{M_{\odot }}\right)}^2{\left(\frac{H}{r}\right)}^{-5}{\alpha }^{-1},$$(13)

where m_pl denotes the planet mass. We followed here the approach by Kanagawa et al. (2018). Equations (12) and (13) imply that a higher viscosity of the disk causes the transition between type-I and type-II migration to occur later. As a result, the distance covered by planets through migration is larger.

Finally, the tidal damping of inclinations and eccentricities was implemented using the prescription by Cresswell & Nelson (2006) and Cresswell & Nelson (2008). Specifically, in the type-I regime the eccentricity and inclination damping timescales τ_e and τ_i are given by $$\begin{array}{c}{\tau }_e=\frac{{\tau }_{\text{wave}}}{0.780}[1-0.14{\left(\frac{e}{H/r}\right)}^2+0.06{\left(\frac{e}{H/r}\right)}^3\\ +0.18\left(\frac{e}{H/r}\right){\left(\frac{i}{H/r}\right)}^2],\end{array}$$(14) $$\begin{array}{c}{\tau }_i=\frac{{\tau }_{\text{wave}}}{0.544}[1-0.30{\left(\frac{i}{H/r}\right)}^2+0.24{\left(\frac{i}{H/r}\right)}^3\\ +0.14{\left(\frac{e}{H/r}\right)}^2\left(\frac{i}{H/r}\right)],\end{array}$$(15)

where both formulae include the typical type-I damping timescale τ_wave derived by Tanaka & Ward (2004); $${\tau }_{\text{wave}}=\left(\frac{M_{\odot }}{m_{\text{pl}}}\right)\left(\frac{M_{\odot }}{{\text{Σ}}_{\text{gas.pl}}{a}_{\text{pl}}^2}\right)h_{\text{pl}}^4{\text{Ω}}_{K,\text{pl}}^{-1},$$(16)

and all quantities are evaluated at the location of the planet, as indicated by the label “pl”. When the embedded planet has carved a gap (Eq. (11)), we used a type-II damping prescription, where the eccentricity and inclination damping timescales are fixed to τ_mig,II/K_damp, with K_damp = 5, motivated by the fact that damping is more efficient than migration (Lee & Peale 2002; Bitsch et al. 2013, 2020).

2.4 Simulation Setup

From analyses of dust gaps and gas kinematics in the disk around HD 163296, the presence of multiple planets has been inferred (Liu et al. 2018; Zhang et al. 2018; Teague et al. 2018, 2019; Rodenkirch et al. 2021). We chose to model the evolution of the three planet candidates with the most circumstantial evidence of their existence. According to recent results by Rodenkirch et al. (2021), the inner planet has an approximate mass of m₁ = 1.0 M_Jup and is located at a₁ = 48 AU. For the outer two planets we follow Teague et al. (2018) who found two planet candidates by comparing measurements of rotation curves of CO isotopo-logue emission with hydrodynamic simulations. The derived semi-major axes of the planets are a₂ = 83 AU and a₃ = 137 AU with corresponding masses of m₂ = 1.0 M_jup and m₃ = 1.3 M_Jup. As a reference, the initial surface density values at these positions are ${\Sigma }_{\text{gas,0}}\left(\left\{a_1,a_2,a_3\right\}\right)=\left\{5.86,2.68,1.07\right\}\frac{\text{g}}{{\text{cm}}^2}$. To account for the uncertainties in semi-major axes and masses inferred for these planets we drew a_1,2,3 and m_1,2,3 from Gaussian distributions centered on the aforementioned nominal values. For semi-major axes and masses we chose standard deviations of 10% and 50% of the mean value in agreement with Teague et al. (2018). In doing so, our initial conditions reflect the spread in planet parameters found in other publications (Liu et al. 2018; Zhang et al. 2018; Teague et al. 2019).

As was done in Bitsch et al. (2019), the initial inclinations and eccentricities of the planets were uniformly distributed between 0.01° and 0.05° and 0.001 and 0.01, respectively. All other orbital angles were chosen randomly between 0° and 360°. The evolution of the planetary system and interactions with the surrounding gas were modeled until the end of the gas lifetime, itself drawn at random as mentioned earlier. Subsequently, the resulting system was further integrated without dissipation until 100 Myr to investigate long-term stability. In total, we ran three sets of simulations with increasing α viscosity parameters of α = 1 × 10⁻⁴,3.16 × 10⁻⁴, and 1 × 10⁻³. Each set consisted of 100 individual setups and the sets differed only by the viscosity parameter α of the disk, thus allowing us to study how the final configurations depend on this parameter.

3 Results

In this section, we first present in detail the results of two selected individual systems that are representative of the general phenomenology, before discussing the ensemble statistics of all simulated planetary systems. We then focus on mass, semimajor axis and eccentricity distributions of the simulated planet population and compare them to the distributions of known exoplanets. To allow for a fair comparison, we only include exoplanets around stars with luminosities exceeding 5L_⊙ in the statistics, corresponding to the luminosity limit of A-type stars and an approximate lower mass boundary of 1.5 M_⊙. This way we take the observed correlation of giant planet occurrence rates with stellar mass (Johnson et al. 2010) into account and ensure that the host stars have similarly sized inner dead zones as HD 163296 (Flock et al. 2019).

3.1 Individual Systems

It can be insightful to examine the evolution of individual systems in order to interpret the over-all statistics of all simulations. Therefore, in this section we illustrate the effects of orbital migration and gas accretion as well as resonant capture using the example of three selected systems.

Figure 2 displays the growth and dynamical evolution of three synthetic systems from here on referred to as systems I, II, and III. In all systems α is equal to 10⁻³. All planets migrate inward where the rate of migration depends on the planetary mass and the aspect ratio of the disk at the position of the planet, see Eqs. (12) and (13). Migration stops either with the arrival of the planet at the migration trap at the inner edge of the disk, or at the end of disk life time, which for systems I, II and III occurs at 1.56, 1.34 and 6.59 Myr, respectively. In all three simulations the initial mass of planet 2 is smaller than for planet 1. Hence, planet 2 migrates faster and the period ratio of the two planets decreases over time, a process known as convergent migration. As the orbits approach each other, the gravitational interaction between the planets drives an increase of their eccentricities (not directly shown but indicated by the shaded region around the semi-major axes which extends form pericenter (a_pl(1 − e)) to apocenter (a_pl (1 + e))).

All the while, the planets grow more massive by accreting gas from the surrounding disk. In the beginning, the rate of gas accretion onto the planets increases because the planets migrate inward toward higher surface densities. At later stages during the disk lifetime, the exponential decrease in surface density becomes noticeable in that the accretion process slows down, which can be seen in the mass evolution of system II. Interestingly, in all exemplary systems the inner two planets enter a 2:1 mean motion resonance within the first Myr. This can be seen from the libration of the mean motion resonance angles ${\psi }_1^1$ and ${\psi }_2^1$which for a general k:(k − 1) mean motion resonance are defined as $${\psi }_1^1=k{\lambda }_2-\left(k-1\right){\lambda }_1-{\varpi }_1$$(17) $${\psi }_2^1=k{\lambda }_2-\left(k-1\right){\lambda }_1-{\varpi }_2.$$(18)

Here, ϖ_i denotes the longitude of the pericenter (which is given by the sum of the longitude of the ascending node and the argument of pericenter of planet i ∈ 1,2,3), λ_i is the mean longitude of the orbit of planet i (defined as the sum of the longitude of pericenter and the mean anomaly) and k is the resonance index, which equals 2 in the case of a 2:1 resonance. In the case of a mean motion resonance the resonance angles librate around a constant value.

Mean motion resonances are quite common among all simulations (see Sect. 3.4) and can occur both for the inner and outer planet pair. In fact, in system I, not only the inner planet pair experiences resonant motion but the outer pair is locked in a 2:1 mean motion resonance as well. This resonant chain is what prevents the closely spaced system from becoming unstable. In contrast, in system II, the inward migration of planet 3 disturbs the inner two planets which forces the pair out of mean motion resonance and results in the collision of planet 1 and 2 before the end of disk lifetime. However, the rare occurrences of resonant chains are not the only chance for the inner planet pair to remain in a stable configuration. Instead, more often than not the outer planet has a sufficiently low migration rate in order not to migrate to regions where it could disturb the inner planets. This is the case for the majority of simulations with α ∈ {1 × 10⁻⁴,3.16 × 10⁻⁴) and for comparably high initial masses of planet 3. System III shows an example of this behavior.

Fig. 2 Evolution of three selected planetary systems I (top), II (center), and III (bottom) over the integration time of 100 Myr. The plots display the development over time of the semi-major axes ai (left), the planetary masses mi (center) and the resonance angles ${\psi }_1^1$, ${\psi }_2^1$ of the inner planet pair (right). The disks in all three systems have α = 10−3. System I and III have gas disk lifetimes of l.56Myr and 1.34Myr, respectively, and remain stable throughout the integration whereas system II has a gas disk lifetime of 6.59 Myr and becomes unstable at 2.36 Myr when the inner two planets collide and merge.

3.2 Planetary Masses and Orbital Distances

We now turn to the full sample of our synthetic systems. Of particular interest to this study are the final locations of the simulated planets in the mass vs. semi-major-axis plane, especially in relation to the known exoplanet population. Figure 3 shows a superposition of all 300 individual systems. We plot the final locations of all planets that survived the 100 Myr integration time. The lines across the planet span the range from a_pl(l − e) to a_pl(1 + e) with e the mean eccentricity averaged over the final 10 Myr.

We note that the observed gaps in HD 163296 indicate that the initial planet masses should exceed the pebble isolation mass at the respective positions. For reference, in Fig. 3 we show the estimated pebble isolation mass (Eq. (6)) for α = 10⁻³ as a dotted line. The initial conditions of all simulations are indicated in Fig. 3 as red dots. Almost all simulations start out above M_iso and the few exceptions reach pebble isolation mass within the first Myr of the integration.

From there, the same evolutionary pattern is evident for all simulations, with rapid gas accretion and orbital migration leading to mass growth and an overall decrease in semi-major axis, visible as a leftward and upward motion in the m−a-plane. Despite this general pattern, there is a large variance among the final positions and the mass and semi-major axis values cover a range of more than two orders of magnitude. The final masses and semi-major axes correlate to some extent with the gas lifetime of the disk, visualized in Fig. 4, as well as with the initial positions of the planets. Figure 4 shows that a disk lifetime of well over 2 Myr in addition to the current system age is required for the formation of planets with masses greater than 10M_Jup. This is in rough agreement with the results of Schlaufman (2018) who suggest that for planets with masses >10M_Jup gravitational instabilities might be the predominant formation pathway. The dependencies of the final configurations on disk lifetime and initial positions are to be expected because for a given mass a closer initial position of a planet with respect to the star results in a smaller final orbital distance and longer gas lifetimes allow more time for the planets to migrate inward (see also Fig. 2). However, it turns out that the crucial factor in determining the end-positions of the systems is the α viscosity parameter of the underlying disk structure. It is apparent from Fig. 3 that the simulations with α = 10⁻³ constitute the most massive and close-in population. A clustering of planets near 0.7 AU can be observed due to the migration trap at the inner edge of the disk. As a consequence of planet-planet interaction, several planets among the α = 10⁻³ set have found their way even further in, to radii as small as 0.3 AU. To compare with the results of Bitsch et al. (2020), for this choice for α the mean planetary mass is (6.8 ± 4.2)M_Jup, the mean semi-major axis is (6 ± 12) AU and the mean separation of planets in mutual Hill radii is 8.7 ± 4.5.

In contrast, the simulations with the lowest viscosity, corresponding to α = 10⁻⁴, have lower migration and accretion rates. The planets in this set have final masses generally below 3 M_Jup with a mean mass of (1.4 ± 0.5) M_Jup. The majority remains outside the inner region of 10 AU and the mean semi-major axis is (55 ± 40) AU. In this case, the mean separation in mutual Hill radii is with 11.1 ± 5.5 much larger than for the high viscosity simulations which contributes to the comparably low rate of planet-planet interaction in this set, see Sect. 3.3.

For the ensemble with intermediate viscosity, α = 3.16 × 10⁻⁴, a large fraction of planets remain outside 10 AU but the set also contains several systems in which the inner planet migrates up to the disk edge, thus linking the high- and low-α simulations nearly seamlessly in terms of semi-major axis. The masses of the planets in this set are in between those of the high- and low-α set with mean mass and semi-major axis values of (2.4 ± 0.9) M_Jup and (26 ± 31) AU. The mean separation of the planets in this set is 12.0 ± 5.7 in units of mutual Hill radii.

To put our results into perspective, we analyze the synthesized giant planet population in the context of the known exoplanets. To this end we use the exoplanet data provided by the NASA Exoplanet Archive as of 03 November 2021. We note that we are not aiming at directly linking the current exoplanets census with our synthetic population of planets, since we are dealing with a single system whose intrinsic properties cannot be generalized to the full sample. Instead, we use the known exoplanet catalog as a gauge to compare our results to the orbital features of known exoplanets. Out of the sample of known exoplanets we select those for comparison to our simulations that have known m_pl sin (i) values from radial velocity (RV) measurements and orbit stars comparable to HD 163296 with luminosities exceeding 5 L_⊙. For the purpose of this paper we consider all simulated planets as detectable that induce RV amplitudes of more than 1 m s⁻¹ and have orbital periods P < 10 yr. This results in a division of the m-a plane into an observable and an RV unobservable region as visualized in Fig. 3 by a dashed black line.

Several conclusions can be drawn from the distribution of semi-major axes and masses of the simulated planets. If we assume that embedded planets such as the presumed planets in the HD 163296 disk evolve to form part of the observed exoplanet population, then higher viscosities with α ≳ 3.16 × 10⁻⁴ are likely required for sufficiently high migration rates. The lower viscosity simulations with α = 10⁻⁴ seem not to produce enough close-in planets to account for the abundance of giant planets detected with RV measurements. Increased sensitivity of RV and direct imaging observations could help to further constrain the value of α, under the assumption that the initial population of wide orbit giants is the norm rather than the exception. If future surveys will uncover a larger number of planets in the region of m-a parameter space with m_pl ≲ 3M_Jup and a_pl≳ 10 AU, this would point toward an intermediate α viscosity parameter of α < 10⁻³. Instead, if Jupiter-mass planets with semi-major axes of several tens of AU will be ruled out by observations, then our results would imply that the viscosity should be even higher, that is α ≈ 10⁻³.

It is apparent in Fig. 3 that we do not fully reproduce the observed distributions of exoplanet parameters. For example, in the observed population of exoplanets, there is a cluster of planets with semi-major axes between 1 and 3 AU and masses of 1 to 3 M_Jup which does not appear in the simulated planet population regardless of the value of the viscous alpha parameter. However, we would like to emphasize that we do not expect to reproduce the observed distribution because our simulations are based on the example of a single disk, albeit with varied initial conditions and disk parameters. Instead, we make statements about the choice for the disk parameter values for which the simulations fit into the general trends of the observed planets.

The best-studied and classic example for direct imaging campaigns to date is the planetary system of the A5V star HR 8799, which provides a benchmark for planet formation models and is one of few discovered systems with multiple directly imaged planets (Marois et al. 2008, 2010; Wang et al. 2018). It is, therefore, interesting to understand whether the planets of HR 8799 could have formed from rings and gaps as we simulate in our work. When trying to answer this question, we have to keep in mind that the results of these simulations rely on the idea that the gaps and rings observed in protoplanetary disks are caused by planets, an assumption that is still under debate.

The approximate positions of the four directly imaged planets of HR 8799 in the m-α-plane are shown in Fig. 3 in relation to our simulated planet population (Wang et al. 2018). The semi-major axes of the four planets cover a range typical of our simulations with low viscosity, that is 10⁻⁴ ≤ α ≤ 3.16× 10⁻⁴. Themasses of the planets, derived from luminosity and dynamical constraints, are, however, more in line with our high viscosity simulations with α = 10⁻³. The apparent incompatibility between masses and semi-major axes could indicate a mismatch between migration and accretion rates but could also be explained if the planets formed on wider orbits. Alternatively, trapping mechanisms located outside the inner edge could have halted migration. For example, a large number of magneto-hydrodynamic simulations have shown that magnetorotational turbulence in disks can evoke radial structuring in the form of zonal flows (e.g., Balbus & Hawley 1991; Flock et al. 2015). Coleman & Nelson (2016) have demonstrated that a radial structuring of the disk can in turn lead to the formation of persistent planet traps at large orbital radii. It is also conceivable that the local gas surface density of the protoplanetary disk might have been higher which would have led to higher accretion rates. Since the high masses and large orbits of the HR 8799 planetary system still pose a challenge for core accretion models, gravitational instabilities of the gas disk cannot be ruled out as an alternative formation mechanism either (Boss 2011). We highlight again that these comparisons are based on the underlying assumption that the HD 163296 disk is representative of the full disk sample which is not known with certainty.

Fig. 3 Mass-distance plot comparing all simulated planets to the known exoplanet population. Small red dots indicate the initial values of the simulations. The pebble isolation mass is shown with a dotted line for α = 10−3. Blue, orange and green dots denote the planets from simulations with α viscosity parameters 10−4, 3.16 × 10−4, and 10−3, respectively. The x-bars extend from the perihelion to the aphelion distances for a given planet. Gray dots denote the exoplanet population detected via radial velocity and black circles mark those around massive stars with luminosities exceeding 5 L⊙. For reference, the planets of HR 8799 are shown as black dots (Wang et al. 2018). A dashed line divides the parameter space into a region observable by current radial velocity measurements with a sensitivity of 1 m s−1 and observation period of up to 10 yr and a currently not observable region. The histograms depict the normalized cumulative frequencies of the final masses (right) and semi-major axes (top). The color scheme is the same as for the main plot, where the bars show the distribution among all simulated planets and the dashed outlines represent the distribution among hypothetically observable planets. The cumulative frequencies among the exoplanets of massive stars with L > 5 L⊙ are plotted for comparison with a black solid outline.

Fig. 4 Mass vs. semi-major axis plot comparing all simulated planets to the known exoplanet population. The plot is in all aspects the same as Fig. 3, except that the simulated planets are colored to reflect the gas lifetime tdisk of the disk in which they formed instead of the α viscosity parameter. Simulated planets with tdisk < 2 Myr, tdisk between 2 and 5 Myr, and tdisk ≥ 5 Myr are shown as colored dots in blue, turquoise and pink, respectively. It is apparent that the dependence of the final orbital parameters on the disk lifetime tdisk is small compared to the dependence on the viscosity parameter α (see Fig. 3).

3.3 Multiplicities and Eccentricity Distribution

For the duration of the disk lifetime the orbital eccentricities and inclinations are tidally damped (see Sect. 2.3). With time, as the disk evolves, the eccentricity and inclination damping become less efficient because of the exponential decrease in surface density. After the dispersion of the disk during the final stages of disk lifetime, the dampening of eccentricities and inclinations ceases altogether. On the other hand, eccentricities and inclinations may increase due to the interaction between multiple planets and the exchange of angular momentum between them. For the observed exoplanets, the eccentricity distribution was found to be consistent with being an outcome of planet-planet scattering (for a review see Davies et al. 2014 and references therein). However, in the absence of other sources to disturb the system, these simulations require a minimum number of three planets to trigger sufficient scattering events (Weidenschilling & Marzari 1996; Jurić & Tremaine 2008; Raymond et al. 2009; Sotiriadis et al. 2017).

These mechanisms – the tidal damping and the effect of planet scattering – shape the cumulative distribution function of eccentricities of our simulated planets, which we show in Fig. for different α parameters. The planetary systems with lower α values typically remain in stable configurations where the planetary orbits are widely separated and no close encounters occur. Accordingly, the eccentricities are comparably low, with a mean value of 0.028 for α = 10⁻⁴.

For higher viscosities the migration rates increase causing more planets to migrate to the inner edge. Therefore, close encounters occur more frequently and the gravitational interactions among the planets lead to mutual excitation of the eccentricities. This is evident in a significant increase of planet eccentricities with increasing α, with mean values of 0.052 and 0.076 for simulations with α = 3.16 × 10⁻⁴ and 10⁻³, respectively. For comparison, the mean of the observed eccentricities among the detected exoplanets, shown in Fig. 5 to be around 0.135, is even higher than for our α = 10⁻³ set.

The moderate eccentricities among the simulated population could be a result of the relatively low number and large initial separations of planets in the disk. A higher mean eccentricity could be achieved with additional planets in the synthetic systems which is not unreasonable, considering the growing evidence for a fourth planet in the HD 163296 disk (Pinte et al. 2018; Teague et al. 2019, 2021). Furthermore, eccentricities could be increased via external perturbations (for instance, by planetes-imal belts such as those invoked in the Nice model to drive the giant planet instability in the Solar System, (Raymond et al. 2010)) which can act as catalysts for dynamical excitation. In addition, observations of exoplanets might be slightly biased toward higher eccentricity values. More precisely, recent observations (Anglada-Escudé et al. 2010; Wittenmyer et al. 2013, 2019; Kürster et al. 2015; Boisvert et al. 2018; Hara et al. 2019) indicate that in cases of sparse data, two giant planets in circular orbits can mimic the radial velocity signature of a single eccentric giant planet, thus, leading to potential misinterpretations of the eccentricity distribution.

For the exoplanet population it was shown that eccentricities correlate inversely with system multiplicities (Bach-Møller & Jørgensen 2021). The same is true for the simulated set of planets, as is apparent in the lower panel of Fig. 5. In the 2-planet systems the eccentricities of the remaining planets are driven toward higher values simply because these are the systems in which strong mutual interaction occurred which lead to dynamical instabilities and scattering events. Therefore, the 2-planet systems have generally higher eccentricities than their 3-planet counterparts. Given that in the type-II regime the migration timescale scales inversely with the α viscosity parameter, it is not surprising that with 4, 10, and 76 ejections and collisions for the α = 1 × 10⁻⁴,3.16 × 10⁻⁴,1 × 10⁻³ sets, respectively, the majority of scattering events occurs for the set with highest α value, α = 1 × 10⁻³.

The above ejection and collision rates correspond to an average of 2.96, 2.9, and 2.24 surviving planets per system for the α = 1 × 10⁻⁴,3.16× 10⁻⁴,1 × 10⁻³ sets, respectively. Considering only the hypothetically observable planets, these values are reduced to system multiplicities of 0.16, 1.04 and 1.71. The apparent discrepancy between observed and actual multiplicities implies that in systems observed to harbor one giant planet, several others might still be hidden. This is in line with findings by Bitsch et al. (2020) whose simulations indicate that giant planets should, on average, occur in multiples rather than singly.

Fig. 5 Cumulative histogram of the eccentricity distribution among the simulated planet population. In the upper panel, the normalized distributions of simulations with α viscosity values 10−4, 3.16 × 10−4, and 10−3, are depicted in blue, orange and green, respectively. Dashed lines indicate the eccentricities among all simulated planets, whereas the eccentricity distributions of only the hypothetically observable planets are shown as solid lines. The lower panel confronts the eccentricity distribution of planets in stable 3-planet systems (purple) with 2-planet systems (red) in which one planet collided or got ejected. For comparison, the eccentricity distribution of observed exoplanets around stars with luminosities exceeding 5 L⊙ is shown as a black solid line in both panels.

3.4 Period Ratios

The period ratios in multiple planet systems reflect the spacing of neighboring planets and the dynamical configuration of the system. In Fig. 6, we show the period ratios of adjacent planet pairs from our simulations. The normalized histograms in the top panel depict the cumulative distribution of period ratios for simulations with different α viscosity parameters. Dashed lines indicate the distribution among all simulated planets, the distributions among hypothetically observable planet pairs are shown as solid lines. Since lower α viscosity values imply slower migration, there are fewer planet pairs with α = 3.16 × 10⁻⁴ in the observable region than for α = 1.0 × 10⁻³. This can also be seen from Fig. 7 where we show the combined mass of the planet pairs as a function of their period ratio. For the lowest α value, only one planet pair enters the observable region.

Mean motion resonances occur frequently among all our simulations and can arise both for the inner and outer planet pair. Independent of the viscosity parameter, the distributions show an accumulation of period ratios near the 2:1 and 3:1 mean motion resonances. The majority of 2:1 resonances arises for the inner planet pair. In these cases, the second planet approaches the inner one either because it is smaller and migrates faster or because the migration of the innermost planet is stopped at the inner edge of the disk while the second planet continues to migrate further in. This resonance can act as a stabilizing factor in cases of closely packed systems. If, however, no resonance occurs or the resonance is disturbed by the third planet, as is the case for system II in Sect. 3.1, the system becomes unstable and at least one planet collides with another planet or the star or it gets ejected.

In the lower panel of Fig. 6, the difference in period ratios between 3-planet and 2-planet systems is shown. Interestingly, there is a pile-up of period ratios near the 3:1 resonance for the remaining planet pair in unstable 2-planet systems.

The clustering of period ratios near the 2:1 and 3:1 resonances that is apparent in Fig. 7 is interesting also in the light of recent findings by Bitsch et al. (2020), whose simulations suggest that mean motion resonances should be rare among giant planets. Bitsch et al. (2020) argue that the lack of resonances might be due to a low α = 10⁻⁴ parameter and hence slow migration rates. The fact that resonances occur frequently in our simulations indicates that different formation pathways can result in different system architectures. For instance, the simulation set up used by Bitsch et al. (2020) differs from ours in that the number of planetary embryos is much higher, ranging from 15 to 60. At the same time, their systems are more compact and the outer embryos are located at a maximum distance of 17 AU. Both factors lead to an increase in the number of scattering events in comparison to our simulations.

A direct comparison with observations is currently not feasible because there are too few detected giant planet pairs around luminous stars to draw meaningful conclusions. Moreover, given that the diversity of disks and initial solid distribution is expected to yield different system configuration, we cannot draw strong conclusions from our single-disk simulations. However, considering all known multiplanet systems, there seems to be a clustering of period commensurabilites for planet pairs with a total mass exceeding 1 M_jup (Wright et al. 2011; Winn & Fabrycky 2015) which would agree with our results. Currently, we are limited to the example of HD 163296 but future discoveries of other protoplanetary disks with strong evidence of multiple planets will help significantly to understand whether the simulated processes can fully reproduce the observed eccentricity and period distributions.

Fig. 6 Normalized cumulative distribution function of the period ratios among simulated planet pairs. In the upper panel, the period ratio distributions of simulations with α viscosity values 10−4, 3.16 × 10−4, and 10−3, are depicted in blue, orange and green, respectively. Dashed lines indicate the period ratios among all adjacent planet pairs, whereas period ratios of those planet pairs that would hypothetically be observable are shown as solid lines. The lower panel compares the period ratio distribution of planets in 2- (purple) and 3-planet systems (red). We observe a significant pile-up of period ratios near the 2:1 resonance for 3-planet systems and around the 3:1 resonance for 2-planet systems.

4 Discussion

In this section, we discuss the underlying assumptions of our modeling and outline some of the implications of our results for exoplanet and protoplanetary disk observations.

4.1 Effect of Changes to the Disk Model

The outcome of our simulations is largely determined by the α viscosity parameter of the underlying disk. Since most relevant quantities, such as the migration and accretion rates scale with α, see Eq. (9), and the α parameter values of the simulations span an order of magnitude, we do not expect qualitative changes in our findings from changes in the disk structure by a factor less than that (such as a factor 2 in surface density as proposed by Calahan et al. 2021).

In this work, we limited ourselves to disk models with constant α viscosity parameter as these are the most commonly applied. The effects of a varying α parameter as suggested by Liu et al. (2018), for example, and applied by Rodenkirch et al. (2021) remain to be tested but are not expected to cause significant changes in our results. In the case of the profile used by Rodenkirch et al. (2021), for example, the α parameter at the location of the innermost planet is set to α(48 AU) ≈ 2 × 10⁻⁴, meaning that the inner planet would likely evolve in a similar manner as the inner planets in our α = 1 × 10⁻⁴ simulations. The increase of α with radius would cause the outer planets to migrate and grow faster than in our simulations until they reach the inner parts of the disk. Therefore, a probable result of these changes would be a slightly more closely spaced final planet population. The decrease in relative distances would likely lead to an increased level of planet-planet interaction and somewhat higher planet scattering rates with respect to our α = 1 × 10⁻⁴ simulations.

4.2 Constraints for Protoplanetary Disk Parameters

Our simulations are based on the premise that the supposed planets in HD 163296 actually exist. If, in addition, we assume that such massive planets on wide orbits are typical of young exoplanet systems, at least for host stars of similar sizes, this enables us to draw conclusions about the parameters of protoplanetary disks. More precisely, if the putative planets of HD 163296 are common examples of young exoplanet systems they should evolve so that their final planetary and orbital parameters match those of the discovered exoplanet population. Our models show that the agreement of the final mass and semi-axis distributions of the simulated planets with the observations depends crucially on the viscosity parameter α of the disk and requires that α ≳ 10⁻⁴. From the masses and semi-major axes of the synthesized planets we constrain the viscosity parameter to an approximate range of 3.16 × 10⁻⁴ ⪅α⪅ 10⁻³. A a parameter toward the lower end of this range is expected to cause the formation of a larger number of planets with masses below m_pl ≲ 3M_Jup and orbits outside a_pl ≳ 10 AU whereas α ≈ 10⁻³ results in a more massive and close-in population. This prediction provides a way to further constrain the value of α with future observations.

Independent of the choice for α, many planets remain outside the region in m-a parameter space that is accessible to RV observations at current measurement sensitivities. Therefore, the multiplicities of the simulated systems that would hypothetically be observed are far lower than the actual system numbers. Still under the assumption that the HD 163296 planets are common examples of young and massive exoplanets, this result suggests that exoplanet multiplicities are likely being underestimated and that there might be hidden companions in systems known to harbor at least one giant planet. Indications of a potential underestimation of multiplicities due to the RV detection limit were also found in previous studies by Wagner et al. (2019) and Bitsch et al. (2020), among others.

Future direct imaging campaigns are a promising alternative to RV measurements for detecting planets in orbits beyond 10 AU. With the exception of HR 8799, direct imaging has so far mostly discovered systems with no or single planets. If this trend continues in future observations, it would suggest larger values of α on the order of α ≈ 10⁻³, because for higher viscosity the giants should migrate to the inner disk regions where they would no longer be detectable by direct imaging. On the other hand, if the viscosities are low, that is 1 × 10⁻⁴ ⪅α⪅ 3.16 × 10⁻⁴, a larger number of Jupiter-mass planets should be detected at distances of a > 10 AU. As for the number of planets per system, the observed low number could also indicate that planets forming in rings similar to our simulations should mostly have only one survivor. These results, of course, all rely on the assumption that the observed gap and ringlike structures in protoplanetary disks are in fact caused by planets, which to date is not entirely clear.

Fig. 7 Final period ratios of all simulated planet pairs plotted against their combined mass. Crosses mark planet pairs in 3-planet systems and circles denote 2-planet systems, where the third planet collided or got ejected. The transparency of the markers indicates whether or not the pair would be observable with current standard radial velocity measurements. The different colors denote the α viscosity value of the protoplanetary disk with blue, orange and green for α = 10−4, 3.16 × 10−4, and 10−3, respectively. The period ratios of known exoplanet pairs around stars with luminosities exceeding 5 L⊙ are marked with black dots. The horn-shaped black solid curves indicate the approximate resonance boundaries $\left(p_{out}/P_{\text{in}}-p/q\right)=0.05{\left(\left(M_1+M_2\right)/M_{\text{Jup}}\right)}^{1/2}$ for a p:q mean motion resonance, in this case for 3:2, 2:1, and 3:1 resonances. It is evident that period ratios close to the 2:1 mean motion resonance occur frequently for all α viscosity values. But since for the lower two α viscosity values the distances to the star are typically larger only few of the resonant pairs would be observable.

4.3 Implications of Mean Motion Resonances

In contrast to Bitsch et al. (2020) whose findings indicate that mean motion resonances among giant planets should be rare, (near)resonant configurations of planets occurred frequently in our simulations. We conclude that the cause of this discrepancy is likely to be the difference in formation pathways, see Sect. 3.4. A conceivable implication for observed exoplanet systems is that those systems exhibiting resonances originated from disks with fewer and or more widely spaced planetary embryos. In this way, the chance of scattering events would have been reduced and newly formed resonant configurations would have been less likely to be perturbed by additional scattered or migrating bodies. However, our simulations presented here are not able to reproduce the eccentricity distribution of the giant planets, indicating that a lower number of initial embryos might be responsible for forming systems in resonance, while systems with many more initial embryos could explain the eccentricity distribution of the giant planets. This actually indicates that there are maybe different planet formation pathways for giant planet systems.

In our simulations, there appears to be a connection between the number of hypothetically observable mean motion resonances among the planets and the underlying viscosity parameter α of the simulated disk. Under our assumptions for the formation pathway detailed above, it is therefore conceivable to use the abundance of observed exoplanet pairs in mean motion resonance to narrow down the magnitude of the α viscosity parameter in disks such as the HD 163296 protoplanetary disk. Since for lower α values the distances covered by migration are smaller, fewer exoplanet pairs will exist for which both planets are observable. This implies that fewer resonant systems would be observed whereas the opposite would be true for high α values. A step in this direction has already been taken by Hühn et al. (2021), who derived information about the viscosity and surface density of the former protoplanetary disk from the observed orbital architecture of the Kepler-223 system.

4.4 Location of Rings

For higher α parameters and, as a result, increased migration rates, one would expect the ring positions to evolve in time as the planets migrate inward, because the ring caused by the planet should follow the planet’s inward migration. For similarly sized disks, the location of the outermost ring should, therefore, be a function of the disk age, still under the assumption that all rings are caused by forming planets. So far, no such correlation has been observed (van der Marel et al. 2019). This could indicate lower viscosities but might also be a consequence of the relatively small number of resolved disks in the sample. Different migration rates should in principle also be reflected in the morphology of the gaps and rings induced by the planet (Meru et al. 2019; Weber et al. 2019). However, using this approach to identify migrating planets has proven to be difficult and has not yet been tried for the disk of HD 163296 (Nazari et al. 2019).

5 Summary and Conclusions

Detection of planets in their early stages of evolution is crucial to test and improve current planet formation theories and to constrain the initial conditions of population synthesis models. To date, however, only one system harboring young planets and a protoplanetary disk has been detected by direct imaging (Keppler et al. 2018) and for few other stars the presence of such embedded planets has been put forth.

In this paper, we investigated one of these rare cases, namely the planet candidates around HD 163296, inferred from hydro-dynamical simulations and the observations of gas kinematics in the protoplanetary disk (Isella et al. 2016, 2018; Huang et al. 2018; Zhang et al. 2018; Liu et al. 2018; Teague et al. 2018, Teague et al. 2019; Pinte et al. 2018; Rodenkirch et al. 2021). We performed N-body simulations that model the evolution and growth of the three putative planets and their interactions with the surrounding disk. In contrast to previous work on this topic, we considered multiplanet systems and the effects of gravitational interaction between planets. Furthermore, we took the evolution of the underlying disk structure including a likely range of disk lifetimes into account. We motivated our simulations from the available observational constraints of the planet and disk parameters.

We find that the final masses and semi-major axes of the simulated population are largely determined by the α viscosity parameter of the underlying disk structure, for which we adopted three different values (α ∈ {1 × 10⁻⁴,3.16 × 10⁻⁴,1 × 10⁻³}). A higher α parameter results in higher migration and accretion rates of the embedded planets, leading to a more massive and close-in population with a clustering of planets near the inner edge of the disk for α = 10⁻³. In contrast, a lower viscosity parameter, α = 10⁻⁴, results in a less massive final population where the simulated planets remain mostly outside the inner region of a ≲ 10 AU and generally have masses below 3 M_Jup. An intermediate value of α ~ 3 × 10⁻⁴ results in a blend of the two outcomes in terms of semi-major axis, while reaching intermediate masses. Only to a lesser extent do the final parameters correlate with the gas life-time of the disk or the initial positions of the planets.

If young planets with wide orbits, such as the putative planets of HD 163296, are the norm rather than an exception, we might expect them to eventually form the population of massive exo-planets that has been observed over the last decades. Under this assumption, we studied the synthesized giant planet population from our simulations in the context of known exoplanets of stars similar in size and luminosity to HD 163296.

From the mass and semi-major axis distribution of the simulated planets we conclude that a high α viscosity parameter on the order of 3.16 × 10⁻⁴ ≲ α ≲ 10⁻³ is most likely to produce a planet population that is consistent with current observations. Simulations with low viscosity, that is α = 10⁻⁴, did not generate a sufficient number of close-in planets to explain the observed abundance of giant planets on short-period orbits. We argue that future direct imaging observations with increased measurement sensitivity could further constrain the value of α, still under the assumption that the employed initial population of wide orbit giants is typical of exoplanet systems. If a larger number of planets with m_pl ≲ 3 M_Jup and a_pl ≳ 10 AU will be uncovered in future observations, their existence would point toward lower migration and accretion rates and an intermediate viscosity parameter, that is α < 10⁻³, whereas the lack thereof would imply even higher viscosities on the order of α ≈ 10⁻³.

Our simulations show an increase of the orbital eccentricities with increasing α but even for the simulations with the highest α values the mean eccentricities of the simulated planets are still below those of the known giant planets. We attribute the moderate eccentricities to the comparably low number of planets and point out the growing evidence for a further planet in the system which would likely change the distribution. Although we do not expect to fully reproduce the eccentricity distribution of the observed exoplanets, since we are only modeling an example system, it is nevertheless intriguing that the eccentricities of the simulated planets are comparatively low. In addition, further studies are needed to determine whether most observed RV and direct imaging planetary systems form from the rings and gaps we observe in protoplanetary disks such as HD 163296.

In accordance with the results by Wagner et al. (2019) and Bitsch et al. (2020), we observe a discrepancy between the hypothetically observable and the actual multiplicities of our simulated systems. This implies that many stars known to have one giant planet might in fact harbor several more that are inaccessible to us at current measurement sensitivities.

(Near) Resonant motion is common among all our simulations and often times acts as a stabilizing factor in closely spaced systems. This result is interesting for two reasons. Firstly it is in contrast to findings by Bitsch et al. (2020) who found no evidence for resonant motion of giant planets. This seems to imply that the occurrence of resonances depends on the formation pathways, including the initial number and spacing of the accreting planetary embryos. Secondly, we suggest that the occurrence rates of resonances among detected exoplanets could be used to constrain the α viscosity parameter of the disks they originated from, under the assumptions that we made for the formation pathway. We base this on the fact that the number of observable resonant pairs in our simulations increases for higher values of the viscosity parameter α.

Additional simulations and observations, especially of disks where several planets are expected, could provide valuable insights to constrain the different planet formation pathways even further. In particular, many future observational campaigns will be aimed at fully characterizing targeted systems of particular interest, rather than performing broad surveys of the sky. In this context, theoretical investigations of individual systems, such as this work on the putative planets around HD 163296, rather than of general representatives of the overall sample, represent an essential tool for making sense of the anticipated observational data and expanding our knowledge of planetary formation.

Acknowledgements

B.B. and G.P. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. We also thank the MPIA summer internship program who helped to finance J.M.-H. stay at the MPIA. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

References

All Tables

All Figures

	Fig. 1 Distribution of gas disk lifetimes tdisk used in the simulations. The 100 samples (blue histogram) drawn randomly from the underlying exponential distribution (red line) are displayed as well as the cumulative frequency (black outline). The distribution is limited to lifetimes between 1 ≤ tdisk ≤ 7.5 Myr additional to the current age of the system of approximately 5 Myr.
In the text

Fig. 2 Evolution of three selected planetary systems I (top), II (center), and III (bottom) over the integration time of 100 Myr. The plots display the development over time of the semi-major axes ai (left), the planetary masses mi (center) and the resonance angles ${\psi }_1^1$, ${\psi }_2^1$ of the inner planet pair (right). The disks in all three systems have α = 10−3. System I and III have gas disk lifetimes of l.56Myr and 1.34Myr, respectively, and remain stable throughout the integration whereas system II has a gas disk lifetime of 6.59 Myr and becomes unstable at 2.36 Myr when the inner two planets collide and merge.

In the text

Fig. 3 Mass-distance plot comparing all simulated planets to the known exoplanet population. Small red dots indicate the initial values of the simulations. The pebble isolation mass is shown with a dotted line for α = 10−3. Blue, orange and green dots denote the planets from simulations with α viscosity parameters 10−4, 3.16 × 10−4, and 10−3, respectively. The x-bars extend from the perihelion to the aphelion distances for a given planet. Gray dots denote the exoplanet population detected via radial velocity and black circles mark those around massive stars with luminosities exceeding 5 L⊙. For reference, the planets of HR 8799 are shown as black dots (Wang et al. 2018). A dashed line divides the parameter space into a region observable by current radial velocity measurements with a sensitivity of 1 m s−1 and observation period of up to 10 yr and a currently not observable region. The histograms depict the normalized cumulative frequencies of the final masses (right) and semi-major axes (top). The color scheme is the same as for the main plot, where the bars show the distribution among all simulated planets and the dashed outlines represent the distribution among hypothetically observable planets. The cumulative frequencies among the exoplanets of massive stars with L > 5 L⊙ are plotted for comparison with a black solid outline.

In the text

Fig. 4 Mass vs. semi-major axis plot comparing all simulated planets to the known exoplanet population. The plot is in all aspects the same as Fig. 3, except that the simulated planets are colored to reflect the gas lifetime tdisk of the disk in which they formed instead of the α viscosity parameter. Simulated planets with tdisk < 2 Myr, tdisk between 2 and 5 Myr, and tdisk ≥ 5 Myr are shown as colored dots in blue, turquoise and pink, respectively. It is apparent that the dependence of the final orbital parameters on the disk lifetime tdisk is small compared to the dependence on the viscosity parameter α (see Fig. 3).

In the text

Fig. 5 Cumulative histogram of the eccentricity distribution among the simulated planet population. In the upper panel, the normalized distributions of simulations with α viscosity values 10−4, 3.16 × 10−4, and 10−3, are depicted in blue, orange and green, respectively. Dashed lines indicate the eccentricities among all simulated planets, whereas the eccentricity distributions of only the hypothetically observable planets are shown as solid lines. The lower panel confronts the eccentricity distribution of planets in stable 3-planet systems (purple) with 2-planet systems (red) in which one planet collided or got ejected. For comparison, the eccentricity distribution of observed exoplanets around stars with luminosities exceeding 5 L⊙ is shown as a black solid line in both panels.

In the text

Fig. 6 Normalized cumulative distribution function of the period ratios among simulated planet pairs. In the upper panel, the period ratio distributions of simulations with α viscosity values 10−4, 3.16 × 10−4, and 10−3, are depicted in blue, orange and green, respectively. Dashed lines indicate the period ratios among all adjacent planet pairs, whereas period ratios of those planet pairs that would hypothetically be observable are shown as solid lines. The lower panel compares the period ratio distribution of planets in 2- (purple) and 3-planet systems (red). We observe a significant pile-up of period ratios near the 2:1 resonance for 3-planet systems and around the 3:1 resonance for 2-planet systems.

In the text

Fig. 7 Final period ratios of all simulated planet pairs plotted against their combined mass. Crosses mark planet pairs in 3-planet systems and circles denote 2-planet systems, where the third planet collided or got ejected. The transparency of the markers indicates whether or not the pair would be observable with current standard radial velocity measurements. The different colors denote the α viscosity value of the protoplanetary disk with blue, orange and green for α = 10−4, 3.16 × 10−4, and 10−3, respectively. The period ratios of known exoplanet pairs around stars with luminosities exceeding 5 L⊙ are marked with black dots. The horn-shaped black solid curves indicate the approximate resonance boundaries $\left(p_{out}/P_{\text{in}}-p/q\right)=0.05{\left(\left(M_1+M_2\right)/M_{\text{Jup}}\right)}^{1/2}$ for a p:q mean motion resonance, in this case for 3:2, 2:1, and 3:1 resonances. It is evident that period ratios close to the 2:1 mean motion resonance occur frequently for all α viscosity values. But since for the lower two α viscosity values the distances to the star are typically larger only few of the resonant pairs would be observable.

In the text