© The Authors 2023

1 Introduction

Recent surveys show that planetary systems that host multiple giants are common (e.g., Lissauer et al. 2012; Fabrycky et al. 2014; Zhu 2022). These planets are believed to have formed in the same protoplanetary disk, where they acquired their massive gaseous atmospheres. Previous hydrodynamical studies investigate either the growth of single planets in the disk (e.g., Ayliffe & Bate 2009; D’Angelo & Bodenheimer 2013; Crida & Bitsch 2017; Schulik et al. 2019; Bergez-Casalou et al. 2020) or the evolution of already formed multiple planets (e.g., Baruteau & Papaloizou 2013; Lega et al. 2013; Pierens et al. 2014; Morbidelli et al. 2018). This dichotomy originates from the difficulty to accurately model each evolution process.

It is vital to understand exactly how gas is accreted onto giant planets. Some studies have described how gas accretion can be approximated in 2D, with or without planet migration (Kley 1999; Crida & Bitsch 2017). Other studies used complex 3D simulations to take the various fluid and thermal processes governing the gas accretion of an embedded planet into account (e.g., Ayliffe & Bate 2009; Machida et al. 2010; D’Angelo & Bodenheimer 2013; Szulágyi & Mordasini 2017; Lambrechts et al. 2019; Schulik et al. 2019). Due to their high computational cost, such simulations are often integrated over short timescales, around 100 planetary orbits, making it impossible to investigate the long-term growth of single accreting planets with the current computing facilities.

Moreover, the gas distribution around embedded giant planets is impacted by gap opening. The formation of these gaps has been observed both in the dust (e.g., ALMA Partnership 2015; Andrews et al. 2018; Benisty et al. 2021) and in the gas (e.g., Huang et al. 2018; Smirnov-Pinchukov et al. 2020). Several studies have investigated the characteristics of the gap (such as its depth and width) as a function of the planet and disk characteristics (Lin & Papaloizou 1986; Crida et al. 2006; Fung et al. 2014; Kanagawa et al. 2015). However, due to computational limitations, these studies neglected gas accretion. Recent works show that gap opening and gas accretion are highly dependent on each other through the viscosity of the disk (Bergez-Casalou et al. 2020; Rosenthal et al. 2020).

As gap opening has a non-negligible impact on the gas structure, the growth of a giant planet must be impacted by the presence of a second accreting planet in the same disk. The goal of this study is to quantify this impact as a function of the disk parameters and planet characteristics (i.e., their radial separation and delay in accretion). We perform 2D isothermal hydrodynamical simulations similar to Bergez-Casalou et al. (2020) to monitor the growth of two planets accreting in the runaway accretion regime from the same disk. Our 2D isothermal setup allows us to integrate the evolution of the planet masses for around 0.5 Myr, which is longer than the majority of the studies that investigate gas accretion onto planets using hydrodynamical simulations.

This paper is structured as follows: The numerical setup is described in Sect. 2. In Sect. 3 we investigate the impact of different disk viscosities on the growth of each planet. A comparison with single accreting planets is presented in Sect. 4. Different planet separations are investigated in Sect. 5. In Sect. 6 we show the mass evolution of planets accreting at different times. Our findings are discussed in Sect. 7, where we elaborate on the impact on the stellar gas accretion and on the dynamical evolution of the planetary systems. A comparison with the exoplanet population is discussed as well before we summarize and conclude in Sect. 8.

2 Numerical setup

In this paper we simulate two accreting planets on fixed circular orbits embedded in their gaseous disk with the hydrodynamical code FARGO-2D1D (Crida et al. 2007). This code is composed of two different kind of grids: a 2D grid where the planets are located and which is surrounded by two 1D grids. The first 1D inner disk ranges from 0.1 to 0.78 AU, the 2D grid ranging from 0.78 to 23.4 AU and the second outer 1D disk spans from 23.4 to 260 AU. With this code, the global viscous evolution of the disk is self-consistently modeled at a reasonable computational cost. As mentioned in Bergez-Casalou et al. (2020), it is important to have a consistent viscous evolution of the full disk in order to accurately describe gap formation and gas accretion.

For computational accuracy, the code uses dimensionless units. We normalized masses with the mass of the central star, M₀ = M_⊙ and lengths with the position of the planet, r₀ = 5.2 AU and we set the gravitational constant as G = 1. The unit of time is therefore based on the orbital period at r₀ with P = 2πt₀ where $t_0={\left(r_0^3/\left(G{M}_{*}\right)\right)}^{-1/2}$.

In order to focus on the impact of gas accretion, dynamical interactions between the planets and planet migration are neglected. This choice is discussed in Sect. 7.3. The planets are fixed on circular orbits, at key positions corresponding to different period ratios. We use the term “period ratio” instead of “resonance”, as they are not dynamically locked in resonance but are forced by the code to stay at their position. Four period ratios are investigated: 2:1, 3:1, 4:1, and 5:1. These positions were chosen such as the planets are far enough from each other to be considered dynamically stable during their growth (Chambers et al. 1996; Raymond et al. 2009). In order to compare to a single-planet case, we also simulate the growth of single planets located at positions corresponding to the investigated period ratios. The planets’ configurations are summarized in Table 1.

Each planet starts with an initial mass of 20 M_⊕, allowing it to directly accrete in the runaway gas accretion regime (Pollack et al. 1996). These initial cores are slowly introduced in the disk with the following mass-taper function: $$m_{\text{taper}}={\sin }^2\left(t/\left(4n_{\text{orb}}\right)\right),$$(1)

making the planet grow from 0 to its initial mass in n_orb = 3 orbits of the inner planet. Except for Sect. 6 where the accreting times are specified, all planets start accreting simultaneously after 100 orbits of the inner planet.

