How planetary gas accretion changes the shape and depth of gaps in protoplanetary discs

It is well known that giant planets open deep gaps in their natal gaseous protoplanetary discs. However, it is unclear how gas accretion onto growing planets influences the shape and depth of their growing gaps. We perform isothermal hydrodynamical simulations with the Fargo-2D1D code in which planets accrete gas within full discs that range from 0.1 to 260 AU. The gas accretion routine uses a sink cell approach, where we use different accretion rates to cope with the large span of gas accretion rates in the literature. We find that the planetary gas accretion rate increases for larger disc aspect ratios and larger viscosities. Our main result shows that gas accretion has an important impact on the gap opening mass: we find that when the disc responds slowly to a change in planetary mass (i.e. at low viscosity), the gap opening mass scales with the planetary accretion rate, with a higher gas accretion rate resulting in a larger gap opening mass. On the other hand, if the disc response time is short (i.e. at high viscosity), then gas accretion helps the planet carve a deep gap. As a consequence higher planetary gas accretion rates result in smaller gap opening masses. Our results have important implications for the derivation of planet masses from disc observations: depending on the planetary gas accretion rate, the derived masses from ALMA observations might be off by up to a factor 2. We discuss the consequence of the change of the gap opening mass on the evolution of planetary systems by taking the example of the Grand Tack scenario. Planetary gas accretion also impacts the stellar gas accretion, where we find only a small influence due to the presence of a gas accreting planet.


Introduction
Recent ALMA observations show protoplanetary discs with numerous different features in the gas (Teague et al. 2018;Pinte et al. 2020) or in the dust (Andrews et al. 2018). An important question is to determine if these features can be explained by the presence of planets. Recent exoplanet statistics show that most stars host planets around them (Mayor et al. 2011;Howard et al. 2012;Suzuki et al. 2016;Mulders et al. 2018), but it can be difficult to observe them directly as they might be extinct by the disc itself (Sanchis et al. 2020). One should rely then on the impact of the presence of these planets on discs. Several studies showed that planets are able to produce the features observed in these discs (Pinilla et al. 2012). For example, Zhang et al. (2018) shows that the gaps and rings from the DSHARP survey could be explained by the presence of gap opening planets in the disc.
The formation of gaps by giant planets in protoplanetary discs in the dust (Paardekooper & Mellema 2006) and in the gas (Lin & Papaloizou 1986;Crida et al. 2006;Fung et al. 2014;Kanagawa et al. 2015) has been studied in the past in hydrodynamical simulations. Different criteria for the gap opening mass (i.e. planetary mass needed to open a gap of a given depth) have been derived as a function of disc parameters. The main motivation for these studies was in fact the migration of the giant planets that can open these deep gaps. However, only recently have these migration studies started to include planetary gas accretion (Dürmann & Kley 2015;Crida & Bitsch 2017;Robert et al. 2018).
Gas accretion is a complicated problem that requires 3D high resolution hydrodynamical simulations. The core accretion model (Pollack et al. 1996) suggests that after forming a solid core, a gaseous envelope is slowly accreted until the planet reaches a critical mass where the mass of the gas envelope is equal to the mass of the core. The planet then enters a runaway growth phase of gas accretion where we can assume that only gas is accreted (Hubickyj et al. 2005;Lissauer et al. 2009). Previous studies described how gas accretion can be modeled in different frameworks: in 2D with and without migration (Crida & Bitsch 2017;Kley 1999) or in 3D (Ayliffe & Bate 2009;Machida et al. 2010;D'Angelo & Bodenheimer 2013;Lambrechts et al. 2019;Schulik et al. 2019). These complex hydrodynamical simulations are needed to understand each phase of gas accretion. Due to their high computational cost, they are often integrated over short timescales. In fact a lot of these recent simulations can only cover 100 planetary orbits (or less), making it unfortunately impossible to study the long term evolution of these systems within the 3D high resolution framework with the current computing systems.
As gas accretion requires a computationally expensive resolution, previous studies on gap opening by giant planets neglected the planetary gas accretion. Our goal in this paper is to apply previously found gas accretion rates (Kley 1999;Machida et al. 2010) in a 2D isothermal disc framework to study the long term behavior of the growth of the planet and its influence on the shape and depth of the created gap.
Article number, page 1 of 17 arXiv:2010.00485v1 [astro-ph.EP] 1 Oct 2020 A&A proofs: manuscript no. aanda This work is structured as follows. In section 2 we describe our numerical set-up with the description of the FARGO-2D1D code ) and our gas accretion routine. In section 3 we study the difference between an accreting planet and non accreting planets, the impact that different accretion rate can have on the gap shape and on the accretion onto the star. Then different disc parameters (h and α) are explored in section 4, as well as the impact of gas accretion on the gap opening mass and migration of other planets that could form in these discs. We finish by some discussion on the consequences for observation and planetary system structures in section 5, before concluding in section 6.

The FARGO-2D1D code
In this paper, we simulate an accreting planet on a fixed circular orbit embedded in a viscous disc with the hydrodynamical code FARGO-2D1D ). This code allows us to simulate the whole viscous evolution of a protoplanetary disc for a reasonable computational cost by surrounding a 2D-grid (r, φ), where the planet is located, by two 1D-grids azimuthally symmetric as shown in figure 1. The extent of the grids is chosen so that the 1D-grids are far enough from the planet, where we can consider the disc axisymmetric. The 1D inner disc spans from 0.1 AU to 0.78 AU. The 2D grid then ranges from 0.78 AU to 23.4 AU, with a 1D outer disc from 23.4 AU to 260 AU. The resolution is such that there are five cells per Hill radius of the planet at t=0, which leads to N r = 802 and N φ = 1158 (dr/r 0 dφ 0.005). Considering that the resolution is fixed in time, the Hill sphere region will be more and more resolved as the planet grows (r H ∝ m 1/3 p ). For computational accuracy, the code uses dimensionless units. We normalized masses with the mass of the central star M 0 = M , lengths with the position of the planet r 0 = 5.2 AU and set the gravitational constant to be G = 1. We end up with a time normalized by the orbital time at r 0 , t 0 = (r 3 0 /(GM * )) −1/2 , meaning that the orbital period of the planet is P = 2πt 0 .
The disc is locally isothermal and its aspect ratio h = H/r is constant. Three different values are investigated in section 4.1, h = 0.03, 0.05 and 0.07. The density profile is defined by Σ(r) = Σ 0 × (r/r 0 ) −1 , where Σ 0 = 3.10 −4 = 93.6 g/cm 2 at r 0 . Its value is chosen so that the total mass of the disc is 0.1M * , corresponding to a heavy disc due to its large radius (Baillié et al. 2019). This large radial extent allows us to neglect self-gravity. The disc is subject to an α-viscosity as described by Shakura & Sunyaev (1973). We investigated the influence of the viscosity in section 4.2 by taking different values for α = 2×10 −2 , 10 −2 , 10 −3 and 10 −4 . Our fiducial setup is shown in the first row of table 1, where we summarize the parameters of the different simulations we compare in this paper.
The accreting planet starts with a mass of m p = 20 M ⊕ or m p = 10 M ⊕ depending on the aspect ratio (see table 1). With this initial mass, we can assume that the planet will be directly in the runaway gas accretion regime (Pollack et al. 1996). The accreting planet is introduced in the disc with a mass taper function making the planet grow from 0 to its initial mass in n orb = 3 orbits. The mass taper function is written as: In order to let the disc adapt to the planet, we wait until the initial mass creates a steady gap. We define the equilibrium by a change of gap depth of less than 1% within 100 orbits. The times needed The dot represent the location of the planet and the star represent the central star. The planet, located at r p = 1 = 5.2 AU, is considered far enough from the boundaries of the 2D-grid for the disc to be considered axisymmetric.
to reach the steady state and steady gap depth are summarize in table 2 as a function of the disc parameters. In principle, the introduction of the planet in the simulation can change the resulting structure of the disc depending on how long the ramp-up time is (Hammer et al. 2017). However, this applies mostly for giant planets and low viscosities. For our initially small planets, we found no significant difference if we let gas accretion onto the planet start after longer times.

