Adapting a gas accretion scenario for migrating planets in FARGO3D
Luiz Alberto DePaula Tatiana Alexandrovna Michtchenko


firstpage–lastpage
Computing the Inversion-Indel Distance
    Eyla Willing, Jens Stoye and Marília D. V. Braga
E. Willing, J. Stoye and M. D. V. Braga are with the Genome Informatics group, Faculty of Technology and Center for Biotechnology (CeBiTec), Bielefeld University, Germany
    December 30, 2023
===============================================================================================================================================================================================================================













FARGO3D has been proposed to investigate numerically the gravitational interaction of the planet with the gas disc, 
providing original outcomes. However, FARGO3D does not consider the gas accretion of the planet that may affect 
the planetary migration process. Thus, the aim of this paper is to include a model for the gas accretion in the 
FARGO3D code. For this task, we choose Russell's scheme, which is an adaptation of Kley's model for the gas 
accretion onto migrating planets. Initially, we study the influence of the gas accretion 
onto the planet on type II migration of giant planets.
For this purpose, we follow the evolution of the planets in a two-dimensional locally isothermal 
disc with a specific accretion rate and different values of viscosity and planetary mass considering 
two scenarios. In the first, the gas mass is withdrawn from the disc and is not added to the planet. In the 
second the planets migrate, while their masses grow due to the gas accretion.
Then, we extend our  study to estimate the time needed for a low mass planet to open a gap in the gas disc, and compare its 
value with the characteristic time of type I migration.





Protoplanetary discs – Planet disc interaction – Planet formation








§ INTRODUCTION

   According to the classical theory of planetary formation, giant planets form in outer regions of the protoplanetary disc where
   there is enough material available <cit.>.    However,  a large number of giant planets with very short period orbits,    known as “hot Jupiters”, has been discovered <cit.>. One of the mechanisms,  proposed to explain the locations   of such planets, is the planetary migration resulting from the gravitational interaction of the planet with the    gas disc <cit.>.    Indeed, this interaction leads to an exchange of the angular momentum between the disc and the planets, changing their orbital parameters and resulting, in general, in the reduction of the planet's semiaxes and damped eccentricities <cit.>.

   Several numerical codes have been employed to investigate the gravitational interaction of the    planet with the gas disc, such as NIRVANA <cit.> and FARGO <cit.>, among others.    However, the majority of the codes does not consider gas accretion onto the planet or deal with    simplified models. This is due to the fact that    the accretion process depends on several aspects of the disc structure and on a variety of physical    processes and their consideration  increases significantly computational costs.    Therefore, little is known about how the accretion process affects the local structure around the    planet, and, consequently, how it affects the planetary migration.
   Moreover, according to the classical theory of planetary formation, the rate of the gas accretion onto    the planet also depends on the mass of the planet. In fact, for low mass planets, the gas accretion process    is regulated by the Kelvin-Helmholtz contraction of the planet's gas envelope. On the other hand, as    the mass increases, the hydrostatic equilibrium breaks up and, as a consequence, the planet enters in a runaway
   accretion regime, at which the planet's mass increases very fast.

There are two main migration regimes known as type I and II migration.    Type I migration occurs for low-mass planets (≲ 50M_Earth),    whose interaction with the disc is weak enough to leave the disc structure almost unperturbed. This type of migration can be approximated analytically by the linear regime <cit.>, when the total torque    on the planet due to the gas is defined by the combined effects of the Lindblad and corotation resonances <cit.>.    <cit.> were the first to formulate
   the disc-planet interaction and derive formulae for the migration process. Later, <cit.>
   extended the linear regime of excitation to three-dimensional isothermal discs and recalculated    the torques of the Lindblad and corotation resonances.

   On the other hand, type II migration occurs for high-mass planets (> 0.5 M_J).    In this regime, a planet opens an annular gap, within which the disc surface density  is reduced  substantially from its unperturbed value; as a consequence, the gas in this region
   contributes less to the torque on the planet, thus reducing the amount of total torque of the system <cit.>.
   The gap is responsible for a coupling between the orbital migration of the planet and
   the viscous evolution of the protoplanetary disc. So, in the case when the gas does not cross the gap,
   the planet migrates inward in a timescale determined by the disc's viscosity <cit.>.
   However, studies show that the gas can cross the gap via horseshoe orbits, implying that the rate of migrates are different than those of the viscous evolution of the disc <cit.>.

   The calculation of the migration rate of the planets of intermediate mass is complicated because this process cannot be fitted by neither the  type I nor the type II migration, thus this issue still remain open <cit.>.

   Numerical simulations allow us to investigate how the accretion of gas onto the planet affects the    flow of gas in the gap, modifying the planet's migration rate. For example, <cit.>,    using the method originally described by <cit.> for gas accretion of the planet, have found that the orbital migration of massive planets is always inwards and that the planet reaches the central star after 10^4 initial orbital
   periods, regardless of the fact that the planet is accreting mass or not.    Later, <cit.>, adding some refinements for the gas accretion to the same model, have   shown that a fast migrating planet is able to accrete more gas than a slowly migrating planet.
   Moreover, they showed that the amount of gas crossing the gap changes during the migration process    as the migration slows down.

   This work investigates type II migration considering the model for gas accretion presented in    <cit.> and adapting it to suit the FARGO3D code <cit.>. Indeed, the FARGO3D code does not include a model for gas accretion, although, this process can affect the planetary migration.   Following <cit.>, we use a timescale of gas accretion that considers    the orbital period of the gas around the planet. Additionally, we 	 using    the minimum value between the Bondi radius and the
   Hill radius. Finally, we introduce a gas accretion rate to deal with small planets, using the Kelvin-Helmholtz contraction time scale for planet's gas
   envelope. In particular, we investigate the influence of Russell's gas accretion model on the type II migration,    and the time for a low mass planet to open a significant gap in the disc.

   This manuscript is organized as follows. In Sect. 2, we describe the gas accretion model and introduce  the numerical model for the disc in Sect. 3.    In Sect. 4, we discuss how the gas accretion model affects the migration of giant planets and compare    our results with those reported in <cit.> and <cit.>. In Sect. 5, we investigate the time it takes    for a accreting planet with initial mass ≃ 20M_Earth, on a fixed orbit, to open a gap in the gas disc. Finally, we discuss the results in Sect. 6.





