The planets accrete following the accretion routine described in Bergez-Casalou et al. (2020). With this accretion routine, the amount of gas accreted by the planets is dictated by Machida et al. (2010; Ṁ_M) and is limited to what the disk can provide by the approach of Kley (1999; Ṁ_K). It is written as follows: $${\dot{M}}_p=\min \{\begin{array}{l}{\dot{M}}_M=<\text{Σ}{>}_{0.9r_H}{H}^2{\text{Ω}}_K\times \min \left[0.14;0.83{\left(r_H\text{/}H\right)}^{9/2}\right]\hfill \\ {\dot{M}}_K={\displaystyle \int {\displaystyle {\int }_{A_{\text{disk}}}{f}_{\text{red}}\left(d\right)}\text{\hspace{0.17em}\hspace{0.17em}Σ}\left(r,\varphi ,t\right)\pi f_{\text{acc}}\text{\hspace{0.17em}d}r\text{\hspace{0.17em}\hspace{0.17em}d}\varphi }\hfill \end{array}$$(2)

where $<\sum {>}_{0.9r_H}$ is the averaged surface density around the planet up to 0.9 r_H with r_H = r_p(m_p/3M_*)^1/3 being the Hill sphere of the planet; A_disk is the disk area; H is the disk scale height; Ω_K is the Keplerian orbital frequency of the planet; d is the distance from the planet; f_acc is the inverse timescale upon which the accretion rate of Kley (1999) is occurring; and f_red is a smooth reduction function that predicts what fraction of gas must be accreted as a function of the distance from the planet. Values for f_acc and f_red are chosen in order to reproduce a realistic accretion efficiency¹. The detailed accretion routine is presented in Appendix A.

Once the amount of gas is determined by the accretion routine, it is removed from the disk and added to the planet’s mass as in Kley (1999). Unless specified, we remain in the regime where Ṁ_M < Ṁ_K throughout, meaning that we always remove the amount of mass suggested in Machida et al. (2010).

The resolution of the 2D grid is such that there are five cells per Hill radius of the inner planet before it starts growing, which leads to N_r = 802 and N_ϕ = 1158. Considering that the resolution is fixed in time, the Hill sphere region will become better and better resolved as the planets grow $\left(r_{\text{H}}\propto m_{\text{p}}^{1/3}\right)$. As we enhance the planet separation, we need to adjust the 2D–1D boundary located at the outer edge of the 2D grid in order to properly take into account the perturbations of the furthest planet. Therefore, in the 4:1 and 5:1 case, this boundary is moved to 36.4 AU, enhancing the radial number of cells of the 2D grid to N_r = 1262. The azimuthal number of cells remains unchanged.

In FARGO-2D1D, the disk is locally isothermal. As in Bergez-Casalou et al. (2020), the surface density profile is chosen such that the total mass of the disk is M_d = 0.1M_*, leading to ∑₀(r₀) = 93.6 g cm⁻² with r₀ = 5.2 AU. Even if this can be considered as a heavy disk (Baillié et al. 2019), its large radial extent allows us to neglect self-gravity. The aspect ratio, h = H/r, of the disk is constant. The disk is subject to an α-viscosity as described by Shakura & Sunyaev (1973). In the following section, different gas kinematic viscosities v = αh²r²Ω_K are investigated by varying the aspect ratio as well as the α-viscosity parameter. The investigated values are h = 0.03, 0.05, 0.07 and α = 10⁻⁴,10⁻³,10⁻², with our fiducial values marked in bold.

3 Influence of the disk viscosity

The flow of gas in the disk is dictated by the kinematic viscosity, v, which depends on the α-viscosity parameter on one hand and on the aspect ratio h of the disk on the other. In this section we investigate the influence of the disk viscosity on the accretion behavior of two planets fixed in a 3:1 period ratio. We start by changing the α-viscosity parameter in Sect. 3.1, then the influence of different aspect ratios is studied in Sect. 3.2. We focus on the influence of the viscosity on the evolution of the planet mass-ratio in a last subsection (Sect. 3.3).

Fig. 1 Planetary accretion rates (top row) and masses (bottom row) as a function of time for different α-viscosities. Here the planets are fixed in a 3:1 period ratio: the inner planet is shown in red and the outer one in blue. As in Bergez-Casalou et al. (2020), the accretion rates slightly increase with increasing viscosity. The oscillations in the accretion rate at low viscosities are due to the presence of vortices. The inner and outer planets display similar accretion rates, leading to similar planet masses.

3.1 Influence of the turbulent viscosity

As mentioned previously, the turbulent viscosity is parametrized by the α-viscosity parameter. Disk turbulence increases with increasing α, leading to faster evolving disks. We show in Fig. 1 the planetary accretion rates (top row) and the resulting masses (bottom row) for planets in disks with different α parameters. From left to right, α increases from 10⁻⁴ to 10⁻². The behavior of the accretion rates is the same as in Bergez-Casalou et al. (2020): reducing the viscosity induces a slightly lower planetary accretion due to a slower flow of gas in the vicinity of the planet.

The difference in accretion rates between the inner (red) and the outer (blue) planet slightly evolves as a function of time and viscosity as shown on the top row of Fig. 1. At the beginning, the planets accrete in the first Machida et al. (2010) accretion regime (dominated by Bondi accretion), leading to a larger inner planet. The flip in accretion rate is due to the switch of accretion regime, from a Bondi to a Hill dominated accretion scheme (Machida et al. 2010; Bergez-Casalou et al. 2020). While the inner planet accretes slightly more in the first regime of accretion, the accretion rates become more similar when the planets accrete in the second accretion regime.

This similarity in accretion rates results in planets of similar masses (bottom row of Fig. 1). Even if the planets are more massive in the high viscosity case than in low viscosity disks, the differences between the inner and outer planet masses does not seem to be highly influenced by the change of α-viscosity. This can be expected as the Machida accretion recipe does not directly depend on this parameter (see Eq. (2)). The influence of the α-viscosity is indirect, as it influences the gas flow around the planet, changing the surface density $<\sum {>}_{0.9r_H}$ from which the planets accrete their gas.

At low viscosity, instabilities are triggered (Klahr & Bodenheimer 2003; Fu et al. 2014). The Rossby wave instability (RWI; Lovelace et al. 1999; Li et al. 2001) produces vortices at the locations of steep pressure gradients, such as produced by the gaps of massive planets. The presence of vortices modifies the flow of gas in the vicinity of the planets, creating the oscillations observed on the top left panel of Fig. 1. In the case of two accreting planets, vortices are produced at three different locations: at the outer edge of the outer gap, in between the planets and even at the inner edge of the inner gap. We show in Appendix B the 2D (r, ϕ) surface density maps of the disk containing the planets in the 3:1 period ratio at low viscosity (α = 10⁻⁴ and h = 0.05) at three different times. The strength of these vortices depend on the growing timescale of the planets (Hammer et al. 2017, 2021). However, here, due to the presence of the second planet, the vortex between the planets quickly vanishes (in less than 4 × 10⁴ yr). The strongest vortex (i.e., with the highest over-density) is the one located at the outer edge of the gap of the outer planet and takes about 10⁵ yr to vanish.

3.2 Influence of the aspect ratio

Another way to study the impact of the disk viscosity is to modify the disk’s aspect ratio. We show the impact of this parameter in Fig. 2. Here, α is fixed to 10⁻³ and h varies from 0.03 (left panel) to 0.07 (right panel). The fiducial value of 0.05 is shown in the middle panel for reference. It should be noted that changing the constant aspect ratio from 0.03 to 0.05 accounts for a reduction of viscosity v by a factor of 2.8 while enhancing the aspect ratio from 0.05 to 0.07 increases the disk viscosity by a factor of 2.0.

Even if the change in the disk kinematic viscosity is weaker than when the α-viscosity is varied in the previous section (change of a factor of 10), the accretion rates behavior are significantly different depending on the aspect ratio. This arises from the direct dependence of the accretion rate recipe on the aspect ratio, while it is indirectly dependent on α. When the gas disk scale height H increases, it requires a larger Hill sphere radius for a planet to switch from the Bondi regulated accretion regime to the Hill regulated one (Machida et al. 2010; Bergez-Casalou et al. 2020). This can be seen on the top row of Fig. 2: the flip in accretion rates occurs at later times for higher aspect ratios. At low aspect ratio (left panel), the planets initial mass is large enough to start accreting immediately in the second regime. In this case, the accretion rates always decrease in time, proportionally to the local surface density.

The timing of the accretion switch between the Bondi and Hill regimes has an important influence on the evolution of the planet masses and their mass ratios. As the inner planet accretes significantly more than the outer one when they are limited to the Bondi accretion regime, a later switch results in more diverse planetary masses. Therefore, the difference between the inner and outer planets increases with an increasing aspect ratio.

Fig. 2 Same as Fig. 1, but for different aspect ratios, increasing from left to right, with α = 10−3. We note on the top of each row the difference in kinematic viscosity caused by the change in the disk’s aspect ratio (v = αh2r2ΩK). The differences in the planetary accretion rates originate from the structure of the accretion formula itself (see Eq. (2)). At a low aspect ratio (h = 0.03, left panel), the accretion starts immediately in the second Machida regime. On the other hand, for h = 0.07, the change of accretion regime occurs at a higher mass and therefore at a later time, leading to a very different planet mass evolution.

3.3 Impact on the planets’ mass ratio

In order to better understand the differences in planetary mass, we analyze in Figs. 3 and 4 the ratios of the planetary masses. We focus first in Fig. 3 on the impact of the œ-parameter on these ratios. The top row represents the mass ratio as a function of time. We arbitrarily decided to show the ratio of the outer planet mass divided by the inner planet mass: this means that when the ratio is decreasing, the inner planet accretes more than the outer one and vice versa. Depending on the viscosity, we see that the ratio shows different flips in time. A first flip ① occurs at each viscosity around 2.5 × 10³ yr. This flip originates from the accretion rate switch discussed in the previous paragraphs: the planets become massive enough to change their accretion regime, resulting in a higher accretion rate for the outer planet (increasing mass ratio).

A second flip ② is observed around 10⁴ yr, at all viscosity except for α = 10⁻² (right panel, which is discussed below). This flip is linked to the formation of deep planetary gaps. In the bottom panels of Fig. 3, we show the perturbed surface densities at the different times marked by the vertical lines in the top panel The perturbed surtace density is defined as the surface density of the disk in the presence of the planets normalized to the surface density of a disk without planets at the same time (Σ_perturbed(t) = ∑_planet(t)/∑_disk(t). These perturbed surface density profiles are used to determine when a gap is opened. We use the definition suggested by Crida et al. (2006): a gap is considered opened when the gas surface density is depleted by 90% compared to a disk without planets. This thresh old is represented by the horizontal gray dotted fine at ∑_planet/∑_disk = 0.1. The second flip ② in the mass ratio occurs when the inner planet reaches this threshold: when the inner planet opens a deep gap, it starts to accrete more than the outer planet.

This link between the gap opening of the inner planet and the decrease of the mass ratio can be explained by the impact of gas accretion on gap opening described in Bergez-Casalou et al. (2020). At low viscosity, gas accretion does not help carve a gap because the disk reaction time is long. The gas is therefore dominantly pushed away from the planets orbit by gap opening, enhancing the surface density in between the planets as well as in the inner and outer regions of the disk. This results in a perturbed surface density larger than one (light gray line in the bottom panels of Fig. 3). When the gap is opened, the planetary gas accretion becomes the dominant process influencing the gas distribution. As the inner planet opens its gap first, it starts by depleting the material present in between the planets and in the inner region of the disk, leading to an inner planet with a higher accretion rate than the outer one. Then, the outer planet also opens a deep gap, helping the inner planet deplete the material located in between them. At this stage, the amount of gas present around the inner planet (at r < r_p,out) is dictated by three different processes. Gas is removed from this region of the disk by (i) the accretion onto the planets and (ii) the accretion onto the star, and is replenished by (iii) the viscous diffusion of the gas from the outer part of the disk through both planet gaps. At low viscosity, only a small amount of gas manages to diffuse through the gaps of both planets. This results in a depletion in gas of the inner disk, resulting in the starvation of the inner planet.

The third flip ③ in the mass ratio, shown in the top panels of Fig. 3, corsesponds to the moment when the outer planet starts accreting more than the starved inner planet, because it is supplied in gas by the outer region of the disk. The depletion of the inner disk can be seen in the corresponding perturbed surface density profiles of the bottom row of Fig. 3: the gray line shows a depleted region in between the planets and in the inner disk. We additionally show the sur face density profiles at 20 000 inner planetary orbits (≃2.4 × 10⁵ yr) after the last mass ratio flip ③: after this time, the inner disk is almost completety emptied in gas, with a perturbed surface density of less than 0.1 within 12 AU. It is clear that after this time, the outer planet will accrete more than the inner one, until it becomes the most massive planet of the system. The only differences between the α = 10⁻⁴ (left panel) and α = 10⁻³ (middle panel) cases are the delay in time of the flips due to different viscous timescales and the presence of vortices at α = 10⁻⁴, influencing the accretion rates or both planets, as mentioned earlier.

The behavior of the mass ratio in a high viscosity disk (right panel) is slightly different than at lower viscosities. As shown in Bergez-Casalou et al. (2020), at high viscosity, gas accretion helps gap formation. Therefore, the material located in between the planets and the inner disk is immediately depleted by gas accretion and viscous stellar accretion. As these two regions are not enhanced in gas by gap opening, the inner planet accretion rate slowly reduces as the inner disk is immediately depleted in gas, meaning that no additional mass ratio flip is observed except for the very first one originating from the accretion recipe. At this viscosity, the gas manages to diffuse efficiently through the gaps of both planets, avoiding a complete depletion of the inner region: this can be seen by comparing the surface density profiles at two different times toward the end of the simulation (marked by the gray and black lines on the right panels). Indeed, the difference in the profiles after 2000 orbits of the inner planet is small, meaning that the disk gas flow is high enough to prevent the total depletion of the inner region, unlike at lower viscosities (middle and left panels). However, in this configuration, even if the material around the inner planet is replenished by viscous evolution, the amount of gas diffused through both gaps is not high enough to allow the inner planet to accrete more than the outer planet.

These behaviors at high and low viscosity are also observed when we vary the disk aspect ratio. In Fig. 4, we show the evolution of the mass ratios as a function of time (top row) as well as the perturbed surface density profiles at given times (bottom row) for the different aspect ratios presented in Fig. 2. As the behaviors described in the previous paragraphs are only dependent on the gas kinematic viscosity v, we recover the same behaviors when we change the aspect ratio: at low viscosity, the mass ratios are highly influenced by gap formation whereas at high viscosity, the mass ratio does not present more than one flip. However, due to the dependence of the accretion recipe on the aspect ratio, the amplitudes of the mass ratios are highly dependent on h. While the mass ratios reached values between 0.8 and 1.1 for the different α parameters (Fig. 3), here the planets show a larger spread in mass ratio at the beginning of the simulations in disks with larger aspect ratios. However, once the gaps are formed in all cases and the inner disk is depleted, all the mass ratios stay quite close to 1, with an outer planet less than 1.2 times more massive than the inner planet.

As a conclusion, the mass ratio between the planets is highly and mainly dependent on the disk kinematic viscosity v. The resulting mass ratios are close to 1 even after 0.5 Myr of evolution, leading to planets with rather similar masses. In all the explored cases, we expect the outer planet to become the most massive planet of the system.

Fig. 3 Mass ratio (top row) and perturbed surface density proflies (bottom row) as a function of time for different α-viscosities. As in Fig. 1, the planets are fixed at the 3:1 period ratio positions, represented by the two vertical dotted lines in the bottom row. The horizontal dotted lines in the bottom panels mark the ∑planet/∑disk = 0.1 threshold, where we consider that a gap is opened (Crida et al. 2006). In the top panels, a decreasing (respectively increasing) mass ratio indicates that the inner (respectively outer) planet is accreting more than the other planet. Different flips are observed in the evolution of the mass ratios, marked in circled numbers. The color of the surface density profiles shown in the bottom row corresponds to the color of the different vertical lines in the top row and represents different times. The last snapshot of the perturbed surface density is taken 2.4 × 105 yr (20 000 orbits of the inner planet) after the last mass ratio flip, shown by the vertical dotted black line. The first flip ① can be explained by the accretion formula itself (see Eq. (2)), while the other two flips correspond to the evolution of the surface density.

Fig. 4 Same as Fig. 3 but for different aspect ratios, with α = 10−3. As in Fig. 2, we note at the top of each row the difference in kinematic viscosity caused by the change in the disk’s aspect ratio. Here the mass ratio evolutions differ in amplitude, but they present the same behavior as in Fig. 3. The difference in amplitude originates from the direct dependence of the accretion routine on the aspect ratio, and not on the α-viscosity (see Eq. (2)).

4 Single accreting planet compared to two accreting planets

The presence of the second planet highly influences the growth of the inner planet. In order to quantify the effect of the presence of the second planet, we compare the growth of the two planets to two separate simulations where the planets are alone in their disk.

4.1 Accretion rate and mass comparisons

It has been shown in Sect. 3.3 that the outer planet has the capacity to starve the inner planet once the gaps are formed and the inner disk is depleted. Consequently, the differences between two accreting planets and single accreting planets should increase with time. We show in Fig. 5 the comparison between three different simulations: the first one, represented by the blue and red colors, considers the simultaneous accretion of two planets in the 3:1 period ratio as in the previous section; in the second simulation, the disk hosts only a single planet located at the position of the inner planet (purple line); in the third simulation, the disk hosts also a single planet located this time at the position of the outer planet (cyan line). Each planet configuration is represented in the top left corner of the figure. As in Sect. 3.1, different α are shown in the different columns, increasing from left to right.

The different planetary gas accretion rates are shown in the top row of Fig. 5. At high viscosity, in the right panel, the accretion rate between the single planets and the two planets are very similar. As expected, only the accretion rate of the inner planet is significantly impacted by the presence of the second planet: after 2 × 10⁴ yr, the accretion rate of the inner planet in the two-planet case (red line) starts to be reduced compared to the single-planet case (purple line). This is particularly visible in the bottom right panel of Fig. 5, where we compare the masses in the two-planet case to the single-planet case: the red line represents the mass ratio of the inner planets (m_2p,in/m_1p,in) while the blue line represents the mass ratio of the outer planets (m_2p,out/m_1p,out). After 2 × 10⁴ yr, the red line continuously decreases, meaning that the planet at the inner position has a reduced accretion rate in the two-planet case compared to the single planet. This is due to its starvation, as discussed in the previous section.

At high viscosity, the outer planet is only slightly impacted by the presence of the inner planet. While their accretion rate seems very similar (blue and cyan line in the top right panel), the mass ratio shows a slight reduction of the mass in the two-planet case (blue line slightly below one in the bottom right panel). This originates also from the depletion of the inner disk: the gas accretion rate of the outer planet only relies on the material located outside the planet’s orbit in the two-planet case, while the single planet accretes material from both the outer and inner disk.

At lower viscosities, the impact of the presence of the second planet occurs earlier and is more significant. Focusing on the intermediate viscosity (α = 10⁻³, in the middle column), we see that the accretion rates (top panel) in the two-planet case are enhanced after ~10⁴ yr for both the inner and the outer planet compared to the single planets from this point in time. This enhancement, absent at high viscosity, originates from the gap opening process: as mentioned in Sect. 3.3, at this viscosity, gap opening pushes material away from the vicinity of each planet, enhancing the surface density around them. In this case, the planets are close enough to each other to push material in the feeding zone of the neighboring planet. Therefore, the inner planet pushes material toward the outer planet and vice versa. Each accretion rate is enhanced until the depletion of the inner disk (r < r_p,out), resulting in lower accretion rates than in the single-planet cases. This enhancement can also be seen in the evolution of the planet masses (bottom and middle panels of Fig. 5 for α = 10⁻³). In the long term, the evolution of the mass ratios between the two-planet and single-planet cases are expected to behave like at high viscosity: the inner planet will be starved in gas, leading to a smaller inner planet in the two-planet case compared to the single-planet case. As for the outer planet, once the inner disk is depleted, it is fed only by the outer disk in the two-planet case whereas the single outer planet accretes also from the inner disk. This leads to a decreasing outer planet mass ratio m_2p,out/m_1p,out until reaching an almost constant ratio, similarly to the high viscosity case.

Regarding the lowest viscosity case (α = 10⁻⁴, in the left column of Fig. 5), the presence of additional vortices located in between the planets significantly alters the accretion rates of the two accreting planets compared to the single planets. Even if the single planets also produce vortices, leading to the oscillations observed in their accretion rate too, the additional vortex present in between the planets quickly enhances the accretion of both planets. Except from this non-negligible influence, the overall evolution of the planet masses follows the trend observed at α = 10⁻³, where we expect the different flips to occur at later times due to the larger viscous timescale.

Finally, we also compare the masses of the outer planets with the inner planets, as in Sect. 3.3. At intermediate viscosity (α = 10⁻³, middle panel of Fig. 6), the mass ratio in the single-planet cases (dashed line) features only one flip ① before reaching a plateau after ~2 × 10⁴ yr. As before, the first flip ① is due to the accretion recipe, and it occurs at the same time in the two-planet case as in the single-planet cases. The second flip ② in the case of the two planets in the same disk, due to gap opening, is absent in the single-planet cases. This is expected because it is the presence of the second planet creating a gap that enhances the accretion of the inner planet, which is not the case in the single-planet simulations. With our disk profile (a surface density power law of −1 and constant aspect ratio), the single planet located at the outer position accretes naturally more than the inner planet, leading to an increasing mass ratio. The plateau observed after ~2 × 10⁴ yr corresponds to the moment when the single planets open a deep gap: then, the accretion of each planet is mainly governed by the flow of gas originating from the outer region of the disk. With our disk profile and due to the close proximity of the planets, the gas flow from outside of the inner planet position is similar to the flow from outside of the outer planet position, leading to similar accretion rates. The same behavior is observed at low viscosity (α = 10⁻⁴, in the left panel of Fig. 6), considering the perturbations and enhancement produced by the vortices.

At high viscosity (α = 10⁻², right panel of Fig. 6), the behavior of the single-planet simulations is again different compared to the simulations with two planets. The absence of additional flips after 2 × 10³ yr in the two-planet case is explained by the relative absence of impact of gap opening at this high viscosity (see Sect. 3.3): the viscosity is high enough to cause the immediate depletion of the inner disk via both viscous accretion toward the star and the accretion of the inner planet, immediately leading to its starvation. However, in the mass ratio of the single planets, we observe another flip ② around 4 × 10⁴ yr. After this time, the inner planet accretes more than the outer one. This originates from the flow of gas reaching the inner disk once each gap is opened: at this viscosity, the gas can significantly flow through the gaps of the planets. However, it is easier to flow through the gap of the single planet located in the inner region as it is less wide than the gap of the single planet located in the outer region. This results in a inner disk that is more depleted in the case of the outer single planet than in the case of the inner single planet, leading to the reduction of the accretion rate onto the outer single planet.

Fig. 5 Comparison between single accreting planets and two simultaneously accreting planets. The planets are fixed at positions corresponding to the 3:1 period ratio. As in Figs. 1 and 2, the top row presents the accretion rates as a function of time and the middle row the planet masses. In the bottom row, we show the ratio of the masses in the single- and two-planet case: the red line represents the ratio of the inner planets (m2p,in/m1p,in), and the blue line represents the outer planets’ ratio (m2p,out/m1p,out). At low viscosity (left panels), the differences with the single-planet case originate from the additional formation of vortices in between the planets, which enhances the accretion rates during the vortex lifetimes. At high viscosity (right panel), while the outer planet is slightly impacted by the presence of the inner planet, we observe the starvation of the inner planet. We expect to see similar behavior at lower viscosities but delayed in time, due to longer viscosity timescales.

Fig. 6 Comparison of the mass ratios in the single-planet and two-planet cases as a function of time and for different viscosities. Again, the planets are fixed at the 3:1 period ratio positions, as shown in the upper-left corner. In solid lines we show the mass ratios of the outer over the inner planet in the two-planet case, whereas the dashed lines represent the mass ratio of the single planets. Different flips are marked by circled numbers and can be explained by the impact that gap opening and depletion of the inner disk have on the accretion behavior of the planets.

4.2 Impact on the accretion onto the star

The presence of a single gap opening planet can alter the gas accretion onto the central star (e.g., Manara et al. 2019; Bergez-Casalou et al. 2020). Here, we investigate the impact of the presence of a second planet on the evolution of the stellar gas accretion. The stellar gas accretion rate is defined by the flow of mass through the inner edge of the disk: Ṁ_* = −2π r_inυ_r,in∑_in, where υ_r,in and ∑_in are the radial velocity and surface density at the inner boundary located at r_in = 0.2 AU. In Fig. 7, we compare the stellar accretion rates of disks that contain two accreting planets in the 3:1 period ratio (solid black line), disks that contain single planets located at the inner (purple dashed line) and outer (cyan dotted dashed line) positions of the two-planet simulation, or disks that host no planets at all (gray dotted dashed line). Again, the α-viscosity increases from left to right. Independent of the viscosity, the stellar accretion rate in the presence of two planets features oscillations: they originate from the periodic overlap of the planets spiral density arms, locally enhancing the surface density of the gas.

At intermediate viscosity (α = 10⁻³, middle panel of Fig. 7), the presence of the planets has two distinct impacts. Before 10⁴ yr, the planets are not large enough to influence the accretion onto the star. Between 10⁴ yr and 10⁵ yr, the disks hosting planets harbor an enhanced stellar accretion rate compared to the disk without planets. This originates from gap opening: material is pushed to the inner region, feeding the star. After 10⁵ yr, planetary gas accretion and gaps prevent part of the gas from reaching the inner region, leading to the decrease in the stellar accretion rate. It takes more time for the disk hosting the single planet located at the outer position to reduce the stellar accretion rate due to the larger inner disk present in this configuration. After a given time (here after 2 × 10⁵ yr), the flow of gas through the gaps reach a quasi equilibrium state, leading to a linear evolution of the stellar accretion rate, with a similar slope compared to the slope of the stellar accretion rate of the disk without planets. At this stage, the stellar accretion rate is reduced by a factor of between 4 and 5 when the disk hosts two accreting planets compared to the disk without planets.

The enhancement produced by the two accreting planets is slightly more pronounced than in the single-planet case as it includes the material pushed by both planets. As one could expect, the decrease of the stellar accretion rate in the two-planet case occurs at the same time as in the disk hosting the single planet located at the inner position. However, at this viscosity, the reduction of the stellar accretion rate is only barely influenced by the presence of the second planet, resulting in a decrease of only ~30% compared to the single inner planet case. This means that the flow of gas reaching the inner disk is mostly governed by the influence of the inner planet.

At high viscosity, the viscous spreading of the gas prevents the enhancement of the stellar accretion rate: the material pushed to the inner regions by the gap opening planets is both less important than at lower viscosity and quickly diffused toward the star by viscous spreading. This results in a reduction of the stellar accretion rate after 10⁴ yr only. The gas is efficiently diffused through the planet gaps, maintaining a high stellar accretion rate, even in the presence of multiple accreting planets. Here, the stellar accretion rate is only reduced by up to a factor of 2.5 compared to the accretion in a disk without planets. Moreover, the presence of the second planet only influences the stellar accretion rate by less than 20% compared to the cases with single planets. Focusing on the single-planet cases, we see that the stellar accretion reaches a similar equilibrium in both cases, independent of the planet location. As expected, with this high viscosity, the presence of accreting planets barely influences the gas disk evolution.

The opposite behavior is observed at low viscosity (α = 10⁻⁴, left panel of Fig. 7). Here, the stellar accretion rate is highly influenced by the vortices and by the gaps formed by the planets. Material is efficiently pushed in the inner regions, enhancing the stellar accretion rate by up to a factor of 10 in the disk hosting two planets compared to the accretion of a disk without planets. The reduction of the accretion occurs at later time than at higher viscosity, which is expected because the planets grow slower and gap opening takes more time (Bergez-Casalou et al. 2020). However, the presence of the second planet in the disk highly influences the stellar accretion rate, because the viscosity is not high enough to diffuse gas through both gaps. Due to the long viscous timescale, we are not able to determine the final reduction factor compared to the higher viscosity cases. At the end of the simulation containing the two planets, the stellar accretion rate is reduced by a factor of 3 compared to the disk without planets, giving a lower estimation of the reduction factor. In Sect. 7.1, we discuss the impact that this reduction factor can have on observations and compare it to other studies.

Fig. 7 Influence of a second accreting planet on the stellar accretion rate at the inner edge of the disk (0.2 AU). As in Figs. 5 and 6, the planets are located in the 3:1 period ratio positions, with an increasing α-viscosity from left to right. The comparison is made between the two-planet case (solid black line), the single-planet cases (dashed purple line for the single inner planet and dotted cyan dashed line for the single outer planet) and a disk without planets (dotted dashed gray line). The oscillations present in the case of the two accreting planets are due to the overlap of their spiral arms, coupled with the presence of vortices at low viscosity. This sometimes results in a negative radial velocity, explaining the dip seen at α = 10−4 due to the logarithmic scale. Compared to a disk with no planet, the accretion onto the star is only reduced by up to a factor of 3 in the presence of multiple giant planets, similar to the scenario with a single planet (Bergez-Casalou et al. 2020).

5 Influence of the planet separation

Our study has focused so far on planets placed at positions corresponding to the 3:1 period ratio. However, the separation between forming planets is not fixed in time as they can dynamically interact with the disk and with each other (e.g., Baruteau et al. 2014; Crida & Bitsch 2017). While we will investigate the impact of the radial evolution of the planets on their growth in a future study, we study in the following section the impact of the planet separation in their growth by placing the planets at different period ratios.

Fig. 8 Mass ratio as a function of time for different α-viscosities and different period ratios. The darker the line, the farther away the planets are located from each other. The separation between the planets has a small impact on the range of their mass ratio, meaning that planets that start accreting at the same time will have similar masses. The behavior of the mass ratio evolution is not impacted by the planet separation: at all viscosities the mass ratio flips are due to the same reasons described in Fig. 3. The exception is for the 2:1 period ratio at high viscosity (right panel): here the planets are close enough and the viscosity is sufficiently high to maintain a significant flow toward the inner disk through the gaps, preventing the inner planet from being starved. This results in an additional mass ratio flip, ⦷, similar to the high viscosity single-planet case of Fig. 6.

5.1 Impact on the mass ratio

As shown in Sect. 3.3, the planet mass ratio evolution reflects the accretion history of the planets. We show in Fig. 8 the mass ratio evolutions of two simultaneously accreting planets located at different period ratios, ranging from 2:1 up to 5:1. These period ratios, represented to scale in the top left corner of the figure, were chosen such that the planets can be considered dynamically stable during their growth, as we neglect their dynamical interactions. Their stability is monitored thanks to their mutual Hill radii (Chambers et al. 1996). As in Fig. 3, the mass ratio presented is the ratio of the outer planet mass divided by the inner planet mass and each panel represents an α-viscosity, increasing from left to right. A darker line depicts a larger planet separation.

Independent of the viscosity, we first observe that the first mass ratio flip ①, originating from the switch of accretion regime in our accretion recipe, is delayed in time when the planets are farther away from each other. Indeed, with our disk profile, a planet located farther away in the disk has a slightly lower gas accretion rate, meaning that it will need more time to reach the mass needed to switch from the Bondi to the Hill accretion regime. This delay affects the similarities between the planets: the farther the planets are from each other, the more extreme are their mass ratios at the beginning of the simulations (i.e., the deeper the first mass ratio flip ① is).

Focusing onthe high viscosity behavior (α = 10⁻², right panel of Fig. 8)), we observe the same behavior as described in Sect 3.3 except for the 2:1 period ratio. While the other period ratios feature only the first flip ① due to the rapid depletion of the inner disk via viscous stellar accretion, a second flip ⦷ is observed at around 10⁵ yr in the mass ratio of the 2:1 period ratio case. In this configuration, the planets are close enough from each other to facilitate the diffusion of gas through both gaps compared to the other period ratios, as they quickly form a common gap. As the flow of gas to the inner disk is higher, more material is present around the inner planet. In this case, the amount of gas diffusing through both gaps, is high enough to prevent the starvation of the inner planet, leading to a decreasing mass ratio.

The behavior of the mass ratio at α = 10⁻³ (middle panel of Fig. 8) is also the same for each period ratio. The second mass ratio flip ②, occurring at around 10⁴ yr is always due to the formation of the planetary gaps. The delay in the flip originates from the time needed for she outer planet to also significantly enhance the surface density in between the planets via gap opennig. The amplitude of the mast ratio between the second ② and third flip ③ (i.e., the difference between the local maximum ② and the local minimum ③) is decreasing with decreasing planet separation. As this moment corresponds to the moment when the inner planet is accreting more thanks to the enhancement of material in between the planets and in the inner disk after gaps are formed, closer planets have less materias in between them by construction leading to a smaller mass ratio amplitude. We note that it also takes more time for planets located farther away from each other to empty the material located between the planets, meaning that the inner planet accretes more than the outer planet for a longer time compared to planets closer to each other. Therefore, the inner planet is starved more easily the closer the planets are.

At low viscosity (α = 10⁻⁴, left panel of Fig. 8)), the planet separation has an additional impact on the formation of vortices. In the 2:1 period ratio simulation, the planets quickly create a common gap, preventing the formation of a strong vortex in between them. As we mentioned in Sect. 3.1, as vortices push material toward the planets, it means that the gas accretion of the planets in the 2:1 period ratio are less impacted by the presence of vortices compared to planets located farther away from each other. In the 5:1 period ratio case, the planet gaps are clearly distinct from each other. Each planet therefore create sharp density gradients at the inner and outer edge of their gaps, resulting in the formation of four vortices. Even though the vortices located in the inner disk (i.e., inner to the outer planet) dissipate quickly, they influence the evolution of the planet masses. Except from the formation of vortices, the global behavior of the mass ratio is due to the same processes as at intermediate viscosity.