Gas accretion routine
Gas accretion is modeled using the Machida et al. (2010) and Kley (1999) principles, like in Crida & Bitsch (2017). The Kley (1999) principle is an arbitrary accretion rate, but limited by how much gas is available within the hill sphere of the planet. The Machida et al. (2010) accretion rate is based on 3D isothermal shearing box simulations and represents the accretion rate of the planet in the runaway gas accretion phase. The accretion rate is different than in the contraction phase predicted by the core accretion model (Piso & Youdin 2014).
Crida & Bitsch (2017) used Machida's accretion rate without reducing the gas surface density from its initial value Σ 0 when computing the accretion rate. To make sure that they were limited by what the disc can provide, they limited the accretion rate to the minimum between Machida's and Kley's accretion rates. As planets create gaps, the surface density profile around the planet is reduced compared to the initial profile. Therefore, in our study, we use the local surface density, defined as the average of the surface density from the position of the planet to 0.9 r H . Our accretion routine is written as follows: (2) where S disc is the disc surface, H the disc scale height, Ω the keplerian orbital period of the planet, d the distance from the planet, r H = r p (m p /3M * ) 1/3 is the Hill sphere of the planet and f acc is the inverse time scale upon which the accretion rate of Kley (1999) is occurring. Here, f acc = 1 to determine which of the accretion rate is the smallest. f red is a smooth reduction function used to predict which fraction of gas must be accreted on the planet as a function of the distance to the planet. It is defined as: This function is based on Robert et al. (2018) where they assumed that close to the planet, gas accretion is 100% efficient (f red = 1). However, it seems that 100% efficiency is not realistic close to the planet. Schulik et al. (2019) have shown that gas accretion does not start from the full Hill sphere but only from a fraction of it. The accreted mass fraction increases closer to the planet and does not accrete 100% of the gas close to it. We choose then to limit ourselves to a maximum value for f red of 2/3. As Machida's formula only gives us an information on the amount of mass that should be accreted, we remove the gas in this regime with the same formalism as for the Kley (1999) method. This means, if the planet is accreting in the regime of Machida et al. (2010), the total amount of gas it accretes is given by eq.2 and the distribution where the gas is taken out is given by eq.3. This way, the removal scheme of the gas is the same for both regimes, but the mass that can be accreted is either limited by the derived accretion rates of Machida et al. (2010) or by the maximum amount the disc can provide given by the Kley (1999) method. Unless specified, we are always in the regime wherė M M <Ṁ K , meaning that we always remove the amount of mass suggested by Machida et al. (2010).

Influence of gas accretion
In this section we first investigate the difference in gap shape between accreting and non accreting planets. Then we focus on the influence of different gas accretion rates, to explore the variety of gas accretion rates found in the literature. We finish this section by studying the influence of planetary accretion on the stellar accretion.

Accretion VS non accretion
In this section, we investigate the effect of gas accretion on the gap shape and the pressure bumps generated by the planet exterior to its orbit. We compare three different simulations: (i) a planet accreting gas following the recipe presented in section 2.2, where the mass is removed from the disc and added to the planet (solid gray line in the following plots), (ii) a mass tapered planet where the mass of the planet is changed via the mass taper function given by equation 1 and no mass is removed from the disc (dashed blue line) and (iii) a planet with a mass fixed to the value where the comparison is made (dotted dashed red line). Each simulation is compared when the planets reach the same mass m comp . For the mass tapered planet, n orb is chosen so that the planet grows from 0 to m comp in the same amount of orbits that the accreting planet needs to reach this mass. We chose to make the comparison before and after reaching gap opening mass for α = 10 −2 and h = 0.05, where the gap opening mass is the mass needed to open a gap of depth Σ/Σ unp = 0.1, as in Crida et al. (2006). We compare the structure of the disc for planets of masses equal to 0.5 M J (figure 2) and 1.5 M J (figure 3). With these parameters, n orb is equal to 580 orbits for m comp = 0.5 M J and to 4236 orbits for m comp = 1.5 M J .
On the top panel of figure 2 we show the time evolution of the planetary masses in the three different cases for m p = 0.5 M J . The mass evolution of the planet growing via the mass taper is very similar to the accreting planet. The dots show at which mass and time the comparison is made. On the second panel we show the gap profiles at the time when the accreting planet reaches 0.5 M J (marked by the dots in figure 2). The surface density profile has been normalized to the surface density of a disc without planet at the same time (Σ unp ), to get rid of the effect of natural viscous evolution of the disc. Our simulations show that the accreting planet generates a slightly deeper gap than the planet with a fixed mass, in line with the results of Dürmann & Kley (2015).
To create a gap, a planet needs to have a strong enough gravitational torque to overcome the pressure torque and viscous spreading that tend to close the gap (Crida et al. 2006). The accreting planet presents a deeper gap as material is removed from the disc compared to the other two simulations. The difference in surface density at the bottom of the gap is very similar to the difference in surface density in the inner disc, due to the fact that mass is removed partially from the inner disc. This will impact the stellar accretion, as we discuss in section 3.3. For these disc parameters (α = 10 −2 , h = 0.05), Dürmann & Kley (2015) found the same behavior (figure 2 of their work). Therefore, relatively to the edges of the gap, the gaps have similar depths, but when compared to a disc without planet, the accreting planet produces a deeper gap, as the non accreting planets only push material away from their orbit, enhancing the surface density in the inner disc while the accreting planet is removing material. On the other hand, the fixed mass planet is opening a deeper gap than the mass tapered planet, because the fixed mass planet exerts a stronger torque on the disc material due to its larger mass until the time of comparison. At later time, when the mass taper fixes the mass of the planet to 0.5 M J and the disc has time to adjust to this planetary mass, the shape of the gaps of the fixed mass and mass tapered planets is the same.
At later time, when the planet becomes heavier (m p = 1.5 M J ), we also see on the middle panel of figure 3 that the inner disc is more depleted for the accreting planet. Gas accretion has therefore an important role in the depletion of the inner disc. One of the consequence of this behavior is discussed in section 3.3. The gap depth is only slightly influenced by gas accretion in this case.
The difference in the gap depth can lead to a difference in the torques acting on the planet and thus for its migration behavior. Kanagawa et al. (2018) found that when the planet creates a gap deep enough (Σ min /Σ unp < 0.6) then the total torque felt by the Article number, page 3 of 17 A&A proofs: manuscript no. aanda Top panel: Time evolution of the planetary mass in the three different simulations: accreting planet (solid gray), mass tapered planet (dashed blue) and fixed planet mass (dashed red). Mass is removed from the disc only for the accreting one. We plotted the mass evolutions after gas accretion is allowed (after t = 100 planetary orbits 1.2 × 10 3 yrs) Middle panel: Perturbed surface density at the time when m = 0.5 M J (dots on the top panel). The horizontal gray dashed line marks the gap opening criterion as defined by Crida et al. (2006). An accreting planet will create a deeper gap than a non accreting one when the gap opening mass is not reached yet. Bottom panel: Pressure gradient in the three different cases for m=0.5 M J . Similar pressure gradients imply that gas accretion does not have a major influence on the pressure structure in the protoplanetary disc. This implies that particles are trapped at a similar position and we can thus not distinguish between accreting and nonaccreting planets via the shape of dust traps in this case.
planet is proportional to the gap depth. Thus, the difference in gap depth between the accreting and the planet with fixed mass should lead to an equivalent difference in total torque. The difference in gap depth for m p = 0.5 M J is of about 24% while the difference in the total torques measured in our simulations is ∼ 17%. Using the approach by Kanagawa et al. (2018) regarding the migration behavior and comparing it to our simulations indicate a deviation of 25% between the actually measured torques and the difference that should occur because of the difference in gap depth. This deviation might rely within the errorbars of the fit derived by Kanagawa et al. (2018). On the other hand, for the more massive case plotted in figure 3, the difference in gap depth is 38% and the difference in the total measured torques is ∼ 37%. Here the differences match and our simulations confirm the results from Kanagawa et al. (2018). Additional simulations with 1 M J planets confirm this trend. Therefore the differences in torques between an accreting planet and a planet with fixed mass can be explained by their difference in gap depth, as long as the gap is deep enough.
The difference in gap depth and shape will have an influence on the drift of dust particles in the disc. Observations (Teague et al. 2018;Andrews et al. 2018) show very often rings and gaps in the dust profile of discs, which could be explained by the presence of planets (Pinilla et al. 2012;Zhang et al. 2018), which generate pressure bumps exterior to their orbits where dust can accumulate (Paardekooper & Mellema 2006). We show on lower panels of figures 2 and 3, the pressure gradients obtained in the three different cases. As expected from the surface density profiles, the pressure gradients are too similar to show differences in observations.
In the following sections of this paper we investigate different gas accretion strength and other disc parameters influence on the depth and shape of gaps.