§ GAS ACCRETION SCENARIO


<cit.> was one of the first works to introduce a model for gas accretion to a planet forming in the disc.  <cit.>  modified the model including Bondi's radius for gas accretion zone, effects of the disc height, Kelvin-Helmholtz gravitational contraction time scale, among others.
The calculation of the Bondi's radius is  a condition for the calculation of the gas accretion zone.
Bondi's radius determines the region around the massive body immersed in the gas cloud, from which the body captures the gas. By definition, in this region, the resulting (thermal + relative to the body) velocity of the gas  is smaller than the escape velocity of the gas <cit.>. This condition provides that


    R_B_ij = 2GM_p/c_s^2 + Δ v_ij^2,


where Δ v_ij is the module of the relative velocity between the planet and gas residing in cell (i,j),
c_s is speed of sound in the cell, G is the gravitational constant and M_p is the mass of planet.

The radius of the accretion zone of  the planet, R_a_ij, is then determined by the condition


    R_a_ij = min (R_B_ij, R_H),


where R_H is the Hill radius of the planet given by


    R_H = a ( M_p/3 M_⋆)^1/2,


with a and M_⋆ denoting the semimajor axis of the planet and the mass of the star, respectively. Defined in this way, R_a_ij determines a minimum distance at which a cell should be from the planet, to be considered within the accretion region of the planet. If the distance from the cell to the planet r_ij is smaller than R_a_ij, then the gas will be removed from this cell during accretion.

The mass of gas removed from a cell with position [i,j] inside the accretion zone in a particular time step dt (intrinsic time of the hydrodynamic code), is calculated as


    δ m_ij = f_ijΣ_ij A_ij dt,


where Σ_ij is the cell surface density, A_ij is the cell surface area and f_ij is the frequency
with which the gas is removed from the cell, defined by


    f_ij = 1/2π√(G M_p/r_ij^3),


where r_ij is the distance between the cell and the planet. Thus, the total mass, Δ m_n, removed from all accreting cells per nth time-step , is a sum


    Δ m_n = ∑_i,jδ m_ij,

and the mass of planet is then increased by this amount.

The linear momentum is also transferred from the gas to the planet during accretion. However, according to results obtained by <cit.>, this process does not change significantly  the migration rate of the planet. So, in this paper, we choose to work only with the increase in planetary mass.

The application of the Bondi's radius to calculate the accretion zone results in a dynamic accretion rate, which varies around the planet. 
The model is adaptable for the runaway regime for gas accretion of giant planets. On the other hand, low mass planets are outside the runaway regime and, 
as shown in <cit.>, do not accrete mass in a rate beyond the Kelvin-Helmholtz contraction time scale. Indeed, the Kelvin-Helmholtz gravitational 
contraction time scale, τ_KH, represents the characteristic growth time of a planet's gas envelope. <cit.> found that 
τ_KH is strongly dependent on the planetary mass and is expressed approximately as


    τ_KH≃ 10^b( M_p/M_⋆)^-c( κ/1 cm^2g^-1) yr,


where κ is the grain opacity with power index b = 8 and c = 2.5.
In the present work, however, we adopt an approach used in <cit.>, setting κ = 1.0 cm^2g^-1, which gives


    τ_KH≃ 10^8( M_p/M_⋆)^-2.5yr.

This expression gives an approximate growth time of the gas envelope for a planet with low mass.

According to <cit.>, the rate which limits gas accretion onto the planet is determined by the minimum
between the potential gas accretion rate (Ṁ_U) and the Kelvin-Helmholtz gas accretion
rate (Ṁ_KH) for a particular time step dt (intrinsic time of the hydrodynamic code):


    Ṁ = min (Ṁ_U, Ṁ_KH) = min(Δ m_n/dt,M_p/τ_KH) ,


where τ_KH is given by Eq. <ref> and Δ m_n by the Eq. <ref>.

The accretion rate given by Eq. <ref> determines the ratio of gas removed from the cells
within the planet's accretion zone. The ratio of gas accretion, μ, is:


    μ = Ṁ_KH/Ṁ_U,


where μ⩽ 1. The mass removed from each cell in the planet's accretion zone is reduced
by μ. Thus,


    δ m_ij = μ f_ijΣ_ij A_ij.


This equation implies that, if μ < 1, Ṁ_KH < Ṁ_U, and,
consequently, the Kelvin-Helmholtz rate limits the gas accretion. Otherwise, μ = 1 and the
accretion of gas occurs in the runaway regime.


In summary, unlike the scenario described in <cit.>, the model described by Russel allows us to work with 
the accretion mechanism for small planetary masses using the Kelvin-Helmholtz contraction time scale for planet's gas envelope. 
Besides that, it is possible to adapt the size of the accretion zone taking into account the radius of Bondi. In addition, 
Russell's model uses, for the runaway regime, a timescale of gas accretion that considers the orbital period of the gas around 
the planet that results in a dynamic range of accretion time scales for each cell where gas is accreted.

In this work, Russel's model described above has been adapted in FARGO3D, for use in GPU programming.