Overall, independent of the planet separation or disk viscosity, the mass ratios stay close to one (0.7 < m_out/m_in < 1.25). We discuss in Sect. 7.5 how does this compare to the observed planetary systems.

Fig. 9 Gap-opening mass as a function of the viscosity for different criteria and our simulations, as in Bergez-Casalou et al. (2020). In the left panel, we show the gap opening masses of the inner planets in the two-planet case compared to the single planet cases located at the inner position. The outer planet cases are shown in the right panel. In each panel, the lines represent the different gap-opening criteria from the literature: Crida et al. (2006) in dashed blue, Fung et al. (2014) in solid orange, and Gyeol Yun et al. (2019) in dashed-dotted red. As shown with the schematics above the panels, the two-planet cases are represented by triangles (upward for the inner planets and downward for the outer planets), and the single planets are represented by circles. The inner single-planet case is shown in gray, and the colored circles correspond to the colors of the outer positions in each configuration. For clarity, in the right panel, the viscosity is divided into discrete intervals, allowing us to compare the gap opening mass at one given viscosity in the different configurations. The gap opening masses of the inner planet are barely impacted by the presence of a second outer planet. On the other hand, at low viscosity, the closer the planets are to each other, the higher the gap opening mass is compared to the single-planet case.

5.2 Impact on the gap opening mass