Influence of different gas accretion rates
Gas accretion rates onto planets are not very precisely constrained. In the runaway gas accretion phase accretion rates range fromṀ cent studies are claiming that the accretion rates could be even smaller: Lambrechts et al. (2019) found that the mass flux through the envelope is different from the gas accretion rate, therefore even if the mass flux is of the order of 10 −5 M J /yr, the gas accretion rate is actually 10 to 100 times smaller. Tanigawa & Tanaka (2016) derived an even smaller accretion rate when taking into account gap opening (Ṁ p 10 −8 M J /yr Fig. 5. Pressure gradient in different gas accretion models for m = 0.5 M J . Even though the different models present different inner disc structures at the same planetary mass, the difference in the pressure gradients is too small to be observationally resolved. The bump located at the planet position is due to the material around the planet location. It vanishes with higher accretion rates as the planet accretes all the material at this location. 3 × 10 −6 M ⊕ /yr). In order to reflect this discrepancy, we investigate different accretion rates in this subsection, by scaling our nominal Machida accretion rates (eq.2) by a factor of 0.1 to 10. In top panel of figure 4 we show the seven different accretion rates obtained when we scale the Machida accretion rate by a factor of 0.1, 0.2, 0.5, 1, 2, 5 and 10. A darker color represents a higher accretion rate. The maximum of the accretion rate is shifted to later time if the accretion rate is reduced. This maximum is reached when the accretion rate (eq.2) switches from a regime dominated by the Hill sphere to a regime proportional to 0.14. The critical mass that is needed for this switch in accretion rate is reached at later times for reduced accretion rates.
As gas accretion is limited by what the disc can provide, the accretion rates in the case of an enhancement by 5 and 10 are the same at the beginning of the simulation: they start in the Kley accretion regime (eq. 2). Both simulations switch to Machida accretion regime and become different after a certain time. This indicates that most of the time the planet is accreting less than what the disc can provide. The values of the accretion rates range from ∼ 6×10 . We note, however, that the accretion rates reduce in time, where at later stages our measured accretion rates are below 2 × 10 −5 M J /yr. This indicates that the gap opening process of planets plays a role in setting the accretion rate.
In the second panel we show the time evolution of the planetary mass in the seven different cases. The biggest planet formed after 6 × 10 4 years is slightly less than 3 M J , obviously in the case of the highest accretion rate. In addition to the assumptions of our simulation setup, the time evolution plays a crucial role. A longer integration time would lead to a larger planetary mass. In reality, a planet would also migrate to the inner disc, where the gas surface density is larger and thus the accretion rate could be higher, even if the Hill sphere of the planet shrinks while getting closer to the central star, depending on the surface density profile slope (Crida & Bitsch 2017).
In the third and fourth panels of figure 4 we show the azimuthally averaged surface density profiles. We normalized the profiles to the profile of a disc without a planet as explained in section 3.1. This way we compare the actual shape of the gap. For this viscosity and aspect ratio, α = 10 −2 and h = 0.05, gas accretion seems to have only a small effect on the gap shape, resulting in a slightly deeper gap for the most reduced accre- Fig. 6. Influence of planetary gas accretion on the stellar gas accretion for our fiducial disc parameters (h = 0.05, α = 10 −2 ). Top panel: Time evolution of the stellar gas accretion rates at the inner disc (0.2 AU). The different lines represent different simulations: with an accreting planet (solid gray), with a mass tapered planet (dashed blue), with a fixed mass planet (red dotted dashed) and without a planet (bold black dashed). We can see here that an accreting planet decreases the stellar accretion rate. Bottom panel: Perturbed surface density at three different times for the accreting case and after 4200 orbits in the mass tapered and fixed mass cases. The dashed vertical line is located at 0.2 AU, where we measure the stellar accretion rate. The three different times correspond to the dots in the upper panel. The time evolution of the perturbed gas surface density shows how the inner disc is slowly depleted by the viscous stellar accretion but also by the accreting planet. In comparison, the mass tapered and fixed mass cases present a less depleted inner disc, due to the absence of planetary gas accretion (see also figure 3). tion rate. The difference in gap depth is due to the depletion of the inner disc: for the reduced accretion rate, the planet reaches 0.5 M J at a later time, giving time to the inner disc to be depleted via viscous spreading towards the star. We discuss this in more details in the next section. As the gap shapes are only slightly influenced by gas accretion, one can expect the pressure bumps to be similar, which is shown in figure 5. The difference in the disc pressure profiles for planets with different accretion rates are very small, making it indistinguishable via observations, in contrast to migrating planets (Meru et al. 2019;Weber et al. 2019). However, in the next section we investigate the influence of planetary gas accretion on another observable: gas accretion onto the star.