§ NUMERICAL SETUP


   To reduce high computational costs, we apply the hydrodynamic code FARGO3D <cit.> in the 2D regime in our simulations. The advantage of using FARGO3D in 2D setup, instead of using its predecessor FARGO, is that FARGO3D allows GPU computation, thus reducing computational time. The application of the three-dimensional model to gas accretion study is left for future works.

   In order to compare results, we choose setup parameters similar to those adopted for the NIRVANA code in <cit.>.    Assuming a thin gas disc, with typical vertical height (H) much smaller than its radius, we introduce a 2D approximation    of the problem.    The disc is considered to be locally isothermal, with viscosity defined by the parameter α from <cit.>.    The star, with M_⋆ = 1M_⊙, is located at the origin of a cylindrical coordinate system, while the disc extends from 1.56 to 15.6 au (0.3 to 3.0 in code units). 
   The resolution of the disc is fixed at 582 × 1346 cells equally spaced and no significant differences (< 5%) in the final position of the planet and the migration    time have been found using higher resolutions. The planet of mass M_p is placed initially
   on a circular orbit at a distance of ⋍5.2 ua, that corresponds to r = 1.0 code units.

   Following <cit.>, we work with a steady accretion flow through the disc, that is, with the constant local    accretion rate ṁ. In this situation, the radial and angular components of the velocity of gas, v_r and v_θ, respectively,   are given by <cit.>:


    v_r = - 3/2α h^2 r Ω_K,
    
        v_θ = √(1 - 3h^2/2) r Ω_K,


where Ω_K = √(GM_⋆/r^3) is the Keplerian orbital frequency of the gas at the position r and the parameter h corresponds to the relative disc thickness and has constant value h = H/r = 0.05. The surface density of the disc is given by
<cit.>:


    Σ = ṁ/3πα h^2√(()G M_⋆) r^-1/2 = Σ_0 r^-1/2,


where Σ_0 is the superficial density at r = 1 code units. These equations can 
be obtained for stationary accretion disks with constant accretion rate.

Table <ref> shows the parameter space used in our simulations. Throughout the paper we will
refer to the disc defined by the of parameters (α, q, ṁ) = (0.003, 0.001, 10^-7 M_⋆/year) as the standard disc, with q = M_p/M_⋆.
It is worth mentioning that in the FARGO3D code the initial profile of the surface 
density of the disc is obtained through the free parameters for the disk aspect ratio (h = 0.05), 
the flare of the disc (γ = 0) and the surface density in r = 1 (Σ_0).
Table <ref> displays the values of Σ_0 obtained for each parameter set
(α,ṁ) according to Eq. <ref>. The values are in kg·m^-2.






We modify the FARGO3D code in such a way that the torque on the planet due to the gas is reduced
when the gas is inside a region with radius equal to 0.8R_H around the planet.
Following <cit.>, we use a Fermi-like tapering function the  given by


    f_taper(r_cell) = 1/1 + exp(-r_cell - 0.8R_H/0.8R_H),


where r_cell is the distance of the cell to the planet and R_H is the planet's Hill radius.

For comparison purpose, we choose the same boundary conditions used by <cit.>. So, we disable the damping of the surface density in the FARGO3D code, to avoid accumulation of gas at the inner boundary and allow gas to flow freely as desired.



§ GIANT PLANET MIGRATION


In order to investigate how Russell's accretion model <cit.> affects the type II migration, we analyse two scenarios for the migration of giant planets. In the first scenario (Sec. <ref>), the planets do not grow during migration, that is, the gas mass is withdrawn from the disc  and is not added to the planet. In the second scenario (Sec. <ref>), the planets migrate, while their masses grow due the gas accretion.

First, we need to know how the gas accretion of the planet affects the profile of the gap, also as the flow of gas through the gap opened by the planet modifying the torque on the planet. Additionally, we analyse the rate of planetary migration for different values of the disc viscosity. In this
work we have three conditions for the viscosity of the disc, we use the term low viscosity when α = 0.001, medium viscosity
when α = 0.003 and high viscosity when α = 0.01.

These analyses require the gas disc to be in equilibrium. Thus, following the approach of <cit.>, we calculate the time needed for the disc to reach the steady state. Figure <ref> shows the profiles of ṁ (r) as a function of the radius of the disc after  different time periods, calculated for the planet of 1.0M_J fixed at r = 1 code units. To calculate the local accretion rate, we include a specific subroutine in the FARGO3D code, which  integrates the azimuthal mass fluxes produced by the code.



As seen in Figure <ref>, although the profile of the accretion rate does not attain a constant
value, it converges to a nearly flat function between 5000 and 5800 orbital periods of the disc.
Thus, we can assume that the disc reaches equilibrium state after 5000 orbits and adopt this value
as a standard timespan.




Figure <ref> shows the azimuthally averaged surface density for non-accreting (solid lines) and accreting (dashed lines) planets, obtained after 5000 orbital periods of the disc,  for the three values of planetary masses and three different viscosity regimes used in our simulations. As expected, we observe the completely developed gaps in all cases. The gaps are increasing in depth with the increasing planet mass and the decreasing viscosity of the disc. The regions of overdensity at the positions of the planets (r = 1 code units) appear in the case when accretion processes are not considered (solid lines). This is due to the accumulation of gas around the non-accreting planet, treated as a point in the hydrodynamic mesh; when the gas accretion into the planet is taken into account, this artefact disappears (dashed lines).  Finally, we have performed the same simulations for different accretion rates of the disc, (10^-8M_⋆yr^-1 and 10^-9M_⋆yr^-1), but we have not observed any significant difference; the same behaviour is reported for the NIRVANA code in <cit.>.




 §.§ Migrating planets without growth