The gap opening mass is an important parameter used both in theoretical models to approximate when a planet switches from the fast type I migration to the slow type II migration (e.g., Ndugu et al. 2018, 2021; Bitsch et al. 2019; Miguel et al. 2020) and in dust observations to indirectly derive the masses of embedded planets from the observed characteristics of gaps (e.g., Zhang et al. 2018; Asensio-Torres et al. 2021). As we showed in Bergez-Casalou et al. (2020), gas accretion has a non-negligible impact on the gap opening mass. In this section we investigate the influence of the presence of a second accreting planet on the gap opening mass of each planet.

We compare in Fig. 9 the gap opening mass of our accreting planets to different gap opening criteria derived in previous studies (Crida et al. 2006; Fung et al. 2014; Gyeol Yun et al. 2019). The gap opening masses of the inner planets are shown on the left panel of the figure while the outer planets are shown on the right panel². As in the previous section, the color represents the different period ratios, as can be seen on the schematics of the disk configurations represented in the top of the figure. Triangles represent the gap opening masses of the planets located in a two-planet disk and single-planet gap opening masses are shown with circles. The gray circles correspond to the gap opening masses of the single planets located at the inner location (they correspond to the gap opening masses presented in Bergez-Casalou et al. (2020) for the fiducial gas accretion rate). The different lines represent the different gap opening criteria compared in this study. In the right panel, the viscosity is divided into discrete intervals for each of the studied viscosity (α = 10⁻⁴, 10⁻³, 10⁻²) to help in visualizing the different configurations.