Influence on the stellar accretion rate
Gas accretion onto the star is an observable feature that can be measured via the UV excess in a star SED (Hartmann et al. 1994;Calvet & Gullbring 1998). We study here the influence that a gas accreting planet has on the stellar gas accretion compared to a disc without a planet. We defined the stellar gas accretion rate as the mass flux through the inner boundary of the disc: M * = −2π r in v r,in Σ in where r in = 0.2 AU and v r,in and Σ in are the radial velocity and surface density at r in respectively. We compare in figure 6, the stellar gas accretion of a disc with an accreting planet (solid gray line), a fixed mass planet (dotted dashed red), a mass tapered planet (dashed blue) and a disc without a planet (bold dashed black). The two non accreting planets have Fig. 7. Influence of different planetary gas accretion rates on the stellar gas accretion for the fiducial disc parameters (h = 0.05, α = 10 −2 ). Top panel: Time evolution of the accretion rates onto the star at the inner disc (0.2 AU) for eight different cases: the seven different planetary accretions rates (enhanced, fiducial and reduced) and a disc without a planet (bold black dashed line). A more efficiently accreting planet decreases the stellar accretion rate more efficiently. This is due to the depletion of the inner disc induced by the planetary accretion. The black dots mark the time at which the surface density profiles are plotted. Bottom panel: Perturbed surface density at time t = 3000 orbits = 36 539 years for the seven different cases. The vertical gray dashed line lies at 0.2 AU which is the location of the inner disc, where the stellar accretion is measured. a mass m p = 1.5 M J , as in figure 3, leading to a mass taper of n orb = 4236 orbits 51600 yrs.
At the beginning of the simulations, only the fixed mass planet shows a significantly higher stellar accretion rate than the others. This is due to the fact that a massive planet of 1.5 M J was introduced in the disc in only 3 orbits. Therefore the disc is highly perturbed and the planet pushes a significant quantity of material in the inner disc. AsṀ * ∝ Σ, the stellar gas accretion is then enhanced. At the same time, the mass tapered and accreting planets are slowly growing into the disc, pushing less material to the inner disc at the beginning of the simulation.
Even though Rafikov (2016) showed that spiral arms can drive accretion via depositing angular momentum, these spiral arms are too weak at this location to make a difference. They also showed that such accretion may occur at the end of the lifetime of the disc, as it corresponds to very short accretion timescales. In consequence, the stellar accretion rates for the mass tapered planet, the fixed mass planet and the disc without a planet are similar.
On the other hand, we see that after ∼ 10 4 years, the accretion rate differs from one simulation to another. First, the fixed mass planet shows an accretion rate slightly lower than the disc without a planet. Even though the gap might prevent some material to reach the inner disc, there is still a large flow of material through the gap for the inner disc to avoid being significantly depleted. This can be seen in the middle panel of figure 3, where the surface density in the inner disc of the fixed mass case is only slightly affected by the presence of the planet. This is due to the high viscosity of this simulations (α = 10 −2 ). Before the end of the simulation, the stellar accretion rate in the case of the mass tapered planet becomes slightly higher than the accretion of the disc without a planet. In the top panel of figure 3 we see that the mass taper function makes the planet grow rapidly after ∼ 10 4 years, resulting in the same effect as for the fixed mass planet: as the planet grows rapidly, it pushes a large quantity of material in the inner disc, enhancing the stellar accretion rate, and then as the gap profile tends to reach an equilibrium, the same flow through the gap as the fixed mass planet will be reached.
On the other hand, the accreting planet reduces the stellar accretion rate. This is due to the fact that the planet, besides accreting material from the outer disc and therefore preventing a large amount of gas to flow through the disc, accretes gas from the inner disc as well. It helps depleting the inner disc quicker than by the viscous accretion naturally present in the disc. The depletion of the inner disc by planetary gas accretion can be seen on the lower panel of figure 6, with the surface density profiles of the accreting planet at different times: Σ/Σ unp < 1, as Σ unp is defined to be the surface density of the empty disc. As in figures 2 and 3 (middle panels), we see also that the depletion is enhanced by the accretion onto the planet compared to the two non accreting planets.
However, we notice that the accretion rate onto the planet can be higher than the stellar accretion rate (top panels of figures 4 and 7). The planet can accrete more than what the disc can supply only at the very beginning of the simulation as it is emptying its horseshoe region. Then, it accretes material from the inner and outer parts of its orbit: there the planetary accretion is limited to what the disc can provide. The reason why the stellar gas accretion rate is smaller than the planetary accretion rate after a longer time is because the stellar accretion rate is measured at the inner edge of the disc (mass flux at r = 0.2 AU) whereas the planet is accreting material from the inner and the outer disc (inside and outside r = 5.2 AU). We compared the stellar flux at a radius located beside the planet (i.e. at 7.5 AU while the planets is located at 5.2 AU) and we observed that the planet accretion rates become limited by the mass flux at this location after a certain time: the time needed to empty its horseshoe region and a large part of the inner disc. In conclusion, the planetary accretion is limited to what the outer disc can provide as the inner disc is depleted by the planet and the viscous accretion towards the star.
To investigate the influence of planetary gas accretion into the depletion of the inner disc, we show in figure 7 the stellar accretion rate for all our enhanced and reduced planetary accretion rates. The top panel shows that larger planetary accretion rates result in lower stellar accretion rates. As the planet accretes more gas, it depletes the inner disc faster, reducing the stellar accretion rate. At t = 36500 yrs = 3000 orbits, time shown with the black dots on the top panel, the lowest stellar accretion rate is roughly three times lower than the highest rate. However, the corresponding planetary accretion rates differ by a factor of four ( figure 4). This means that the stellar accretion does not directly scale with the planetary accretion. It is due to the fact that the depletion of the inner disc also highly depends on the viscosity of the disc. Here, the viscosity is high: with α = 10 −2 , planetary gas accretion and viscous spreading through the gap act together to deplete the inner disc. Therefore, if the planetary gas accretion rate is not high enough, then the inner disc is replenished by viscously spreading gas, whereas when the planetary accretion rate is high, the inner disc does not have time to be replenished by viscous spreading. In section 4.2 we study the influence of the viscosity on the stellar accretion and discuss the consequences in the section 5.1.
In conclusion, at high viscosity, planetary gas accretion rate has a strong influence on the stellar accretion. In the following section we investigate the influence of two disc parameters, the aspect ratio and the viscosity, that can have an influence on the accretion onto the planet and onto the star.

Influence of disc parameters
Disc parameters have a strong influence on how a planet will carve a gap (Lin & Papaloizou 1986;Crida et al. 2006;Fung et al. 2014;Kanagawa et al. 2015). Therefore we can expect them to have a strong influence on gas accretion too. Here we study the influence of the aspect ratio and the viscosity on the planetary gas accretion rate before studying the influence on the gap opening mass. Finally, as migration depends on the disc parameters too, we investigate how planetary gas accretion can alter the migration of outer planets in the disc by studying migration maps.