Soon after the gas disc reaches equilibrium, the planets are allowed to migrate. We followed the migration
tracks of planets for the parameter sets shown in Table <ref>. The measurements were performed for accreting and non-accreting planets during 1200 orbital periods of the disc; in the accreting case, the mass withdrawn from the disc was not added to the planet. 



The migration tracks obtained are shown Fig. <ref>, where the parameters of the standard set are changed one-by-one, in order to analyse the behaviour of the planet.  We show the dependence of the rate of migration on the viscosity of the disc in the left frame of Fig. <ref>. In the case of high viscosity (α = 0.015), there is a notable difference between the behaviour of accreting (dashed line) and non-accreting (solid line) planets. With decreasing α, the difference also decreases and we note no difference in the case of low viscosity (α = 0.001). 

The middle frame of Fig. <ref> shows that the behaviour of the planet for different accretion rate of disc is only slightly different for accreting (dashed line) and non-accreting (solid lines) planets.

According to the right frame of Fig. <ref>, the effect of accretion process is much more significant for low-mass planets.  This is due to the fact that it accelerates the opening of the gap in a region where the gravitational effect of the planet is less intense.

As shown in Fig. <ref>, in general, the orbital decay of the planet is slower for accreting planets than for non-accreting planets. This difference is less significant with the increase of the planetary
mass (right frame) and for the decreasing accretion rate of disc (middle frame). This is due to the fact that the density along the disc decreases with the decrease of the rate of accretion (see Eq. <ref>); as a consequence, the total mass of the disc also decreases. Moreover, the gas mass in the nearby region of the planet also decreases when planets are accreting. In fact, this decrease in the gas mass near the planet causes the effect of the torque on the planet to decrease, thus, large planets in low mass discs have a smaller effect on their migratory process due to gas accretion. Similar results were obtained by <cit.> using a different model for gas accretion.

It should be noted that some simulations were excluded from our analysis. This is the case when the planets escape from the gap, what occurs for the non-accreting planets with the masses of 0.2 M_J and 0.5 M_J, in the low (α = 0.001) and medium (α = 0.003) viscosity regime. This also occurs for the planet with the mass of 0.2 M_J for accreting planets in all viscosity regimes considered.

In what follows, we analyse how the gas accretion model influences both the gap profile and the mass fluxes in the gap.



  §.§.§ Gap profiles for accreting and non-accreting planets


The gap created by the planet in the disc has a significant impact on the migration process because the decreasing superficial density of gas in the gap region reduces the contribution to the torque on the planet <cit.>.

Figure <ref> shows the time evolution of the gaps created by non-accreting planets (left panel) and  accreting planets (right panel). We compare the plots when the planet is in the same radial position, consequently, time is different in both plots. Looking at the planet at given positions helps to better understand the difference in local physical processes. 
For accreting planets, we note a larger lost of mass at the inner border of the disc due to the boundary conditions. When the planet is far away from the inner border of the disc, this effect does not influence on the profile of the gap.
A similar analysis was done for different parameter sets from Table 1 and the conclusion made is that effect is stronger in the higher viscosity regime and for lower accretion rate of disc.



Figure <ref> shows the local disc structure in the gap region obtained for the standard parameter set. The gap profile suffers no significant changes during the migration process of accreting and non-accreting planets.
In addition, in Fig. <ref>, we observe that the density bump (mass concentration around the position of the planet) is eliminated in the case of accreting planets. This feature remains during planetary migration and the same effect can be observed  for different values of the planetary mass and accretion rate of the disc.



Considering different values of the disc mass (i.e., accretion rate of the disc), we have found that, for fixed α-value, the gap structure near the planet does not change considerably, in both cases of accreting and non-accreting planets. These results agree with the results obtained with the NIRVANA code in <cit.>, where only non-accreting planets were investigated.

In summary, for accreting planets, our results show that, despite the accretion of gas generating deeper gaps,
the gas accretion does not alter significantly the gap profile around the migrating planet provided that
the gap is established.



  §.§.§ Flow across the gap for accreting and non-accreting planets


Figure <ref> shows the averaged azimuthal torques, which act on non-accreting (left panel) and accreting (right panel) planets, in their close vicinity. The torques were calculated for the standard parameter set. The magnitude of each torque is normalized with respect to the maximum value obtained along the radial distance. In both cases, the torques are stabilised in a few hundred orbits and reach the highest values in the vicinity of the planet (from -5 to +5 in the x–axis in Fig. <ref>). For both accreting and non-accreting planets, the plots shows that the torques  depend only slightly on the position of the planet in the disc.

Additionally, Fig. <ref> shows that the torque produced by the internal disc is proportionally stronger (of the order of 10%) for accreting planets than for non-accreting planets. The similar behaviour was observed in the low viscosity regime. In the case of the high viscosity regime, we have noticed a positive torque which was closer to the accreting planets. The gray area in Fig. <ref> shows the region of attenuation defined by Eq. <ref>. For accreting planets, the torque within this region is smaller that is related to the decreasing mass due to the gas accretion by the planet. 



In Fig. <ref>, we show the mass accretion rate in the neighbourhood of the planet, ṁ(r), at different positions of the migrating planet,  for non-accreting (left frame) and accreting (right frame) planets. For non-accreting planets (left panel), there is a large positive flow (contrary to the movement of the planet), which arises due to the fact that the migration velocity of the planet is higher than the viscous accretion speed. The same behaviour was observed in the case of low viscosity regime. On the other hand, in the high viscosity regime, with low accretion rates for the disk,  we have found some cases for which the gas flow is negative in the region close to the planetary gap. In this cases the migration velocity of the planet is lower that the viscous speed. These results are very similar to those obtained by <cit.> for NIRVANA code, for non-accreting planets.