In Bergez-Casalou et al. (2020), we conclude that gas accretion has a different impact on the gap opening mass depending on the disk viscosity: at high viscosity, gas accretion helps carve deeper gaps, resulting in a lower gap opening mass for an accreting planet while at low viscosity, gap formation is not helped by gas accretion, resulting in a higher gap opening mass for an accreting planet. When a second planet is added in the disk, the accretion of each planet is impacted by the gap opening of the neighboring planet, as shown in Sect. 4. However, this impact depends on the viscosity.

At high viscosity, gas accretion helps carve a deeper gap. Here, with α = 10⁻² and h = 0.05, the disk is at the intersection between the high viscosity regime and the low viscosity regime described above, meaning that the gas accretion has no important impact on the gap opening mass. The gap opening masses are therefore solely dependent on the amount of gas diffusing through the gaps. As the gas diffuses efficiently in this case, gap formation is not impacted by the presence of a second planet. Indeed, the inner disk is more depleted by viscous accretion than by the accretion of both planets, meaning that the material pushed away by the gap forming planets is dissipated via viscous spreading. This can be seen in Fig. 9, where the gap opening masses are the same in the single or two-planet case, for both the inner and outer planets.

At lower viscosities, gas accretion does not help gap formation (Bergez-Casalou et al. 2020). Therefore, the gap opening mass depends on the accretion rate of the planet. As presented in Sect. 4.1, the accretion rates of the two accreting planets are slightly enhanced compared to the single planets due to the formation of the planetary gaps pushing material in the feeding zones of the neighboring planet. This slight enhancement of the accretion rate of the planets leads to slightly higher gap opening masses. This impacts both the inner and the outer planets, as shown in Fig. 9. As the gap opening mass then relies on the amount of material pushed toward the neighboring planet, planets that are close enough from each other enhance their gap opening mass more. Intuitively, when the planets are farther from each other, they tend to behave as if they are isolated and have gap opening masses closer to the single-planet simulations.

Overall, the gap opening masses are barely impacted by the presence of a second planet in the disk. The small differences originate from the differences in accretion rates as the planets push material toward each other. Considering the conclusions of Bergez-Casalou et al. (2020), it seems that the gas accretion rate on the planet itself has a stronger impact on the gap opening mass than the presence of a simultaneously accreting companion.

Fig. 10 Mass ratio as a function of time for different accretion delays on the inner (purple lines) or the outer planet (green lines). The planets are located in the 2:1 period ratio configuration. The two panels represent two different viscosities: α = 10−3 in the left panel and α = 10−2 in the right. Darker lines represent shorter delays. Note that this time the mass ratio is displayed on a logarithmic scale (previous plots are on a linear scale) due to the extreme mass ratios induced by our initial setup here. The dots mark the moment when both planets have reached at least 0.3 MJ and can be considered gas giant planets. The gray rectangle represents the region where the mass ratio lies between 0.25 and 0.5 and shows the region where Jupiter and Saturn meet the required conditions to enter common outward migration, for α ≲ 10−3 (Masset & Snellgrove 2001; Pierens et al. 2014). Independent of the viscosity, all the mass ratios quickly tend toward the mout/min = 1 line, leading to similar planet masses in both cases (0.8 < mout/min < 2 after 105 yr).

6 Influence of delayed accretion

Giant planet formation models have very few constraints on the timing at which runaway gas accretion occurs (e.g., Paardekooper & Johansen 2018; Raymond et al. 2020). So far, we only considered the simultaneous accretion of both giants. However, depending on the disk local properties and on the formation mechanism, giant planets located in the same disk could start accreting at different times. We investigate here the influence of the delayed accretion on the evolution of the planetary growth. Different time delays are considered, on the outer and on the inner planet. We base our delays on the mass of the neighboring planet: the accretion on the second planet is allowed when the other planet reaches 0.3 M_J, 0.5 M_J and 1 M_J. We chose to investigate the impact of the accretion delay on the 2:1 period ratio configuration at high and intermediate viscosities as they reach these masses in a reasonable computational time.

All the resulting mass ratios are shown in Fig. 10. The mass ratio of planets simultaneously accreting is shown and corresponds to the 2:1 mass ratio presented in Fig. 8. We note the difference in the mass ratio scale: while before the mass ratios are shown on a linear scale, here the scale is logarithmic for readability. The large spread in mass ratio is induced by our initial choice for the different delays. Indeed, as we wait for the neighboring planet to reach a given mass before accreting, this sets the maximal and minimal mass ratio reached by the planets: the maximal value that is reached is 16 (318/20 M_⊕) and the minimal one is 0.06 (20/318 M_⊕). We therefore expect the planets to reach a final mass ratio located in between these initial values.

The evolution of the different mass ratios is very different from the simultaneously accreting planets. When the outer planet accretion is delayed (green lines in Fig. 10), it allows the inner planet to accrete slightly more gas before being starved by the growth of the outer planet. At high viscosity, the gas diffuses efficiently through the planet gaps, depleting the whole disk in gas and leading to a high stellar accretion rate (see Sect. 4.2). Therefore, a longer delay results in accretion in a more depleted disk. As a consequence, the mass ratio is lower for longer delays.

At lower viscosity, the effect of viscous spreading as described above can be perturbed by the formation of the inner planet gap. However, the gap opening mass at this kinematic viscosity is around 0.5 M_J at the location of the inner planet. For a delay of 300 orb. (~4.8 × 10³ yr), corresponding to an inner planet mass of 0.3 M_J, the outer planets starts accreting while the inner planet did not create a deep gap yet. This leads to a very similar final mass ratio evolution as in the simultaneous case. However, we know from Sect. 3.3 that when the inner planet creates its gap, it starts to deplete the material located in the inner disk and in between the planets. Therefore, when the outer planet starts accreting even later (e.g., when the inner planet has reached its gap opening mass), the inner planet already depleted part of the material in between them, leading to a lower gas accretion rate of the outer planet compared to the simultaneous case. This results in lower mass ratios for longer delays.

The mass ratio of the planets in the case of the delayed accretion of the inner planet (purple lines of Fig. 10) highly depends on the depth of the gap of the outer planet and on the viscosity of the disk. Indeed, at low viscosity, if the inner planet starts accreting before the formation of the gap of the outer planet, then the inner region is not depleted in gas yet. The behavior of the mass ratio then quickly tends to be the same as in the simultaneous case. However, if the outer planet already opened its gap, then a longer delay of accretion results in a more depleted inner disk. The inner planet has therefore less material to accrete, leading to higher mass ratios. At high viscosity, the same behavior occurs, with the inner disk being efficiently depleted by stellar accretion.

Even if different mechanisms influence the mass ratio in the case of delayed accretion, the giant planets always end up with rather similar masses: 0.8 < m_out/m_in < 2 after 10⁵ yr. We discuss in Sect. 7.5 how does this compare to the observed mass distributions in the exoplanet population and what kind of constrains on planet formation can be derived.

7 Discussion

7.1 Accretion onto the star

In Sect. 4.2 we investigate the influence of the presence of multiple gas accreting planets on the stellar accretion at different viscosities. The results are compared to the stellar accretion in disks hosting single planets and in disks without any planet. We find that the presence of the second planet only has a significant effect when the viscosity of the disk is low (α ≲ 10⁻⁴). For higher viscosities, the presence of the second planet only influences the stellar accretion rate by up to 30% compared to the case with single planets.

We compare our results with two planets to the different stellar accretion rates obtained with different planetary accretion rate in Bergez-Casalou et al. (2020). At high viscosity, the presence of the second planet has less impact on the stellar accretion rate compared to a significant enhancement of the planet accretion rates. Due to the uncertainties in gas accretion rates (e.g., spread in the values found in the following studies: Kley 1999; D’Angelo et al. 2003; Tanigawa & Ikoma 2007; Ayliffe & Bate 2009; Machida et al. 2010; Tanigawa & Tanaka 2016; Crida & Bitsch 2017; Schulik et al. 2019; Lambrechts et al. 2019), it is impossible to use the stellar accretion rates to determine if the protoplanetary disk hosts a single fast accreting planet or multiple planets accreting at a lower rate.

As discussed in Bergez-Casalou et al. (2020), our results are quite different from the study derived by Manara et al. (2019). In their models, they find that the stellar gas accretion rates can be reduced by over two orders of magnitude when the disk is hosting accreting giant planets. We showed in our previous paper that these large spreads of stellar accretion rates could only be reached by widely changing the disk viscosity (over several orders of magnitude). With this study, we additionally show that the presence of a second accreting companion cannot explain such a large reduction in the stellar accretion rate either. Indeed, while reducing the disk viscosity to α = 10⁻⁴ enhanced the impact of the two accreting planets on the stellar accretion rate, we expect the reduction to be of a factor of 10 at most compared to a disk without planets. These discrepancies between our study and the work done by Manara et al. (2019) are the same as mentioned in our previous paper: their model simulates a 1D gas disk while our study is performed in 2D, allowing us to more accurately determine the flow of gas through the gaps of the planets (e.g., Lubow & D’Angelo 2006). Moreover, their planetary gas accretion rates might be overestimated as they rely on the unperturbed surface density (Mordasini et al. 2012), while we showed here that the depletion of the inner disk leads to the starvation of the inner planet and consequently a reduction of its accretion rate.

Again, our simulations indicate that planetary gas accretion might have a smaller impact than expected on the stellar gas accretion rates, even in the presence of multiple accreting planets.