Different aspect ratios
The aspect ratio of the disc has two influences here. The first is related to gas accretion (eq. 2), which is for some part proportional to (r H /H) 9/2 . At the beginning of the simulation, the planet is in this regime until the Hill sphere reaches 2/3 of the disc scale height and then switches to the regime whereṀ ∝ 0.14. In Machida et al. (2010), the different regimes are explained by the switch from a regime regulated by Bondi accretion to a regime governed by Hill accretion. This switch is visible in our work in the time evolution of the gas accretion rates in figures 4, 6, 8 and 9. The aspect ratio has a direct influence on the mass at which the switch occurs. In top panels of figure 8 we show the resulting planetary accretion rates for different aspect ratios, with h = 0.03, 0.05 and 0.07. As the aspect ratio increases, the switch occurs at higher masses and therefore at later time. For h = 0.03, the simulation directly starts in the regime ∝ 0.14 (eq. 2) as the mass at which r H 2H/3 is smaller than our initial mass of 20 M ⊕ (m switch = 7.8 M ⊕ ).
We show in the lower panel of figure 8 the perturbed gas surface density profiles for the different aspect ratios and compare accreting planets with planets that have a fixed planetary mass of 0.5 M J . Gas accretion has a different impact on the gas surface density profile depending on the aspect ratio. For h = 0.07, the gap is deeper with accretion than without accretion. For h = 0.05, we find back the same result from section 3.1 where gas accretion has only a slight impact on the gap depth. Finally, for h = 0.03, the planet reached gap opening mass, making the gap look really similar in the accreting and non accreting cases. Comparing a planet with a fixed mass of 0.2 M J for h = 0.03 (below the gap opening mass defined by Crida et al. (2006)) with an accreting planet of the same mass (not shown), reveals a deeper gap in the case of the fixed mass planet. Our results show three different behaviors depending on the aspect ratio: gas accretion helps creating a deeper gap (large h), has almost no impact (h = 0.05) or prevents from creating a deeper gap (low h). These three behaviors are explained in section 4.3, where we explore the impact of gas accretion on the gap opening mass. In the low aspect ratio case, one should note that a 20 M ⊕ mass planet already creates a deep gap at low viscosities. We choose to keep this initial mass in this section to compare with the higher aspect ratios as the comparison is made at high viscosity. However, we reduce the initial mass down to 10 M ⊕ in section 4.3 to study the impact on the gap opening mass as a function of the disc parameters.
As the accretion rates are very different from one aspect ratio to the other, the planet masses are diverse too. The middle panel of figure 8 shows that a later switch in accretion rates results in As the time needed for the initial mass to reach an equilibrium is higher for lower aspect ratios, gas accretion starts later for h = 0.03 (see table 2). Second panel: Evolution of the accretion rates as a function of the time after gas accretion is turned on. Third panel: Time evolution of the planetary mass. Lower aspect ratios result in more massive planets. Bottom panel: Perturbed surface density at the time where m p = 0.5 M J (dots on previous panels). Two effects add up here making the gap shapes really different: the influence of the aspect ratio on the shape of the gap and the time at which the planet reaches 0.5 M J . We compare the gap shapes of accreting planets with gaps created by fixed mass planets: depending on the aspect ratio, gas accretion has a different influence. The horizontal gray dashed line marks the gap opening criterion as defined by Crida et al. (2006). more massive planets. As a result, our simulations indicate that planetary gas accretion is more efficient in hotter discs. As it is harder to create a deep gap in hotter discs, the surface density at the location of the planet is larger (lower panel of figure 8) and allows a more efficient planetary gas accretion. On the other hand, if the disc is hotter, it is harder for a planet to reach the pebble isolation mass (Lambrechts et  bump exterior to its orbit preventing solid accretion. In addition, if the disc is hotter, it is more difficult to accrete pebbles efficiently, because they are less concentrated in the disc's midplane (Youdin & Lithwick 2007), resulting in a less efficient formation of planetary cores (Ndugu et al. 2018). Therefore, it is easier to form giant planets via runaway gas accretion when their core is already formed but it is harder to initially form these cores in hotter discs.

Different viscosities
Another important disc parameter is the viscosity of the disc, which determines how a planet opens a gap (Crida et al. 2006;Fung et al. 2014;Kanagawa et al. 2015) and dictates how gas flows in the vicinity of the planet. Both have a strong impact on gas accretion. In order to determine the influence of viscosity, we run five different simulations with the following alpha parameters: α = 2 × 10 −2 , 10 −2 , 10 −3 and 10 −4 (table 1). As mentioned in section 2.1, we need to wait until the initial mass of 20 M ⊕ creates a steady gap, defined as a change of less than 1% of the gap depth within 100 planetary orbits. This waiting time is dependent on the disc parameters and is summarized in table 2. One can expect the planetary accretion rate to decrease with decreasing viscosity because viscosity dictates how fast gas removed by accretion can be replenished from the viscously spreading disc. This behavior can be seen on the top panels of figure 9: on the first panel the gas accretion rates are plotted as a function of time. One can see that at lower viscosities the initial mass needs a longer time to reach a steady gap. This means that at lower viscosity, gas accretion starts later. On the second panel the accretion rates are shifted so that they are plotted as a function of the time after gas accretion is turned on. The accretion rates plotted here show that a lower viscosity results in a lower planetary accretion rate.
At low viscosity, instabilities in the disc are generated (Klahr & Bodenheimer 2003;Fu et al. 2014), resulting in a threshold value for the gas surface density in the vicinity of the planet. If this threshold exists, we would expect a limit of the planetary accretion rate as well. However, in order to resolve these instabilities adequately, high resolution simulations are needed, which Fig. 9. Influence of different viscosities for h = 0.05. Top panel: Time evolution of the accretion rate onto the planet. As the time needed for the initial mass to reach an equilibrium is highly dependent on the viscosity, gas accretion starts later for lower viscosities. Second panel: Evolution of the accretion rates as a function of the time after gas accretion is turned on. As the viscosity is lowered, the Rossby Wave Instability is triggered: for α 10 −4 , vortices are formed and influence the planetary gas accretion rates, explaining the oscillations in the accretion rate curves. Third panel: Evolution of the planetary mass as a function of the time after gas accretion is turned on. The dots represent the time at which the gap opening mass is reached. Bottom panel: Perturbed surface density at the time where gap opening mass is reached (dot on the other panels). The gap opening mass is defined by the mass needed to reach Σ/Σ unp = 0.1 (Crida et al. 2006). It is represented here by the horizontal gray dashed line. are computationally very expensive. We thus limit ourselves to α ≥ 10 −4 .
For α 10 −4 , the Rossby Wave Instability is triggered (Lovelace et al. 1999;Li et al. 2001). This instability generates vortices at the location of steep density gradients. Such steep gradients are induced by planetary gaps formed at low viscosity (Hammer et al. 2017;Pierens et al. 2019), as we can see on the radial density profiles in the lower panel of figure 9. Vor- Fig. 10. Influence of different accretion rates on stellar gas accretion for α = 10 −4 and h = 0.05. Top panel: Planetary gas accretion rate for the seven reduced, enhanced and fiducial accretion rates. The presences of vortices create oscillations in the accretion rate: the larger the vortex is, the larger the oscillations are. Middle panel: Stellar gas accretion at the inner edge of the disc (0.2 AU). At this low viscosity the trend is flipped compared to the high viscosity case: a more efficiently accreting planet will lead to an enhanced accretion onto the star. This is due to the fact that at low viscosity the inner disc is mostly perturbed by the presence of the planet and the viscosity is too low to compensate it. Bottom panel: Perturbed surface density at time t = 4000 orbits = 48 718 yrs. The vertical line shows the location at which the stellar accretion rates are measured (0.2 AU). tices will modify the flow of material into the gap, changing the amount of gas available for accretion by the planet. Such impact can be seen in the accretion rates. In the top panel of figure 9, for α = 10 −4 , the accretion rate is oscillating due to the presence of a vortex at the outer edge of the gap. We show in figure 11 the 2D (r, φ) surface densities at three different times. The fiducial case is plotted in the middle row. The asymmetrical overdensity (in yellow) located at the outer edge of the gap is characteristic of the presence of the vortex. At later times, the vortex vanishes. In the fiducial case, it completely vanishes after t 8.1 × 10 4 years, inhibiting the oscillations of the accretion rate (top panel of figure 9).
As the presence of vortices is linked to the steepness of the density gradient, the characteristics of the created vortices will depend on how fast the planet will grow and create its gap. Hammer et al. (2017) find that a slowly growing planet will create a weak vortex at the outer edge of the planet gap. We find a similar result when we apply our different accretion rates to simulations with α = 10 −4 . In figure 11 we show the 2D surface density snapshots for three different accretion rates (reduced by 10, fiducial and enhanced by 10 from top to bottom) at three different times. The overdensity characteristic of the vortex evolves dif- Fig. 11. Perturbed surface density at three different time (from left to right) for three planetary accretion rates, with 0.1 (top), 1 (middle) and 10 (bottom) times the nominal accretion rate, in discs with α = 10 −4 , h = 0.05. Red dots show the position of the planet. First row: density snapshots for the reduced by 10 accretion rate. As the planet is growing really slowly, the gap edges are not very steep, creating a very weak vortex. Second row: density snapshots for the fiducial accretion rate. The edges of the created gap are steep enough to trigger the RWI, creating a vortex. The vortex vanishes after ∼ 8.1 × 10 4 years. Third row: density snapshots for the enhanced by 10 accretion rate. The planet is growing so fast that the edges of the created gap are steep, triggering the formation of a vortex stronger than in the previous cases. The vortex vanishes after ∼ 5.3 × 10 4 years.
ferently whether the planet is accreting quickly or not. One can see that the created vortex at the outer edge of the gap is stronger when the planet features a higher gas accretion rate. As stated before we expect to see oscillations in the gas accretion rate, with higher oscillations for the enhanced case as the created vortex is stronger.
In top panel of figure 10 we show the accretion rates reached with the enhanced and reduced accretion rates presented in section 3.2 for α = 10 −4 . Like the fiducial case, all accretion rates show oscillations, with different amplitudes and duration in time. As expected, larger accretion rates result in stronger oscillations. The small increase that can be seen on the time evolution of the enhanced accretion rates from t 4.1 × 10 4 years to ∼ 5.1 × 10 4 years, is also due to the presence of vorticies. It corresponds to the moment when the vortices are the strongest. As they vanish after ∼ 5.3 × 10 4 years, the oscillations in the accretion rates vanish too. Our results are therefore in agreement with Hammer et al. (2017) and the vortices push material in the vicinity of the planet as their presence is linked to an enhancement in the gas accretion rate.
The presence of these vortices has a strong impact on the density profiles of the disc. The large viscous timescale at low viscosity results in a long time for the disc to adjust to the presence of the planet, resulting in an highly perturbed disc structure. The perturbed density profiles are shown in the lower panel of figure 10. Perturbations are particularly high in the inner disc.
This has a strong influence on the stellar accretion, shown on the middle panel of figure 10. Compared to what was found in section 3.3, the scaling of the stellar accretion with the planetary accretion is reversed here: a more efficiently accreting planet results in a larger stellar accretion rate. We expect this effect to flip at later time, when the inner disc will be depleted by viscosity: the viscous time at α = 10 −4 is 100 times longer than for α = 10 −2 , needing t ν = 8 Myrs for material to move from the position of the planet to the inner disc, too long for computational integration. Therefore the depletion of the inner disc is less affected by smaller viscosities: instead the inner disc is mostly influenced by the presence of the growing planet. Larger planets induce larger perturbations: the enhancement 1 created by the planet in the inner disc, at r < 1 AU in the lower panel of figure 10, is larger for the largest planet, i.e. for the enhanced by 10 case. The presence of this enhancement at the inner disc is the reason why the stellar accretion is higher compared to the disc without a planet. As soon as the disc will start to deplete its inner disc via viscous spreading, the accretion rate onto the star will behave like in figure 7, where we observe a higher decrease in stellar accretion rate with fast accreting planets. This behavior was confirmed by investigating the evolution of the stellar  Fung et al. (2014) study. None of these works include planetary gas accretion. Red dots represent our fiducial simulations. Results from the different accretion rates are represented by upward triangles for the enhanced accretion rates and downward triangles for the reduced accretion rates. The color coding is the same as for figure 4: darker colors represent larger accretion rates. The two different panels are for the two high aspect ratios: h = 0.05 (top), 0.07 (bottom). The initial mass of the planet is 20M ⊕ . Gas accretion has different impact as a function of disc parameters. This discrepancy can be explained by the time needed to the disc to react to a change of planetary mass. accretion rate for an intermediate viscosity: for α = 10 −3 (not shown), the stellar accretion rate is larger compared to a disc without a planet at the beginning of the simulation, but then becomes smaller when the inner disc starts to deplete via viscous spreading.
Viscosity has an impact on how the planet will open a gap, as can be seen on the lower panel of figure 9: depending on the viscosity, a gap is opened at different mass and the gap has a different shape. In the next section we investigate the impact of gas accretion on the gap opening mass.

Gap opening mass
Criteria for gap opening mass have been derived in several previous studies. In Crida et al. (2006), they defined the gap opening mass as being the mass at which Σ gap = 0.1 Σ unp , which we also use in our work. The derived gap opening criterion depends on the aspect ratio and viscosity: where q is the planet to star mass ratio and a planet opens a gap when P 1. This criterion has been tested for viscosities 10 −5 α 3 × 10 −1 and aspect ratios 0 h 0.3 by Crida et al. (2006). Here we focus on the gap opening mass, i.e. at which q we satisfy P = 1 as a function of the viscosity and aspect ratio. We plot this formula in the q-α space in figures 12 and 13 for the three different aspect ratios studied in this work.
In other studies, equations for Σ gap were derived. In Kanagawa et al. (2015), the depth of the gap is stated as: where K = q 2 h −5 α −1 . This criterion is valid for K 10 4 . For the mass range given on each axis, this gives a validity range from α = 2×10 −6 to α = 1.2×10 −3 for h = 0.05 and from α = 6×10 −7 to α = 5 × 10 −4 for h = 0.07, therefore valid for low viscosities here. It was derived for a uniform disk, meaning that the initial surface density was constant in radial direction. Gyeol Yun et al.
(2019) added a correction to this work by studying gap opening in a non uniform disk, having a power law initial surface density profile with an exponential cut off. They found a very similar result to Kanagawa et al. (2015), with Σ gap /Σ 0 = 1/(1 + 0.046K).
To be able to compare this corrected criterion with Crida's, we looked for which mass Σ gap /Σ 0 = 0.1. We plot it as a function of the viscosity for different aspect ratios in figure 12. Fung et al. (2014) also derived a gap opening criterion for planetary masses of the range 10 −4 q 5 × 10 −3 , for 10 −3 α 10 −1 and for 0.04 h 0.1: In figure 12, we see that Fung and Kanagawa's criteria show a linear evolution for the gap opening mass as a function of the viscosity, unlike Crida's criterion that levels off at low viscosity. This leveling off is caused by the fact that even when the viscosity vanishes, the planet is still in competition with the pressure of the disc to create a gap, creating a lower threshold for the gap opening mass. Moreover, there should be a threshold at low viscosity as the presence of vorticies generate a certain background level of turbulence, independently on the prescribed alpha viscosity.
In order to determine the influence of planetary gas accretion on the gap opening mass within our simulations, we investigated the gap opening masses of fixed mass planets in our simulations (squares on figure 12). Our planets with fixed mass seem to be opening gaps following Fung's criterion for h = 0.05 and high viscosity only. For the other aspect ratios, all gap opening criteria seem to be failing. We suspect that our simulations do not match directly any criteria as they are all derived from fits to simulations, inducing errors. Using the gap opening masses derived from the Crida et al. (2006) criterion (eq.4) we actually find deeper gaps than 0.1 Σ unp . It is also important to mention that we are investigating gap opening mass in non equilibrium discs: the mass accretion rate onto the star is evolving with time, as material is accreted onto the planet and the star. The previous studies (Crida et al. 2006;Fung et al. 2014;Kanagawa et al. 2015) were made in equilibrium discs, meaning that the stellar accretion rate is constant with time. Having a viscously evolving disc may change the gap opening mass as the accretion onto the star will create a depletion of the inner disc, influencing the material around the planet and therefore influencing how the gap is opened (especially at high viscosities). To keep consistency, we therefore compare gap opening masses within our simulations.
We can now identify three behaviors depending on the level of viscosity: -At low viscosity (for example at α = 10 −4 for h = 0.05), the gap opening mass is highly dependent on the gas accretion rate. Simulations with faster accreting planets result in larger gap opening masses, while simulations with slower accreting planets result in lower gap opening masses. -At high viscosity (for example at α = 2 × 10 −2 for h = 0.05), we find the opposite behavior: a simulation with a high accretion rate will open a gap at a smaller mass than a slowly accreting planet. -As the behavior flips between high and low viscosity, it exists a peculiar viscosity, for example at α = 10 −2 for h = 0.05, for which the gap opening mass is insensitive to the different gas In order to explain this behavior, we investigate the time needed for a fixed mass planet to open a gap as a function of its mass, at low and high viscosities. Results are shown in figures 14 and 15 for low (α = 10 −4 ) and high (α = 2 × 10 −2 ) viscosities respectively. The crosses represent the time needed for fixed mass planets of different masses to open a gap. To be consistent in the comparison, in this case the fixed mass planets are introduced exactly like the accreting planets: after the 20M ⊕ core is introduced and creates its initial gap, the mass of the planet is switched to the final mass of interest. The time needed to open a gap corresponds then to the time needed to reach the Σ/Σ unp = 0.1 threshold after the planet mass is switched to its final mass.
As expected, we see that more massive planets open gaps more rapidly than low mass planets (Lin & Papaloizou 1986). In consequence, the way the planet is introduced in the disc has a strong influence on the gap opening mass (more precisely on the time at which the 0.1 Σ unp threshold is reached). To compare with accreting planets, we show on the same figures the time evolution of the masses of the accreting planets. The dots and triangles represent the time and mass at which each accreting planet is observed to open a gap (same dots and triangles of figure 12). On the lower panels of figures 14 and 15 we show the surface density profiles of the non accreting planets at the time and mass used for the gap definition. We observe two different behaviors at low and high viscosities that can explain the gap opening masses of our accreting planets (figure 12): -At low viscosity (figure 14), the gap opening masses of the accreting planets and the one of the fixed mass planets are both inversely proportional to the planet mass (t ∝ q −X ). In this case, the gap opening phenomenon is governed by the disc response to the change of mass of the planet. Gas accretion only changes the value of the slope compared to the fixed mass planets (X −1.47 for the fixed mass planets and X −3.35 for the accreting planets). It means that the process of gas accretion is dominated by the disc response time, i.e. the time needed to the disc to react to a change of planetary mass. In this case, when the planet is slowly introduced in the disc, it will open a gap at a lower mass than a planet that is rapidly introduced. This explains the behavior seen in figures 14 and 12 at low viscosity: slowly accreting planets are slowly changing in mass, resulting in smaller gap opening masses. -At high viscosity (figure 15), the behavior is completely different. The gap opening masses of the fixed mass planets is still inversely proportional to the planet mass (t ∝ q −X ) but the gap opening mass of the accreting planets is proportional to q (t ∝ q +X ). In this case, the disc response is fast enough for planetary gas accretion to help carving a deeper gap. Gap opening is therefore dominated by gas accretion and not by the disc response time.
The key parameter here is the disc response time to a change of planetary mass. A formula for the gap opening time in an inviscid disc has been derived by Lin & Papaloizou (1986). They find that τ gap ∝ q −2 . As plotted in figures 14 and 15, we find τ gap ∝ q −1.47 for α = 10 −4 and τ gap ∝ q −12.81 for α = 2 × 10 −2 , meaning that at higher viscosity, the gap opening time is more dependent on the planetary mass. It is important to mention that the discrepancy between our result and Lin & Papaloizou (1986) comes from the importance of viscosity in the gap opening process and the definition for the gap opening time. Indeed, we define the gap opening mass as the mass needed to create a gap of a certain depth (Σ/Σ unp = 0.1) whereas they define it by equilibrating the forces applied by the planet on the disc and vice versa.
In conclusion, gas accretion has a different impact on the gap opening mass depending on the disc response time. It is a dominating phenomenon at high viscosities. These differences in gap opening mass can have an important impact on the transition to type II migration, which is investigated in the next section.

Migration maps
Planets embedded in discs interact with the gas by exchanging angular momentum, resulting in planet migration (for recent reviews see Kley & Nelson 2012;Baruteau et al. 2014). Migration is a robust phenomenon and a central ingredient in global models of planet formation that attempt to reproduce known exoplanet systems as well as our own Solar System (for a review, see Raymond et al. 2018). In order to study the type I migration of planets that could be present in the same disc as our gas accreting planets, we employ the formulae from Paardekooper et al. (2011) to estimate the torque acting on the smaller bodies. In their work, they derived a formula based on the temperature gradient T ∝ r −β T and surface density gradient Σ ∝ r −α Σ . The total torque is composed of five components coming respectively from the Lindblad torque, the barotropic and entropic horseshoe drag torques and the barotropic and entropic linear corotation torque: where γ = 1.4 is the adiabatic index and Γ 0 = (q/h) 2 Σ p r 4 p Ω 2 p is the normalizing torque depending on q the mass ratio between the planet and the star, h the aspect ratio, Σ p the surface density at the location of the planet, r p the distance star-planet and Ω p the planet's angular velocity. The total torque Γ tot is a linear combination of these torques where the coefficients F, G and K A&A proofs: manuscript no. aanda depend on the saturation parameters p ν and p χ : where p ν = 2/3 (r 2 p Ω p x 3 s )/(2πν p ) is the viscous saturation parameter depending on the horseshoe half width x s as described in Masset et al. (2006) and the viscosity ν p at the location of the planet and p χ = 3p ν /2 χ p /ν p is the thermal saturation parameter depending on the thermal diffusion χ p , which depends on the opacities in the disc, where we use the Bell & Lin (1994) opacities.
A negative torque would result in an inward migration, towards the star. On the other hand, a positive torque reflects outward migration, toward the outer parts of the disc. Inward migration is mostly due to the Lindblad component of the torque. In order to have outward migration in our locally isothermal case, where the temperature gradients are quite shallow, sharp positive radial gradients in the surface density are needed Masset et al. (2006). These gradients appear at the outer edge of the gap created by the accreting planet.
In figures 16 and 17 we show the migration maps for different viscosities and different aspect ratios when the accreting planet reaches gap opening mass. Regions of inward (blue) and outward (red) migration are represented as well as the potential equilibrium position of the low mass migrating planet: the zerotorque position (black line) and the resonances (dashed white vertical lines) with the accreting giant planet, whose position is marked by the dashed blue line.
An outer lower mass planet would stop migrating either at the zero-torque location or if it is trapped in resonance with the inner planet, depending on its migration speed (Thommes 2005;Pierens & Nelson 2008). In our simulations, if an inward migrating planet can jump the 2:1 resonance and continue to migrate further in, it would reach the position of the 3:2 resonance for the lowest viscosities, when the torques saturate. At high viscosity, for all aspect ratios, the outward migration zone is wide enough to overlap over the 3:2 resonances positions and the 2:1 resonances positions if the mass of the second planet is small enough. It these cases, the migrating planet is prevented to be locked in 3:2 resonance by the zero torque line, specifically for low mass planets. When the planet mass is small enough, the 2:1 resonance is also unreachable due to the extend of the outward migration zone. Depending on the shape of the gap created and the migration speed of the low mass planet, the capture in certain resonances might therefore be avoided. The capture in resonance is important for the Grand Tack scenario (Walsh et al. 2011), where Jupiter and Saturn get locked in 3:2 or 2:1 resonance and start migrating outward (Masset & Snellgrove 2001;Pierens & Nelson 2008;Raymond et al. 2011;Pierens et al. 2014;Chametla et al. 2020). We discuss in section 5.3 the impact of resonance capture on the structure of our solar system. Pierens et al. (2014) investigate which disc parameters (Σ 0 and α) and capture in resonance are needed to allow outward migration of Jupiter and Saturn. For h=0.05 (middle row in figure  16) and our surface density (Σ/Σ MMSN = 1.5) 2 with α = 10 −3 , outward migration for the Jupiter/Saturn pair would happen in the 3:2 resonance. Still according to Pierens et al. (2014), the Jupiter/Saturn system would show divergent migration for α = 10 −2 and be captured in the stable 2:1 configuration for α = 10 −4 . In figure 16, we see that the capture in the 3:2 resonance is possible for α = 10 −3 and h = 0.05 as for these disc parameters, Saturn would not open a deep gap (m S = 95M ⊕ < m gap = 117M ⊕ ). Therefore our migration maps agree with the results found in (Pierens et al. 2014), assuming that Saturn's migration speed allows it to be locked in the 3:2 resonance. In addition our results show that N-body simulations of growing planets should take the disc profile into account to assess the migration behavior of the planets correctly.
As a planet is growing and carving a deep gap, it is transitioning from type I to type II migration. In figure 18, we show the migration maps for h = 0.05 and α = 10 −4 for the fiducial, enhanced and reduced by 10 planetary gas accretion rates. If the planet accretes slowly, a gap is opened at a smaller mass than when the planet accretes fast for these disc parameters (see section 4.3). The mass at which the planet will switch from type I to type II migration (depicted by the gray area in figures 16, 17 and 18) is therefore dependent on the planetary accretion rate, assuming that the second potential planet in the system accretes at the same rate and would thus have the same gap opening mass. At these disc parameters, the difference in gap opening mass does not have an impact on the capture in resonance: the region of outward migration is small and does not overlap with resonances for planets more massive than 15 M ⊕ .
The dynamics of multiple planetary systems can be highly impacted by planetary gas accretion, via the influence on the migration type and on the potential trapping in resonance.
2 Note that their surface density profile is different from our setup: Σ ∝ r −3/2 as opposed to Σ ∝ r −1 in this work. Therefore the outcome of planet migration in our case might be different from their results. Fig. 18. Migration maps for α = 10 −4 and h = 0.05 for low mass planets orbiting near our accreting planets. The accreting planet reached gap opening mass. The three panels represent different accretion rates: the reduced by 10 (left), fiducial (middle) and enhanced by 10 (right). The plotted information are the same as in figure 16. As gas accretion has an influence on the gap opening mass, the migration maps look different due to the switch from type I to type II migration when a gap is opened. Here the Σ/Σ unp = 0.5 is at the same mass for all the different accretion rates as the initial depth for these discs parameters is smaller than 0.5 (see table 2).

Accretion onto the star
In this work, we studied the influence of planetary gas accretion on the stellar accretion (section 3.3 and 4.2). We show that the stellar accretion is reduced with increasing planetary accretion rates. Even though at low viscosity our results show an enhancement of the stellar gas accretion rate, we expect this trend to flip and follow the high viscosity case after reaching the viscous time needed for material to reach the inner disc from the planet position.
In both previous cases, the stellar accretion was only decreased or increased by a factor of up to three compared to the case of a disc without any planet. These results are quite different from what Manara et al. (2019) find. In their models, they find that stellar gas accretion rate can be reduced by over two orders of magnitude due to the presence of giant planets. In our simulations these large spreads of stellar accretions could only be achieved by changing the disc viscosity over orders of magnitude, however, this parameter is fixed in the simulations used in Manara et al. (2019). The difference might come from the difference of gas flow between a 2D and a 1D disc: indeed, in Manara et al. (2019) approach, the gas accretion onto the star is derived from the viscous spreading of a 1D disc containing a giant planet. As Lubow & D'Angelo (2006) concludes in their work, in 2D, mass can flow through the gap even in the presence of a giant planet. Moreover, Manara et al. (2019) planetary gas accretion rates are higher than ours as it is proportional to the unperturbed surface density in their accretion routine (Mordasini et al. 2012). Therefore the 1D model from Manara et al. (2019) might be over estimating the efficiency of the planet to block material from the inner disc.
In addition, Lubow & D'Angelo (2006) showed that the modelling of the inner disc and how gas is accreted onto the planet is of crucial importance to calculate the stellar accretion rate. Our simulations indicate that planetary gas accretion might have a smaller impact than expected on the stellar gas accretion rates.

Implications for observations
Observations of protoplanetary discs revealed structures in the dust distribution (Pinilla et al. 2012;Andrews et al. 2018;Zhang et al. 2018). However, these dust structures are ultimately caused by the drifting dust grains trapped at pressure bumps in the gas discs. Therefore as mentioned in previous sections, looking at the pressure bumps created by planets in the gas can give an insight on how the dust will be trapped and how gas accreting planets might influence the interpretation of the observations. When we compared the pressure bumps in sections 3.1 and 3.2 for our simulations with the fiducial parameters, we found that the pressure bumps were too similar to show distinguishable features. This result is expected, as in section 4.3 we showed that for these disc parameters gas accretion has almost no influence on the gap depth. On the other hand, at low viscosity the gap opening mass changed dramatically with the planetary accretion rate.
In figure 19 we show the surface density profiles of discs containing accreting planets at gap opening mass and the corresponding pressure bumps for α = 10 −4 and h = 0.05 and for the different planetary accretion rates studied. Profiles plotted in the top panel show similar gap width for the same gap depth. The created pressure bumps are therefore very similar to each other. Even if the difference in planetary mass is of almost a factor 2, the created gaps are too similar to show features that can be distinguished by observations. Our simulations thus indicate that for the same gap depth, gas accretion has a negligible impact on the gap width. However it has an impact on the mass creating such a gap: depending on the gas accretion rate, a very similar gap can be created by planets of different masses as we can see on figure 19. This implies that the planetary masses derived from 2D simulations designed to match the observations might be off by up to a factor 2 for α = 10 −4 and h = 0.05. This factor is dependent on the viscosity ν and is larger for lower viscosities, as one can see in figure 13 where h = 0.03.
Recent studies (Teague et al. 2018;Pinte et al. 2020) show that CO velocity field observations can be used to observe planets in discs. As gas accretion has an impact on how the gap is formed, it has an impact on the velocity field of the gas around the planet. In Teague et al. (2018), hydro-simulations with non accreting planets were run in 2D. The planetary masses were derived by comparing the deviation to the rotational velocity in the simulations and the observations. Their derived masses can be altered by gas accretion if the disc parameters make the gap formation sensitive to gas accretion.

Implications for the Grand Tack scenario
The gap opening mass is an important parameter for the structure of planetary systems. When multiple planets are taken into account, like in the Grand Tack scenario with Jupiter and Saturn (Walsh et al. 2011) or in recent N-body planet formation simulations , the time scales for gas accretion are important when investigating the growth and migration of multiple planets in the same disc, as required by the Grand Tack scenario. Our study here supports the fact that gas accretion is an important factor not only for planetary growth, but also for the migration behavior of other planets and the gap shape and should be taken into account in future simulations and interpretation of observations of protoplanetary discs where planets are suspected.