For accreting planets (Fig. <ref> right), we note a sharp fall of the gas flow in the accretion region. This is due to the fact that the planet is treated as a massive point in the hydrodynamic mesh, despite that the mass is withdrawn from the disc in the accretion region around the planet.
Analysing the mass flow around the planet in the case of high viscosity regime, we have noted that the fall of the gas flow is sharper due to the fact that the gas mass near the planet is smaller in this case. The observed asymmetry of the gas flow in the case of the accreting planet results in a distribution of gas around the planet, which is different from that around the non-accreting planet. As a consequence, the differential torque is smaller for accreting planets when compared to that for non-accreting planets.

Finally, in the case of low accretion rate of the disc and the high viscosity regime, the positive torque due to the internal disc is proportionally higher in regions close to the planet than in the case of the medium viscosity regime. However, the effect of attenuation of torques on the planet given by Eq. <ref> does not allow the reverse migration of accreting planets.



 §.§ Migrating planets with mass growth


In this section we analyze the growth of a planet with the initial mass of ≃ 1M_J. Since, in this case, we do not need to wait for the disc to attain the equilibrium state,  we allow the planet to start migrating after only 300 orbits, which is the time necessary to the disc stabilize in the planet's presence. In the first 300 orbits, the planet is fixed on circular orbit,  the gas is removed from the disc and the mass of the planet is kept constant.

To understand the influence of the planetary mass increase on the planet's migration process, we compare different scenarios of gas accretion onto the planet. In the first scenario, described in <cit.>, the accretion rate is adjusted by the parameter k_a, which is a factor of increase or decrease of accretion. The second scenario corresponds to Russel's model <cit.> described in Sect. <ref>.
Two cases are observed: in case R, the mass is removed from the disc, but does not added to the planet, while in case RA the mass is removed and subsequently added to the planet. The linear moment of the accreted gas is not added to the planet, since its change does not interfere significantly  on the resulting process <cit.>. Finally, for the sake of comparison, we include a model  N, which does not consider for gas accretion, thus no mass is removed from the disc by the planet.




The time evolution of the semimajor axis of the planet with initial mass equal to ≃ 1M_J, in the different scenarios, is shown in Fig. <ref>. Comparing the N model (black curve) with the R (blue curve) and RA (cyan curve) models, we can observe that the removal of mass from the region around the planet decreases significantly the migration rate. Moreover, when the removed mass is added to the planet, the planetary inertia increases and, consequently, the the migration rate decreases. The decay of the semimajor axis in the model RA follows closely Kley's model with larger value of the parameter k_a = 2.0. According to <cit.>, this value is a saturated upper limit of the accretion rate.



The growth of the planetary mass initially with ≃ 1M_J is shown in Fig. <ref>, where the continuous cyan curve corresponds to Russell's model RA, while the dashed curves correspond to Kley's model, with different values of the parameter k_a. We observe that, according to the model RA,  the planet reaches the mass of more than 3.5M_J in 1000 orbits, which is very similar to what is observed using Kley's model with k_a = 2.0. This shows that, indeed, the runaway regime for gas accretion in <cit.> is equivalent to the model in <cit.> for high accretion rates.



The surface density of the gas around the planet when the planet crosses r = 0.85R_0, is shown in Fig. <ref>, for the same models presented in Fig. <ref>. As expected, the gap is deeper when the planetary mass is higher at that instant, that happens for higher gas accretion rates. According to Fig. <ref>, the larger values of the mass correspond to the scenarios RA and k_a=2.0. Also, the higher accretion rates are responsible for the decreasing of the overdensity effect in the vicinity of the planet. Once again, we notice that the result obtained for the runaway regime with model RA is very similar to that obtained with k_a = 2.0.

Fig. <ref> shows gas accretion rates for a growing giant planet. The behavior for the gas accretion rate in runaway regime is qualitatively very similar to the one obtained by <cit.>, that uses a more complex model for gas accretion. 
However, in our model the planet can migrate and the disc model is different, so the results can not be compared in quantitative terms. A more detailed analysis that includes an adiabatic disk will be carried out in future works.



The results obtained show that Russell's model <cit.> and Kley's model <cit.> provide similar results when work with the giant planets in the runaway regime for gas accretion. However, the model in <cit.> is not adapted to simulate the accretion of the low mass planets, because it does not account for  the Kelvin-Helmholtz scale. In contrast, Russell's model can be used to investigate the growth of the low mass planets since it adopts this approach.




§ GROWTH OF LOW MASS PLANETS


In order to study the accretion process in the case of low mass planets, we place a planet of mass ≃ 20M_Earth in a fixed orbit (r = 1 code units) and analyze the increase of the planetary mass using the model for gas accretion described in section <ref>.



Figure <ref> shows the in situ growth of the planet in the disc described by the standard parameter set from Table 1.  We observe that initially the planetary mass increases monotonously and slowly, but this process accelerates significantly when the planet enters in the runaway regime. This occurs  after approximately 1800 orbits, when the planetary mass reaches ∼ 0.2M_J, this behaviour is in agreement with other models for gas accretion <cit.>. It should be noted that any model which considers only runaway accretion regime, would not be able to reproduce this behaviour, since, in this case, the formation time would be  very short. The runaway accretion decreases the time for the planet to establish type II migration, since the gas along the planetary orbit is quickly exhausted.

Fig. <ref> shows gas accretion rates for the planet fixed in r = 1 code units. 
The behavior for the gas accretion rate is similar to the one obtained by <cit.>, 
that uses a more complex model for gas accretion. The only difference observed is the scale of the accretion rate, which is 
due to the fact that our model for the gas disk has a much higher surface density, thus resulting in higher accretion rates. 
In <cit.> it is also possible to find a complete scenario that 
includes the formation of the solid core. Here, we already start with a mass of ≃ 20M_Earth. 
A complete scenario including a solid accretion model within 
the FARGO3D is under development and will be analyzed in future work.