7.2 Dependence on the accretion behavior

The results from our simulations rely on our accretion recipe, based on the results from Machida et al. (2010) and on the amount of gas available around the planets (i.e., the total mass of the disk and the gas surface density profile). However, one can argue that our understanding of planetary gas accretion in the runaway phase can be improved. In the current literature, the derived accretion rates can range over several orders of magnitude, ranging from around 10⁻⁸ M_J yr⁻¹ (Tanigawa & Ikoma 2007; Tanigawa & Tanaka 2016) up to around 10⁻⁴ M_J yr⁻¹ (D’Angelo et al. 2003; Schulik et al. 2019). Furthermore, some studies (e.g., Lambrechts et al. 2019) show that the actual accretion rate has to be separated from the flux that goes through the atmosphere as the planet is embedded in the gas disk. With this definition, the flux through the atmosphere is on the order of 10⁻⁴ M_J yr⁻¹ but the actual accretion rate is lower, around 10⁻⁶ to 10⁻⁵ M_J yr⁻¹.

However, in our study, the global evolution of the planet’s mass ratio relies mostly on the evolution of the global gas distribution in the disk, via gap opening and the depletion of the disk. Therefore, the way the disk is evolving has a larger impact than the actual accretion rate.

As mentioned in Sect. 2, we start our simulations with a relatively massive disk with a surface density profile following Σ ∝ r⁻¹. Implanting the planets early implies that a large amount of gas is still available for the planets to accrete. If we wait for a longer time before introducing the planetary cores, the disk has time to be accreted onto the star and therefore less material would be available around the planets. This will lead to a lower planetary accretion rate on one hand, and on less material between them on the other hand. Then it is just a question of timing again between the formation of the gaps and the depletion of the inner disk and the material between the planets: the global evolution of the planet mass ratio is expected to be the same but on longer timescales as the gap opening mass is independent of the disk’s surface density (Crida et al. 2006; Fung et al. 2014; Kanagawa et al. 2015). Therefore, the timing of formation (and therefore the total mass of the disk) will be of importance regarding the total mass reached by the planets but will not have a huge impact on the evolution of their mass ratio. Regarding the importance of the surface density profile, the evolution of the planets’ mass ratio will also be shifted in time but the qualitative evolution is expected to be similar.

7.3 Impact of planet dynamics

In order to determine the impact of the gas accretion on two planets embedded in the same disk, we neglected both the dynamical interactions between the planets and their migration. Regarding migration, different studies investigate how it impacts gas accretion (e.g., Dürmann & Kley 2015; Crida & Bitsch 2017; Dürmann & Kley 2017). The main results of these studies are that the evolution of the planet characteristics (i.e., mass and semimajor axis) highly depends on the timescales of each process: a fast migrating planet tends to accrete more gas as it quickly moves toward regions with high surface densities. However, as we showed in Bergez-Casalou et al. (2020), gas accretion has an impact also on the gap opening mass, influencing the migration speed of the planet as it transitions from a fast type I to a slow type II migration. This effect was observed in Crida & Bitsch (2017), where their accreting planet slowed down its migration speed earlier compared to a non accreting planet. Therefore, the resulting planetary systems highly depend on the timescales of migration, gap formation and gas accretion.

Another dynamical effect can play an important role in the evolution of the planets. Previous hydrodynamical studies show that multiple planets can be captured in resonant chains during the gas phase of the disk (e.g., Baruteau & Papaloizou 2013; Pierens et al. 2014; Kanagawa & Szuszkiewicz 2020). The capture in resonance can have an important impact on the migration behavior of the planets. For example, in the Grand Tack scenario (e.g., Masset & Snellgrove 2001; Walsh et al. 2011; Pierens et al. 2014), Jupiter and Saturn are believed to migrate inward and then outward due do their capture in resonance. This outward migration is occurring for precise disk parameters and mass ratios, as shown in Pierens et al. (2014). As migration can be altered by the capture in resonance, gas accretion and gap formation will also be indirectly altered by the planet radial motion.

The capture in resonance can lead to a slight increase of the planets eccentricity. As an eccentric planet will open a less deep gap for the same mass (e.g., Hosseinbor et al. 2007; Sánchez-Salcedo et al. 2023), the planet will end up with higher mass accretion rates than our planets on circular orbits. Having a higher accretion rate will just result in a shift in the different timing but this impact is negligible for low eccentricities. Since the eccentricity of the planets is damped while the amount of gas is significant in the disk (e.g., Moorhead & Ford 2009; Bitsch & Kley 2010), and since our planets do not grow massive enough to excite their eccentricity (more than 5 M_J; e.g., Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch & Kley 2010; Bitsch et al. 2013), we expect the impact of realistic low eccentricities on our results to not be significant.

We plan on implementing the impact of both migration and capture in resonance on the growth of our two planets in follow-up studies. We expect that, if the planets start accreting simultaneously, then the structure of the resulting system highly depends on the timing at which gap opening will occur because it will slow down the migration of the planets and determine their accretion behavior. A potential outward migration can delay the depletion of the inner disk, altering the mass ratio behavior discussed in Sect. 3.3. However, due to the high interdependence of each mechanism, it is difficult to precisely predict how the planets will behave.

7.4 Implications for the Grand Tack scenario

In the Grand Tack scenario, if Jupiter and Saturn have a mass ratio between 0.25 and 0.5, then they can migrate outward to their current locations (e.g., Masset & Snellgrove 2001; Crida et al. 2009; Pierens et al. 2014). In Sect. 5.1, we show that if the planets start accreting simultaneously, they reach mass ratios that are between 0.7 and 1.3. Therefore, in order to trigger outward migration, the planets have to start accreting with a delay. It also requires that the inner planet is more massive than the outer planet; otherwise, the torques arising from the outer disk would be too large, leading to inward migration. This also implies that Jupiter must have started to accrete gas efficiently before Saturn.

Pierens et al. (2014) find that outward migration depends also on both the period ratio of the planets and the disk parameters. In order to trigger outward migration in a low mass disk, a capture in a 2:1 resonance is needed. If the disk is more massive, then the planets need to reach the 3:2 resonance. Both scenarios require a relatively low α-viscosity (α ≲ 10⁻³).

Within our current parameter study, we only investigated the impact of delayed accretion in the 2:1 period ratio configuration. At low viscosity (α = 10⁻³, left panel of Fig. 10), the conditions are barely met for the outward migration to occur: in all cases, the mass ratios quickly evolve in the disk, making the planets barely stay in the needed mass ratio range (marked by the gray area). Therefore, with our current results, it seems that outward migration of the two giant planets is very challenging to reach as this occurs during a very short timescale. However, we plan to expand our parameter space study in the near future, allowing us to better analyze if and how a planetary system like Jupiter and Saturn could have formed via the Grand Tack scenario.

7.5 Comparison to exoplanets

Considering our current parameter space study, planets accreting from the same disk end up with very similar planet masses. Delaying the accretion of the respective planets allowed us to slightly broaden the mass ratio range reached; however, in 10⁵ yr, the final mass ratios obtained are still between 0.8 and 2. These mass ratios are quite different from the mass ratios observed in different planetary systems. In Fig. 11, we compare the evolution of our mass ratios to different exoplanetary systems. The data originate from the NASA exoplanet archive³. We selected the planetary systems as follows: first of all, we are interested in systems containing exactly two giant planets (i.e., with m_p > 0.3 M_J) as we investigate the accretion of two planets in the runaway gas accretion phase. Our simulation considers planet formation in a disk orbiting the Sun; therefore, the selected planetary systems orbit Sun-like stars (4700 K < T_eff < 6500 K and log(g_*) > 4). Each panel represents the ratio of the outer planet mass divided by the mass of the inner planet like in previous figures as a function of the planet period ratio. Vertical dashed lines represent the investigated period ratios. The colorbar shows the sum of the planet masses in M_J.

From the top panel in the figure, it is clear that the observed planetary systems hosting two giant planets have a broad range of mass and period ratios. We highlight the systems hosting at least one hot Jupiter (i.e., planets with periods shorter than 10 days) with a thick black contour. These planets might be formed after a very efficient inward migration or via a scattering event, leaving them very close to their host star. This results in a system where the period ratio of the planet is very large. As we do not implement migration in this study, we focus on the planets located closer to each other. Interestingly, it appears that planets located closer to each other seem to have more similar masses (except for 3 systems with mass ratios higher than 8). For better readability and comparison with our simulation, a zoom on the planets placed in the gray rectangle is shown in the lower panel of the figure.

Due to computational constraints, we could not investigate the evolution of the mass ratio until the end of the disk lifetime. Therefore, in order to make the comparison with fully formed planets as the ones observed in the different planet surveys, we show the maximal and minimal mass ratios obtained by our simultaneously accreting planets with the pink, red and purple vertical lines. As we selected the observed systems by considering planets with m_p > 0.3 M_J, we show the mass ratio spread once both planets reached 0.3 M_J. The ratios obtained for different viscosities are slightly offset from the exact period ratio for visibility. At low viscosity and small period ratios, the mass ratio range was too small to be represented by a line; therefore, the final mass ratio is shown with a square. Over the 45 observed planetary systems (including Jupiter and Saturn and hot Jupiters), 5 have a mass ratio lower than 0.7, 11 systems have a mass ratio lying between 0.7 and 1.3 and 29 systems have a mass ratio higher than 1.3. The majority of the systems cannot be explained by the simultaneous accretion of the planets. We also note that very few systems (only 5 here) feature a mass ratio as seen in our Solar System with Jupiter and Saturn. While this might be explained by the difficulty of our current facilities to see low mass planets, it also raises the question of the peculiarity of the Solar System among other systems, namely whether our planetary system common or an outlier.