We analyze the effect of gas accretion on the movement of the gas around the planet by performing a streamline analysis, shown in Fig. <ref>. In Fig. <ref>a, we observe that, at the beginning of the accretion process, the accreted gas comes from regions near the planet. Moreover, we notice a high gas density region around the planet. In Fig. <ref>b, the gravitational perturbations due to the planet tend to clear the gas along the planetary orbit, at the same time increasing the gas density in the region close to the planet. The accretion rate of gas is still too slow to decrease the gas density close to the planet. As the planetary mass increases, the gas mass tends to accumulate in the spiral arms. The mass transfer from the outer disc to the inner disc in the region close to the planet is still small.



In Fig. <ref>c, the planet starts to accrete the gas in the runaway regime; we see the mass depletion in the region around the planet, and the gas density around the planet tends to decrease. Following, the planet begins to accrete gas from the connection of the inner disc with the outer disc. In Fig. <ref>d, we see a stronger mass transfer between the inner and outer disc. The gas tends to surround the planet and the high rate of gas accretion tends to clear the nearby region of the planet. Finally, in Fig. <ref>e, an amount of gas surrounds the planet, where it is accreted to the planets, while the rest of the gas is accelerated by the planet and driven to the inner disc. Thus, we can observe an intense mass flow between the inner and outer discs.

In summary, much of the mass during the early stages of accretion comes from regions along the orbit of the planet. However, with the increase of the planetary mass and the clearing of the orbital region, a connection between the inner and outer discs is established, from where the mass is accreted now. In any case, despite the accretion of gas, the time needed for the planet to clear its orbit is still very long. Thus, type I migration would be sufficient to drive the planets close to the star until the planet enters a type II migratory process.

In the model investigated in the paper, the gas accretion process for low-mass planets does not depend substantially on the position of the planet in the disc. It should be stressed that we do not consider solid core growth due to accretion of planetesimals that is a process which is strongly dependent on the density of solids along the disc. Therefore, it is necessary to consider accretion of solids to increase the rate of mass accretion. In addition, the accretion rate of solids is important in determining the critical mass for which gas accretion begins <cit.>. A model for accretion of solids will be elaborated in future work.



§ CONCLUSIONS


In this work we studied the migration of massive planets in locally isothermal disc in equilibrium using the FARGO3D code. Initially, the position of the planet was fixed and the structure of the gap was analyzed 
for different values of the planetary mass and disc viscosity (Fig. <ref>). We obtained that the depth of the gap increases with planet mass and reduces with increasing viscosity. This result is in agreement 
with that obtained by <cit.> using the NIRVANA code for non-accreting planets. For accreting planets, the gaps are deeper, while the overdensity regions close to the planets are vanishing. 
The results above are independent of the initial accretion rate of the disc.

Analyzing the flux through the gap during the migration process, we observe rapid inward migration in the case when the planet decays faster than the typical viscous inward drift of the disc material. We find that, 
for smaller disc mass and higher viscosity, the planet decays more slowly than the typical viscous inward drift of the disc material. 

In this paper, we performed simulations for both non-accreting planets and accreting planets. In latter case, we adapt a model for gas accretion from <cit.>, to calculate the rate of accretion. Our main goal is to use a more rigorous and robust gas accretion scheme that is valid and consistent for a larger range of planetary masses. In this model, the size of the accretion zone is obtained taking into account the planet's Bondi radius and the Hill radius; in this way, physical characteristics of the disc, 
such as its thermal velocity, become important to determine which cells belong to the accretion region. The timescales involved in accreting gas onto a planet are also calculated cell-by-cell using the orbital period about the planet instead of the predefined orbital period in the predecessor code. This modification allows us to obtain an accretion rate in the runaway regime that is closer to the Bondi accretion limit. 

Also, we take into account the Kelvin-Helmholtz collapse time-scale in approximation of the planet's gas accretion rate for small masses. This modification allows us to study the gas accretion process for planets which have not yet reached the runaway regime.

As expected, the reduced mass in the vicinity of the planet reduces migration rate for accreting planets. The difference in the values of the migration rates becomes less significant with the increase of the planetary mass (Fig. <ref>, right) and with the decrease of disc accretion rate (Fig. <ref>, center).
Also, the effect of the gas accretion is more significant in higher viscosity regimes (Fig. <ref>, left). In fact, in low viscous discs, the mass tends to accumulate faster in the regions around the planet and accretion can not account for cleaning the region in order to  change significantly the differential torque. Minor planetary masses are more affected from the accretion process, because the smaller the mass of the planet, the smaller the gravitational effect that prevents the mass from entering the gap, therefore, the process of accretion tends to clean the region of the planet favoring the appearance of the planetary gap.





In this work, we study only 2D discs. For massive planets, the averaged gap
profile is identical for 2D and 3D discs and we may expect very
similar effects for the planetary migration. Nevertheless, it may be interesting to perform
3D simulations and analyze the migration properties because little is known yet 
about the mechanism of gas accretion in the 3D case.

The size of the accretion zone depends on the thermal velocity of the gas. 
Recent works show the importance of considering the thermodynamical effects involved in the gas accretion process. Changes in entropy can change the luminosity and the geometry of the accretion zone, thus modifying the gas accretion rate of the planet. Recent 3D hydrodynamic simulations show that supersonic shock fronts are important for analyzing the thermodynamics of the accretion process <cit.>. 
Thus, in future works, it would be interesting to explore the changes
that radiative discs could provoke in the accretion zone. Models with larger numerical
resolutions would be needed in this case.

We test the model of gas accretion for the runaway regime by considering a 1M_J mass planet migrating
towards the star and accreting mass during its motion. We note that the runaway regime
corresponds to that used by <cit.> for parameter k_a = 2.0. This agrees
with the expected result given that the accretion of gas is close to the saturated regime for this parameter.

Also, we simulated the growth of a low mass planet fixed at position r = 1 code units and we evaluated
the time to enter the runaway regime and open the gap needed for type II migration. We found a time of
approximately 1800 orbits, which is very high considering the rapid migration of type I.
Thus, the planet would fall in the star before being able to open the gap. However, this has to be
yet further investigated. In fact, new parameters for the gas disc may reveal a faster gap
opening process, specially on heat transfer discs. In addition, a complete model that includes the
accretion of planetesimals could give us an indication of when the gas accretion process begins.
These factors will be explored in future work.



§ ACKNOWLEDGEMENTS


The authors are grateful to P. Benítez-Llambay for meaningful discussions and suggestions.
This work was supported by FAPESP (Brazil) through the grants 2014/00492-3.
This work has used the computing facilities of the Laboratory of Astroinformatics
(IAG/USP, NAT/Unicsul), that were purchased thanks to the Brazilian
agency FAPESP (grant 2009/54006-4) and the INCT-A.






