In this study, such a large spread in mass ratio was only reached when the planets accrete in the runaway gas accretion phase with different delays. This accretion delay can be justified by the dependence of the beginning of the runaway gas accretion phase on the disk characteristics. Here, we followed the classical core accretion model, where it is assumed that runaway gas accretion is triggered when the planetary core reaches a mass of 20 M_⊕, with a solid core of 10 M_⊕ surrounded by a first gaseous atmosphere of 10 M_⊕ (Pollack et al. 1996). However, more recent studies show that the initial total mass of the core can vary depending on the local properties of the disk. Runaway gas accretion can be triggered at different core masses depending on the local disk temperature and opacity (Ikoma et al. 2001; Piso & Youdin 2014; Bitsch & Savvidou 2021). Moreover, in the pebble accretion scenario, the pebble isolation mass corresponds to the mass at which the core is shielded from the pebble flux by the pressure bump created by its own gap (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018). Then, the atmosphere of the planetary core is not heated anymore by the accretion of solids and cools down, entering the runaway gas accretion phase. The pebble isolation mass depends on the disk’s aspect ratio, α-viscosity, pressure profile and turbulent diffusion of the particles (Bitsch et al. 2018). Therefore, the delay of accretion time between two giants highly depends on the local properties of the disk. Protoplanetary disks can be flared, featuring large aspect ratio variations (e.g., Bitsch et al. 2015; Pierens & Raymond 2016), and different hydrodynamical properties can lead to important radial variations of viscosity (e.g., Flock et al. 2011; Dullemond & Penzlin 2018; Delage et al. 2022). The disk properties can lead to the delay of either the inner or the outer planet.

Even when we applied different accretion delays, the long-term trends could only reproduce mass ratios lying between 0.8 and 2. With these parameters, the only remaining way to reach large (or small) mass ratios, is to stop the accretion at a given mass ratio. The timing of the dissipation of the gas disk can be crucial here. For example, photoevaporation can dissipate the gas disk from inside out by creating a inner hole separating the inner disk (r < 1 AU) from the outer disk and quickly depleting it (for a review, see Ercolano & Pascucci 2017). Such a depletion of the disk might have the capacity to starve the giant planets, influencing the evolution of their mass ratio.

Some differences between our simulations and the observations might originate from dynamical events that are not yet included in our study. For example, dynamical interactions such as collisions might change the planets mass ratio (e.g., Juric & Tremaine 2008; Raymond et al. 2009; Sotiriadis et al. 2017; Bitsch et al. 2020). We also note that planets might also be ejected from the systems after the dispersal of the gaseous disk: while the observed final system would host two giants, the gaseous protoplanetary disk where they formed could have hosted three giants, changing the accretion behavior compared to what we simulate in this paper. However, accurately describing the hydrodynamical evolution of the gas together with the dynamical evolution of more than two planets for a long time (over several Myr) is quasi-impossible nowadays. Therefore, further studies should slowly including more dynamical effects during the gas phase or slowly improving the hydrodynamical assumptions of N-Body simulations.

To summarize, the simultaneous runaway gas accretion of fixed planets cannot explain the distribution of exoplanetary masses observed as it leads to planets with very similar masses. To increase the difference in planet masses, the accretion between the planets have to be delayed and efficient disk dispersal mechanisms are required to end the growth of the planets at given mass ratios. These two last points highly depends on the disk local properties of the gas. Moreover, taking dynamical interactions between the planets themselves and the planets and the disk might improve our understanding of the giants distribution. A future study including dynamical interactions and migration is planned to determine how the results can be impacted.

Fig. 11 Mass ratio evolution of the two accreting planets as a function of their period ratio compared to exoplanetary systems. The investigated period ratios from our simulations are marked by vertical dashed gray lines. The data are taken from the NASA exoplanet archive, for which we selected the systems as follows: the system contains exactly two detected planets, both of them larger than 0.3 MJ. They orbit a single Sun-like star (4700 K < Teff < 6500 K and log(g*) > 4). The color of the dots represents the sum of the planet masses in Jupiter masses, without error bars. A black contour surrounds planets that are considered hot Jupiters (i.e., with a period of less than 10 days). As we expect their formation to be highly influenced by the dynamics of the system, which we do not model here, we focus the comparison with the exoplanets on the planets marked by the gray area in the top part of the figure. A zoomed-in view of this region is shown in the second panel. The vertical red and pink lines represent the maximum and minimum mass ratio reached in each of our simulations once both planets reach 0.3 MJ. Darker colors represent lower viscosities. For visibility, the lines corresponding to the different α-viscosities are slightly offset from the period ratio line. The extent of the α = 10−4 was so small at low period ratios that we represent it with squares. The cross marks the Jupiter and Saturn couple. Simultaneously accreting planets lead to planets that are very similar in mass compared to the exoplanet population. Some systems seem to be consistent with simultaneous accretion; however, another mechanism is needed to explain the existence of the other systems. Delayed accretion as shown in Fig. 10 coupled with different disk lifetimes could explain the difference in planetary masses.

8 Conclusions

In this paper we investigate the mass distribution of two accreting planets located in the same disk. Using 2D hydrodynamical simulations, we monitor the evolution of the planetary mass ratio for different disk viscosities, different planet configurations, and different accretion timings. Our main conclusions can be summarized as:

1. The evolution of multiple accreting planets is mainly governed by the viscosity of the disk. The mass ratio evolution of simultaneously accreting planets depends on the balance between the gas accretion and gap opening timescales. As shown in Bergez-Casalou et al. (2020), at high viscosity, when gas accretion acts in favor of gap formation, the inner planet is rapidly starved by the viscous accretion onto the star and the outer planet accretes more until becoming more massive than the inner planet.

   However, at lower viscosities, when gap formation is only dependent on the disk reaction time (i.e., when gas accretion does not help gap formation), the evolution of the mass ratio of the planets follows a different behavior: the outer planet accretes more gas until the inner planet forms its gap. Then, the inner planet starts depleting the inner disk and the material present in between the planets, resulting in a higher accretion rate for the inner planet than for the outer one. When the amount of material located in these two regions is significantly depleted, the inner planet becomes starved by the outer planet.

2. Simultaneously accreting planets always end up with similar masses, independent of the disk viscosity. In order to reach more extreme mass ratios, we simulated a delayed accretion of the inner or outer planet in one configuration at high and intermediate viscosities. While the initial mass ratios are large by construction, the planets quickly tend toward similar mass ratios (0.8 < m_out/m_in < 2 in 10⁵ yr).

3. Via comparisons with the observed exoplanet population, we conclude that gas accretion occurring at the same time can explain the characteristics of only a few planetary systems. Delayed accretion coupled with different disk lifetimes leads to mass ratios that are more consistent with the observations. While core formation timescales and different disk dissipation mechanisms can explain the possibility of a delayed accretion and depletion of gas, our study shows that the majority of the observed systems of multiple gas giant planets should have started to accrete at different times.

Understanding how material is distributed between multiple planets is crucial to better understanding the dynamical evolution of the forming system. As discussed in Sect. 7.3, the radial evolution of multiple planets is governed by the migration and capture in resonance of the planets, themselves dependent on the gas distribution in the disk, which in turn is governed by gap formation and by planetary and stellar gas accretion. Future studies investigating the growth of multiple giant planets should both consider that the gas accretion of the planets is impacted by the presence of neighboring planets and provide mechanisms that explain the spread in mass ratios observed in distant exoplanetary systems as well as our own Solar System.

Aacknowledgements

C.B-C. and B.B. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. S.N.R. thanks the CNRS’s PNP program.

Appendix A Gas accretion routine

As in Crida & Bitsch (2017) and Bergez-Casalou et al. (2020), gas accretion is modeled using the Machida et al. (2010) and Kley (1999) principles. The Machida et al. (2010) accretion rate is derived by fitting 3D isothermal shearing box simulations and corresponds to the runaway gas accretion phase of the giant planet. This accretion rate is different than in the first slow contraction phase, as predicted by the core accretion model (Piso & Youdin 2014). It can be written as $${\dot{M}}_M=<\text{Σ}{>}_{0.9r_H}{H}^2{\text{Ω}}_p\times min\left[0.14;0.83{\left(r_H\text{/}H\right)}^{9/2}\right],$$(A.1)

where r_H = r_p(m_p/3M_*)^1/3 is the Hill sphere of the planet, $<\sum {>}_{0.9r_H}$ is the averaged local surface density from the planet location up to 0.9r_H; H is the disk scale height and Ω_p is the Keplerian orbital frequency of the planet.

On the other hand, in order to make sure that the planetary accretion is limited by what the disk can provide, we consider that the maximal accretion rate of the planet is given by the Kley (1999) principle. This arbitrary accretion rate considers that the planet can accrete a given fraction of the gas present in its Hill sphere. This fraction of gas is dependent on the distance to the planet and is given by f_red. This accretion rate can be written as $${\dot{M}}_K={\displaystyle \int {\displaystyle {\int }_{A_{disk}}{f}_{\text{red}}\left(d\right)}\text{\hspace{0.17em}\hspace{0.17em}Σ}\left(r,\varphi ,t\right)\pi f_{\text{acc}}drd\varphi },$$(A.2)