99
@urlcharsothermakeother $&#_%
@doi@urlcharsother ifnextchar [ @doi@
  @doi@[]
@doi@[#1]#2tempa#1tempaempty http://dx.doi.org/#2 doi:#2http://dx.doi.org/#2 #1
@eprint#1#2@eprint@#1:#2::nil
@eprint@arXiv#1http://arxiv.org/abs/#1 arXiv:#1
@eprint@dblp#1http://dblp.uni-trier.de/rec/bibtex/#1.xml
  dblp:#1
@eprint@#1:#2:#3:#4niltempa #1tempb #2tempc
  #3tempc empty tempc tempb tempb tempa tempb empty tempb arXivifundefined
  mn@eprint@tempbtempb:tempcmn@eprint@tempbtempc

[Alibert, Mordasini, Benz  &
  WinisdoerfferAlibert et al.2005]Alibert2005
Alibert Y.,  Mordasini C.,  Benz W.,   Winisdoerffer C.,  2005, @doi
  [] 10.1051/0004-6361:20042032, http://adsabs.harvard.edu/abs/2005A

[ArmitageArmitage2010]Armitage2010
Armitage P. J.,  2010, Astrophysics of Planet Formation

[Baruteau & MassetBaruteau &
  Masset2013]Baruteau2013a
Baruteau C.,  Masset F.,  2013, in Souchay J.,  Mathis S.,   Tokieda
  T.,  eds,  Lecture Notes in Physics, Berlin Springer Verlag Vol. 861, Lecture
  Notes in Physics, Berlin Springer Verlag. p. 201 (@eprint arXiv
  1203.3294), @doi10.1007/978-3-642-32961-6_6

[Baruteau & PapaloizouBaruteau &
  Papaloizou2013]Baruteau2013b
Baruteau C.,  Papaloizou J. C. B.,  2013, @doi []
  10.1088/0004-637X/778/1/7, http://adsabs.harvard.edu/abs/2013ApJ...778....7B 778, 7

[Benítez-Llambay, Masset  &
  BeaugéBenítez-Llambay et al.2011]Benitez2011
Benítez-Llambay P.,  Masset F.,   Beaugé C.,  2011, @doi
  [] 10.1051/0004-6361/201015774, http://adsabs.harvard.edu/abs/2011A

[Benítez-Llambay, Masset,
  Koenigsberger  & SzulágyiBenítez-Llambay
  et al.2015]Benitez2015
Benítez-Llambay P.,  Masset F.,  Koenigsberger G.,
  Szulágyi J.,  2015, @doi [] 10.1038/nature14277, http://adsabs.harvard.edu/abs/2015Natur.520...63B 520, 63

[BondiBondi1952]Bondi1952
Bondi H.,  1952, @doi [] 10.1093/mnras/112.2.195, http://adsabs.harvard.edu/abs/1952MNRAS.112..195B 112, 195

[Bryden & LinBryden &
  Lin1999]Bryden1999
Bryden G.,  Lin D. N. C.,  1999, in Bulletin of the American Astronomical
  Society. p. 1130

[Dürmann & KleyDürmann &
  Kley2015]Durman2015
Dürmann C.,  Kley W.,  2015, @doi []
  10.1051/0004-6361/201424837, http://adsabs.harvard.edu/abs/2015A

[Dürmann & KleyDürmann &
  Kley2017]Durmann2017
Dürmann C.,  Kley W.,  2017, @doi []
  10.1051/0004-6361/201629074, http://adsabs.harvard.edu/abs/2017A

[EdgarEdgar2004]Edgar2004
Edgar R.,  2004, @doi [] 10.1016/j.newar.2004.06.001, http://adsabs.harvard.edu/abs/2004NewAR..48..843E 48, 843

[Fortier, Alibert, Carron, Benz  &
  DittkristFortier et al.2013]Fortier2013
Fortier A.,  Alibert Y.,  Carron F.,  Benz W.,   Dittkrist K.-M.,
  2013, @doi [] 10.1051/0004-6361/201220241, http://adsabs.harvard.edu/abs/2013A

[Goldreich & TremaineGoldreich &
  Tremaine1980]Goldreich1980
Goldreich P.,  Tremaine S.,  1980, @doi [] 10.1086/158356, http://adsabs.harvard.edu/abs/1980ApJ...241..425G 241, 425

[Ida & LinIda &
  Lin2004]Ida2004a
Ida S.,  Lin D. N. C.,  2004, @doi [] 10.1086/381724, http://adsabs.harvard.edu/abs/2004ApJ...604..388I 604, 388

[Ikoma, Nakazawa  & EmoriIkoma
  et al.2000]Ikoma2000
Ikoma M.,  Nakazawa K.,   Emori H.,  2000, @doi []
  10.1086/309050, http://adsabs.harvard.edu/abs/2000ApJ...537.1013I
  537, 1013

[KleyKley1999]Kley1999
Kley W.,  1999, @doi [] 10.1046/j.1365-8711.1999.02198.x, http://adsabs.harvard.edu/abs/1999MNRAS.303..696K 303, 696

[Kley & NelsonKley &
  Nelson2012]Kley2012
Kley W.,  Nelson R. P.,  2012, @doi []
  10.1146/annurev-astro-081811-125523, http://adsabs.harvard.edu/abs/2012ARA

[Lin & PapaloizouLin &
  Papaloizou1986]Lin1986
Lin D. N. C.,  Papaloizou J.,  1986, @doi [] 10.1086/164653, http://adsabs.harvard.edu/abs/1986ApJ...309..846L 309, 846

[Lin & PapaloizouLin &
  Papaloizou1993]Lin1993
Lin D. N. C.,  Papaloizou J. C. B.,  1993, in Levy E. H.,  Lunine
  J. I.,  eds, Protostars and Planets III. pp 749–835

[LissauerLissauer1993]Lissauer1993
Lissauer J. J.,  1993, @doi [] 10.1146/annurev.aa.31.090193.001021,
  http://adsabs.harvard.edu/abs/1993ARA

  [Marley, Fortney, Hubickyj, Bodenheimer
   & LissauerMarley et al.2007]Marley2007
Marley M. S.,  Fortney J. J.,  Hubickyj O.,  Bodenheimer P.,
  Lissauer J. J.,  2007, @doi [] 10.1086/509759, http://adsabs.harvard.edu/abs/2007ApJ...655..541M 655, 541
  
[MassetMasset2000]Masset2000A
Masset F.,  2000, @doi [] 10.1051/aas:2000116, http://adsabs.harvard.edu/abs/2000A

[Meyer-Vernet & SicardyMeyer-Vernet &
  Sicardy1987]Meyer1987
Meyer-Vernet N.,  Sicardy B.,  1987, @doi []
  10.1016/0019-1035(87)90011-X, http://adsabs.harvard.edu/abs/1987Icar...69..157M 69, 157
  
[Mordasini, Alibert, Georgy, Dittkrist,
  Klahr  & HenningMordasini et al.2012]Mordasini2012
Mordasini C.,  Alibert Y.,  Georgy C.,  Dittkrist K.-M.,  Klahr H.,
  Henning T.,  2012, @doi [] 10.1051/0004-6361/201118464, http://adsabs.harvard.edu/abs/2012A

[Nelson, Papaloizou, Masset  &
  KleyNelson et al.2000]Nelson2000
Nelson R. P.,  Papaloizou J. C. B.,  Masset F.,   Kley W.,  2000,
  @doi [] 10.1046/j.1365-8711.2000.03605.x, http://adsabs.harvard.edu/abs/2000MNRAS.318...18N 318, 18

[Ogilvie & LubowOgilvie &
  Lubow2006]Ogilvie2006
Ogilvie G. I.,  Lubow S. H.,  2006, @doi []
  10.1111/j.1365-2966.2006.10506.x, http://adsabs.harvard.edu/abs/2006MNRAS.370..784O 370, 784

[Pollack, Hubickyj, Bodenheimer,
  Lissauer, Podolak  & GreenzweigPollack et al.1996]Pollack1996
Pollack J. B.,  Hubickyj O.,  Bodenheimer P.,  Lissauer J. J.,
  Podolak M.,   Greenzweig Y.,  1996, @doi []
  10.1006/icar.1996.0190, http://adsabs.harvard.edu/abs/1996Icar..124...62P 124, 62

[RussellRussell2011]Russell2011
Russell J. L.,  2011, PhD thesis, Canada: Univerty of Guelph, 
 <https://atrium.lib.uoguelph.ca/xmlui/bitstream/handle/10214/3212/thesis.pdf?sequence=11>

[Shakura & SunyaevShakura &
  Sunyaev1973]Shakura1973
Shakura N. I.,  Sunyaev R. A.,  1973, , http://adsabs.harvard.edu/abs/1973A
  
[Szulágyi & MordasiniSzulágyi &
  Mordasini2017]Szulagyi2017
Szulágyi J.,  Mordasini C.,  2017, @doi [Monthly Notices of the Royal
  Astronomical Society: Letters] 10.1093/mnrasl/slw212, 465, L64

[Tanaka, Takeuchi  & WardTanaka
  et al.2002]Tanaka2002
Tanaka H.,  Takeuchi T.,   Ward W. R.,  2002, @doi []
  10.1086/324713, http://adsabs.harvard.edu/abs/2002ApJ...565.1257T
  565, 1257

[Udry & SantosUdry &
  Santos2007]Udry2007
Udry S.,  Santos N. C.,  2007, @doi []
  10.1146/annurev.astro.45.051806.110529, http://adsabs.harvard.edu/abs/2007ARA

[WardWard1986]Ward1986
Ward W. R.,  1986, @doi [] 10.1016/0019-1035(86)90182-X, http://adsabs.harvard.edu/abs/1986Icar...67..164W 67, 164

[WardWard1997]Ward1997
Ward W. R.,  1997, @doi [] 10.1006/icar.1996.5647, http://adsabs.harvard.edu/abs/1997Icar..126..261W 126, 261

[ZieglerZiegler1998]Ziegler1998
Ziegler U.,  1998, @doi [Computer Physics Communications]
  10.1016/S0010-4655(98)00022-8, http://adsabs.harvard.edu/abs/1998CoPhC.109..111Z 109, 111