where A_disk is the area of the disk, f_red is a smooth reduction function predicting the fraction of gas accreted by the planet as a function of the distance from the planet d, Σ (r, ϕ, t) is the gas surface density around the planet and f_acc is the inverse timescale upon which the accretion is occurring. f_red is the smooth reduction function used to predict what fraction of gas can be accreted on the planet as a function of the distance to the planet d. It is defined as $$f_{\text{red}}=\{\begin{array}{ll}2\text{/}3\hfill & \text{if\hspace{0.17em}}d<0.45r_H\hfill \\ 2\text{/}3\times {\cos }^4\left(\pi \left(\frac{d}{r_H}-0.45\right)\right)\hfill & \text{if\hspace{0.17em}}0.45r_H<d<0.9r_H\hfill \end{array}.$$(A.3)

This function is based on Robert et al. (2018), where the authors assume that close to the planet gas accretion is 100% efficient (f_red = 1). However, a 100% efficiency is not realistic in such a case, as shown by Schulik et al. (2019). The accreted mass fraction increases closer to the planet but it does not accrete 100% of the gas in the vicinity. Thus, we chose to limit our study to a maximum f_red value of 2/3.

By combining these two methods, the final amount of gas that the planet can accrete is determined by Ṁ_p = min(Ṁ_M, Ṁ_K). The gas is removed from the disk and added to the mass of the planet. To do so, we remove the gas in this regime with the same formalism as for the Kley (1999) method. This means that if the planet is accreting in the regime of Machida et al. (2010), the total amount of gas it accretes is given by Eq. A.1 and the distribution for where the gas is removed is given by Eq. A.3.

This way, the removal scheme of the gas is the same for both principles, but the mass that can be accreted is limited either by the derived accretion rates of Machida et al. (2010) or by the maximum amount the disk can provide, given by the Kley (1999) method. In this study, we remain in the regime where Ṁ_M < Ṁ_K throughout, meaning that we always remove the amount of mass dictated by Machida et al. (2010).

Appendix B Surface density maps

In Sect. 3.1, we investigate the influence of the disk viscosity on the evolution of the planets growth. At low viscosity (α = 10⁻⁴ and h = 0.05), the RWI (Lovelace et al. 1999; Li et al. 2001) is triggered at the edges of the different planet gaps, creating vortices. In Fig. B.1, we show the 2D perturbed surface density maps of the disk hosting two accreting planets in the 3:1 period ratio, at three different times, increasing from left to right. Polar plots of the density maps presented in the top row are shown in the bottom row. The vortices produce asymmetric over-densities that can be used to trace them (in yellow in Fig. B.1).

At the beginning of the simulation (at 500 inner planet orbits, left panel), we see that vortices are produced at three different locations in this configuration: at the outer edge of the outer planet gap, in between the planets and interior to the inner planet gap. Their presence impacts the flow of gas in the vicinity of the planets, creating oscillations in the planetary accretion rates (see Sect. 3.1). Quickly, the vortices located in between the planet and in the inner disk vanish (middle panel). The strongest vortex (i.e., with the largest over-density) is the one located at the outer edge of the outer planet gap. It takes longer to dissipate and is completely vanished after 10⁵ yr (right panel). The strength of the vortices depends on the growth timescale of the planets (Hammer et al. 2017): if the planets accrete faster, the vortices would be stronger but would also vanish faster.

Fig. B.1 2D perturbed surface density maps at three different times: t = 500 (left), 3 000 (middle), and 9 000 (right) inner planet orbits. The over-densities (yellow asymmetries) are representative of vortices.

References

All Tables

All Figures

Fig. 1 Planetary accretion rates (top row) and masses (bottom row) as a function of time for different α-viscosities. Here the planets are fixed in a 3:1 period ratio: the inner planet is shown in red and the outer one in blue. As in Bergez-Casalou et al. (2020), the accretion rates slightly increase with increasing viscosity. The oscillations in the accretion rate at low viscosities are due to the presence of vortices. The inner and outer planets display similar accretion rates, leading to similar planet masses.

In the text

Fig. 2 Same as Fig. 1, but for different aspect ratios, increasing from left to right, with α = 10−3. We note on the top of each row the difference in kinematic viscosity caused by the change in the disk’s aspect ratio (v = αh2r2ΩK). The differences in the planetary accretion rates originate from the structure of the accretion formula itself (see Eq. (2)). At a low aspect ratio (h = 0.03, left panel), the accretion starts immediately in the second Machida regime. On the other hand, for h = 0.07, the change of accretion regime occurs at a higher mass and therefore at a later time, leading to a very different planet mass evolution.

In the text

Fig. 3 Mass ratio (top row) and perturbed surface density proflies (bottom row) as a function of time for different α-viscosities. As in Fig. 1, the planets are fixed at the 3:1 period ratio positions, represented by the two vertical dotted lines in the bottom row. The horizontal dotted lines in the bottom panels mark the ∑planet/∑disk = 0.1 threshold, where we consider that a gap is opened (Crida et al. 2006). In the top panels, a decreasing (respectively increasing) mass ratio indicates that the inner (respectively outer) planet is accreting more than the other planet. Different flips are observed in the evolution of the mass ratios, marked in circled numbers. The color of the surface density profiles shown in the bottom row corresponds to the color of the different vertical lines in the top row and represents different times. The last snapshot of the perturbed surface density is taken 2.4 × 105 yr (20 000 orbits of the inner planet) after the last mass ratio flip, shown by the vertical dotted black line. The first flip ① can be explained by the accretion formula itself (see Eq. (2)), while the other two flips correspond to the evolution of the surface density.

In the text

Fig. 4 Same as Fig. 3 but for different aspect ratios, with α = 10−3. As in Fig. 2, we note at the top of each row the difference in kinematic viscosity caused by the change in the disk’s aspect ratio. Here the mass ratio evolutions differ in amplitude, but they present the same behavior as in Fig. 3. The difference in amplitude originates from the direct dependence of the accretion routine on the aspect ratio, and not on the α-viscosity (see Eq. (2)).

In the text

Fig. 5 Comparison between single accreting planets and two simultaneously accreting planets. The planets are fixed at positions corresponding to the 3:1 period ratio. As in Figs. 1 and 2, the top row presents the accretion rates as a function of time and the middle row the planet masses. In the bottom row, we show the ratio of the masses in the single- and two-planet case: the red line represents the ratio of the inner planets (m2p,in/m1p,in), and the blue line represents the outer planets’ ratio (m2p,out/m1p,out). At low viscosity (left panels), the differences with the single-planet case originate from the additional formation of vortices in between the planets, which enhances the accretion rates during the vortex lifetimes. At high viscosity (right panel), while the outer planet is slightly impacted by the presence of the inner planet, we observe the starvation of the inner planet. We expect to see similar behavior at lower viscosities but delayed in time, due to longer viscosity timescales.

In the text

Fig. 6 Comparison of the mass ratios in the single-planet and two-planet cases as a function of time and for different viscosities. Again, the planets are fixed at the 3:1 period ratio positions, as shown in the upper-left corner. In solid lines we show the mass ratios of the outer over the inner planet in the two-planet case, whereas the dashed lines represent the mass ratio of the single planets. Different flips are marked by circled numbers and can be explained by the impact that gap opening and depletion of the inner disk have on the accretion behavior of the planets.

In the text

Fig. 7 Influence of a second accreting planet on the stellar accretion rate at the inner edge of the disk (0.2 AU). As in Figs. 5 and 6, the planets are located in the 3:1 period ratio positions, with an increasing α-viscosity from left to right. The comparison is made between the two-planet case (solid black line), the single-planet cases (dashed purple line for the single inner planet and dotted cyan dashed line for the single outer planet) and a disk without planets (dotted dashed gray line). The oscillations present in the case of the two accreting planets are due to the overlap of their spiral arms, coupled with the presence of vortices at low viscosity. This sometimes results in a negative radial velocity, explaining the dip seen at α = 10−4 due to the logarithmic scale. Compared to a disk with no planet, the accretion onto the star is only reduced by up to a factor of 3 in the presence of multiple giant planets, similar to the scenario with a single planet (Bergez-Casalou et al. 2020).

In the text

Fig. 8 Mass ratio as a function of time for different α-viscosities and different period ratios. The darker the line, the farther away the planets are located from each other. The separation between the planets has a small impact on the range of their mass ratio, meaning that planets that start accreting at the same time will have similar masses. The behavior of the mass ratio evolution is not impacted by the planet separation: at all viscosities the mass ratio flips are due to the same reasons described in Fig. 3. The exception is for the 2:1 period ratio at high viscosity (right panel): here the planets are close enough and the viscosity is sufficiently high to maintain a significant flow toward the inner disk through the gaps, preventing the inner planet from being starved. This results in an additional mass ratio flip, ⦷, similar to the high viscosity single-planet case of Fig. 6.

In the text

Fig. 9 Gap-opening mass as a function of the viscosity for different criteria and our simulations, as in Bergez-Casalou et al. (2020). In the left panel, we show the gap opening masses of the inner planets in the two-planet case compared to the single planet cases located at the inner position. The outer planet cases are shown in the right panel. In each panel, the lines represent the different gap-opening criteria from the literature: Crida et al. (2006) in dashed blue, Fung et al. (2014) in solid orange, and Gyeol Yun et al. (2019) in dashed-dotted red. As shown with the schematics above the panels, the two-planet cases are represented by triangles (upward for the inner planets and downward for the outer planets), and the single planets are represented by circles. The inner single-planet case is shown in gray, and the colored circles correspond to the colors of the outer positions in each configuration. For clarity, in the right panel, the viscosity is divided into discrete intervals, allowing us to compare the gap opening mass at one given viscosity in the different configurations. The gap opening masses of the inner planet are barely impacted by the presence of a second outer planet. On the other hand, at low viscosity, the closer the planets are to each other, the higher the gap opening mass is compared to the single-planet case.

In the text

Fig. 10 Mass ratio as a function of time for different accretion delays on the inner (purple lines) or the outer planet (green lines). The planets are located in the 2:1 period ratio configuration. The two panels represent two different viscosities: α = 10−3 in the left panel and α = 10−2 in the right. Darker lines represent shorter delays. Note that this time the mass ratio is displayed on a logarithmic scale (previous plots are on a linear scale) due to the extreme mass ratios induced by our initial setup here. The dots mark the moment when both planets have reached at least 0.3 MJ and can be considered gas giant planets. The gray rectangle represents the region where the mass ratio lies between 0.25 and 0.5 and shows the region where Jupiter and Saturn meet the required conditions to enter common outward migration, for α ≲ 10−3 (Masset & Snellgrove 2001; Pierens et al. 2014). Independent of the viscosity, all the mass ratios quickly tend toward the mout/min = 1 line, leading to similar planet masses in both cases (0.8 < mout/min < 2 after 105 yr).

In the text

Fig. 11 Mass ratio evolution of the two accreting planets as a function of their period ratio compared to exoplanetary systems. The investigated period ratios from our simulations are marked by vertical dashed gray lines. The data are taken from the NASA exoplanet archive, for which we selected the systems as follows: the system contains exactly two detected planets, both of them larger than 0.3 MJ. They orbit a single Sun-like star (4700 K < Teff < 6500 K and log(g*) > 4). The color of the dots represents the sum of the planet masses in Jupiter masses, without error bars. A black contour surrounds planets that are considered hot Jupiters (i.e., with a period of less than 10 days). As we expect their formation to be highly influenced by the dynamics of the system, which we do not model here, we focus the comparison with the exoplanets on the planets marked by the gray area in the top part of the figure. A zoomed-in view of this region is shown in the second panel. The vertical red and pink lines represent the maximum and minimum mass ratio reached in each of our simulations once both planets reach 0.3 MJ. Darker colors represent lower viscosities. For visibility, the lines corresponding to the different α-viscosities are slightly offset from the period ratio line. The extent of the α = 10−4 was so small at low period ratios that we represent it with squares. The cross marks the Jupiter and Saturn couple. Simultaneously accreting planets lead to planets that are very similar in mass compared to the exoplanet population. Some systems seem to be consistent with simultaneous accretion; however, another mechanism is needed to explain the existence of the other systems. Delayed accretion as shown in Fig. 10 coupled with different disk lifetimes could explain the difference in planetary masses.

In the text

	Fig. B.1 2D perturbed surface density maps at three different times: t = 500 (left), 3 000 (middle), and 9 000 (right) inner planet orbits. The over-densities (yellow asymmetries) are representative of vortices.
In the text