© O. Schib et al. 2022

1 Introduction

Many exoplanets are observed in locations where they likely did not form (Kley & Nelson 2012). The dynamical evolution of planets embedded in protoplanetary discs is a key process in understanding planet formation (Ward 1997). For example, the existence of giant planets in mean-motion resonance implies that migration takes place (Baruteau & Masset 2013). In addition, the detection of hot Jupiters like 51 Pegasus b (Mayor & Queloz 1995) has raised the question of how they got to their present day locations as in situ formation is very challenging. However the formation of hot Jupiters is debated. Gas disc migration may explain the formation of parts of the observed population. Orbital migration is also expected to be important for the formation of warm Jupiters, giant planets close to their host star, yet further away than hot Jupiters (e.g. Anderson et al. 2018; Dawson & Johnson 2018). Fast migration, however, is not sufficient to explain the observed exoplanet locations. Processes that can slow down or halt migration are equally important, otherwise planets would tend to migrate all the way inwards and merge with their host star. One important process that can slow down orbital migration is gap formation. By interacting tidally with the disc a planet can reduce the surface density in its orbit substantially, which in turn decreases its migration rate. Broadly speaking, migration can be divided into two major regimes. As long as the planet has a mass lower than ≈10 M_⊕−100 M_⊕ (the exact value depends on the disc’s properties and the stellar mass), it does not perturb the disc significantly. These low masses correspond to type I migration, which can be very fast (timescales of 10³−10⁴ yr). More massive planets perturb the disc, pushing disc material away from their orbits, and therefore open a gap in the disc. Gap opening results in slower migration, a regime known as type II migration. In this regime, migration proceeds more slowly, typically on a 10⁵ yr timescale.

This is, however, a strongly simplified picture. A deeper understanding of planetary migration can only be gained by studying further processes. In the type I regime the total torque can be written as a sum of Lindblad and corotation torques. The Lindblad torque results from the exchange of angular momentum at Lindblad resonances in the disc and usually leads to inward migration (Goldreich & Ward 1973; Tanaka et al. 2002).

Corotation torques are a result of the gravitational interaction of the planet with disc material in the corotation region. They can be directed both inwards and outwards, and can even be strong enough to cause a net outward migration of the planet. Here the thermodynamics of the disc are important. If the disc does not cool efficiently, the net torque can be outward (Paardekooper & Mellema 2006). This is relevant only for lower mass planets, however, since the corotation torque saturates at sufficiently high masses and the outward migration cannot be sustained (Paardekooper et al. 2011). These effects can be studied by calculating torques in isothermal (Tanaka et al. 2002) or non-isothermal (Paardekooper et al. 2010) discs, and the resulting torques can also be applied in 1D models. It is necessary, however, to confirm these results in 2D and 3D hydrodynamic simulations. In the following we give some examples. Masset & Casoli (2010) calculate formulae of the saturated type I torque and perform 2D simulations to verify their results. Paardekooper et al. (2011) provide a formula for the non-isothermal type I torque that includes the effects of diffusion and show that the torque agrees well with their 2D hydrodynamic simulations. Kley & Crida (2008) perform 2D simulations of embedded planets and include heating, cooling, and radiative diffusion to study the magnitude and direction of migration. Baruteau et al. (2011) perform 2D hydrodynamic simulations of migrating planets in self-gravitating discs. They show that planets formed by fragmentation in the outer disc are likely to migrate inwards very rapidly in type I. Paardekooper & Mellema (2006) perform 3D hydrodynamic simulations to study the effect of a proper energy balance on the interaction of a low-mass planet with its disc. They demonstrate that migration can be directed outwards if the disc’s opacity is high enough. Another detailed 3D numerical study of the effect of outward migration can be found in Kley et al. (2009). The authors demonstrate that this effect can be even stronger than in 2D when the same opacity law is used.

There is a wealth of literature on these topics. Reviews can be found in Lubow & Ida (2010), Kley & Nelson (2012), Baruteau & Masset (2013), Baruteau et al. (2014, 2016).

Migration is an inherently three-dimensional process and it has been investigated extensively including 2D and 3D hydrodynamic simulations, as described. The various analytic formulae for the migration rate that have resulted from these studies can be used in 1D models. In the type I regime, these typically depend on the slopes of surface density and temperature in the disc. In a disc with multiple planets carving deep or partial gaps, these quantities are very hard to determine. Furthermore, the torque exerted on the disc material by the planet is not included in this picture. While this is a good approximation for low-mass planets, it does not allow for a gap to be formed. In our work we model the gap formation explicitly by modelling the evolution of the disc’s surface density in reaction to the presence of a growing massive planet.

In the classical picture of type II migration, planets migrate on a viscous timescale. This is not always the case, however. Mass can also cross the gap and migration may be faster or slower than the viscous inflow velocity (Dürmann & Kley 2015). It is still very challenging to model the transition between type I and type II migration. While conditions for gap opening have been derived (e.g. Crida et al. 2006), if the planet’s migration timescale is shorter than the gap-opening timescale, a gap may never open. In massive discs, an additional very rapid runaway type of migration may arise: type III migration (Masset & Papaloizou 2003, see also Malik et al. 2015).

A key feature of the observed exoplanet population is its diversity. Any successful planet formation theory must explain this diversity, which is only possible statistically. For computational reasons, the necessary parameter studies and population synthesis calculations can only be performed using low-dimensional models. One class of such models that is often applied are the 1D vertically integrated models (e.g. Nakamoto & Nakagawa 1994; Hueso & Guillot 2005). The advantage of these models is their low computational cost; the evolution of protoplanetary discs can be studied from formation to dispersal in a large parameter space. Modelling migration in 1D has its limitations. Dynamical processes, like the co-orbital dynamics, cannot be modelled directly in 1D since they are inherently three-dimensional. Therefore, it may be necessary to introduce additional approximations or parametrisations that are valid only in some regimes. At the same time, the location in the disc where gaps form is paramount for the predicted population of planets. A few simple prescriptions, for example for the depth and shape of the gap (Kanagawa et al. 2017; Fung et al. 2014) and the migration process (Ida et al. 2020) have been suggested. The difficulty with these approaches is that they often rely on the knowledge of the unperturbed surface density near the planet. In practice, this property is only known when the disc evolution starts from well-defined initial conditions. After a period of evolution, or when the disc formation is included through infall from the molecular cloud core, or when several gap-opening planets perturb the disc, the unperturbed surface density is meaningless or difficult to determine.

The goal of this paper is to derive a prescription for migration and accretion in 1D disc models, calibrated by the results of 3D hydrodynamic simulations. In order to overcome the difficulties described above, this prescription uses torque densities and relies only on directly accessible disc parameters, and self-consistently models the exchange of mass and angular momentum. We do not attempt to reproduce the exact gap shapes found in hydro-dynamic simulations. To what degree this is possible when considering torque densities in the context of accreting migrating planets is currently unclear. We hope to address this important topic in the future.

The paper is organised as follows. In Sect. 2 we describe the numerical set-up. Section 3 discusses the investigated parameter space and our initial conditions. In Sect. 4 we present our results in comparison with those from hydrodynamic simulations. In Sect. 5 we perform a comparison with a different prescription found in the literature. Section 6 contains a discussion of the results and the limitations of the model. In Sect. 7 we summarise and conclude.

2 Model description

We apply the disc formation and evolution model from Schib et al. (2021). It describes the temporal evolution of a rotationally symmetric 1D disc of gas, described by the vertically integrated surface density Σ ≡ Σ(r, t). Here we introduce a planet on a circular orbit. We use cylindrical polar coordinates, with r denoting the radial direction. The disc’s mid-plane is located at z = 0 au. In this section we briefly review the model in the form in which it is applied in this work. We note that we use the term ‘planet’ throughout this work for simplicity, even though ‘protoplanet’ or ‘proto-brown dwarf’ may be more appropriate in some cases.

2.1 Disc evolution

The evolution of a protoplanetary disc is described by the viscous evolution equation (Lüst 1952; Lynden-Bell & Pringle 1974), adding the effects of angular momentum injection by a planet and mass sinks $$\frac{\partial \Sigma }{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\left(v\text{Σ}{r}^{1/2}\right)-\frac{2\text{ΛΣ}}{\text{Ω}}\right]-S_{\text{acc}},$$(1)

where S_acc is a sink term for the mass accreted by planets (other sink terms such as photoevaporation are not considered here) and Ω denotes the Keplerian angular frequency of the disc. The expression for S_acc is given in Sect. 2.2. The term 2ΛΣ/Ω describes the gravitational interaction between the planet and the disc, leading to angular momentum exchange. The torque density distribution Λ ≡ −dT/dm used in this work is discussed in Sect. 2.3. In our convention, T = 2π ${\displaystyle {\int }_0^{\infty }\frac{d\mathcal{T}}{dm}\left(r\right)\sum \left(r\right)r}$ dr is the total torque on the planet. Equation (1) is solved on a grid of 2800 logarithmically spaced cells extending from 0.03 to 30 000 au. We performed a resolution test where we increase or decrease the resolution by a factor of two and found that it has a negligible effect on our results. We use the solver from Birnstiel et al. (2010).

We assume hydrostatic equilibrium in the vertical direction and define the pressure scale height H through $$\rho \left(r,z\right)={\rho }_0\left(r\right)\exp \left(-\frac{z^2}{2H{\left(r\right)}^2}\right).$$(2)

Here, ρ is the density in the disc, and its midplane value ρ₀(r) is related to Σ(r) through $$\Sigma \left(r\right)={\rho }_0\left(r\right)H\left(r\right)\sqrt{2\pi }.$$(3)

In this paper we assume a constant aspect ratio H/r. The disc’s scale height is related to the isothermal sound speed $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu \text{u}}}$$(4)

through $$H=\frac{c_{\text{s}}}{\text{Ω}}.$$(5)

In Eq. (4) k_B is the Boltzmann constant, µ ≡ 2.3 the mean molecular weight, and u the atomic mass unit.

2.2 Gas accretion

We study gas accretion in the disc-limited regime (accretion is limited by the disc’s mass reservoir, but not the planet’s ability to accrete). For low-mass planets embedded in a disc, their Kelvin-Helmholtz contraction (Ikoma et al. 2001) limits gas accretion (for example Sect. 3.8.1 in Mordasini et al. 2012). In this phase our model represents an upper limit.

Our accretion model is based on the Bondi and Hill accretion described in D'Angelo & Lubow (2008, 2008). The authors study migrating and accreting planets using 3D nested-grid hydrodynamical simulations of locally isothermal discs. We note that the locally isothermal assumption is likely to influence the gas accretion. However, it is still unclear whether a departure from isothermality increases or decreases the accretion rate. For example, Ayliffe & Bate (2009) find that accretion rates are highest in a locally isothermal disc, while Machida et al. (2010) find the opposite, although both of these studies perform radiation hydrodynamics simulations of accretion. Based on their numerical simulations, D'Angelo & Lubow (2008) parametrise the accretion rate onto a planet with mass M_p in their Eq. (9): $${\dot{M}}_{\text{p}}\sim \frac{\text{Σ}}{H}\text{Ω}{R}_{\text{f}}^3.$$(6)

Here R_f is the feeding zone radius, taken to be either the Bondi radius $R_{\text{B}}=G{M}_{\text{p}}/c_{\text{s}}^2$ or the Hill radius R_H = r_p [M_p / (3M_*)^1/3, where the sound speed is evaluated at the planet’s semi-major axis r_p. We assume circular orbits. The inverse growth timescales M_p/M_p (equal to ${\tau }_{\text{B}}^{-1}$ in the Bondi regime, and ${\tau }_{\text{H}}^{-1}$ in the Hill regime) then become $$\frac{1}{{\tau }_{\text{B}}}=C_{\text{B}}\text{Ω}\frac{\text{Σ}{r}^2}{M_{*}}{\left(\frac{r_{\text{p}}}{H}\right)}^7{\left(\frac{M_{\text{p}}}{M_{*}}\right)}^2,$$(7) $$\frac{1}{{\tau }_{\text{H}}}=\frac{1}{3}{C}_{\text{H}}\text{Ω}\frac{\text{Σ}{r}^2}{M_{*}}\left(\frac{r_{\text{p}}}{H}\right),$$(8)

where C_B and C_H are dimensionless coefficients of order unity. DL08 found C_B = 2.6 and C_H = 0.89 to agree best with their hydrodynamic simulations. Since our prescription of accretion differs from theirs, we have to use different values for C_B and C_H, as discussed later. The overall inverse growth timescale is defined as $$1/{\tau }_{\text{G}}=\{\begin{array}{cc}1/{\tau }_{\text{B}}& M_{\text{p}}<M_{\text{t}}\\ 1/{\tau }_{\text{H}}& M_{\text{p}}\ge M_{\text{t}},\end{array}$$(9)

where $$M_{\text{t}}=\frac{M_{*}}{\sqrt{3}}\sqrt{\frac{C_{\text{H}}}{C_{\text{B}}}}{\left(\frac{H}{r_{\text{p}}}\right)}^3$$(10)

is the transition mass when τ_H = τ_B.

Zhu et al. (2012) perform 2D hydrodynamic simulations of self-gravitating protostellar discs with infall. They find that the accretion rate M_C onto their clumps agrees reasonably well with their Eq. (14): $${\dot{M}}_{\text{c}}=4\text{ΣΩ}{R}_{\text{H}}^2.$$(11)

This expression agrees with Eq. (8) up to a factor 4/C_H when R_H = H. We can therefore use this accretion scheme for clumps and for compact planets, and irrespective of whether the disc is self-gravitating and/or subject to infall. We note, however, that the size of the feeding zone is likely to be smaller in the self-gravitating regime (Shabram & Boley 2013).

We next apply Eq. (6) to our 1D model, refining the approach. Instead of using global values of Σ and Ω, we calculate the contributions from each grid cell inside the feeding zone separately. The accreted mass is then removed self-consistently from the disc at the correct location. We obtain the following accretion rate: $$\begin{array}{ll}{\dot{M}}_{\text{gas}}\hfill & =2C_{\text{B,H}}{\displaystyle {\int }_{r_{\text{p}}-R_{\text{f}}}^{r_{\text{p}}+R_{\text{f}}}{\displaystyle {\int }_0^{\sqrt{R_{\text{f}}^2-{\left(r-r_{\text{p}}\right)}^2}}\rho \left(r,z\right)\text{d}z{\upsilon }_{\text{rel}}\text{\hspace{0.17em}d}r}}\hfill \\ \hfill & =\sqrt{2\pi }{C}_{\text{B,H}}{\displaystyle {\int }_{r_{\text{p}}-R_{\text{f}}}^{r_{\text{p}}+R_{\text{f}}}{\rho }_0\left(r\right)H\left(r\right)\text{erf}\left(\frac{\sqrt{R_{\text{f}}^2-{\left(r-r_{\text{p}}\right)}^2}}{\sqrt{2}H\left(r\right)}\right){\upsilon }_{\text{rel}}}\text{\hspace{0.17em}d}r.\hfill \end{array}$$(12)

In the last step we substituted ρ(r, z) from Eq. (2). The factor 2 at the beginning takes into account that we only integrate over half of the circular feeding area. The coefficient C_B,_H is equal to C_B in the Bondi regime and C_H in the Hill regime. The expression $\sqrt{R_{\text{f}}^2-{\left(r-r_{\text{p}}\right)}^2}$ denotes the distance from the mid-plane to the top or bottom of the feeding zone in z-direction.

A schematic of the geometry is given in Fig. 1. The quantity $v_{\text{rel}}=\left|r\Omega -r_{\text{p}}{\Omega }_{\text{p}}\right|$ is the gas’ velocity relative to the planet. The presence of the planet introduces a perturbation to the Keplerian flow of the gas. Therefore, using Ω in the expression for ν_rel is potentially problematic near the planet (e.g. Lubow et al. 1999). However, the form of Eqs. (6) and (11), and the agreement of these expressions with hydrodynamic simulations, suggests that the total accretion rate should not be very different. When the planet mass is low there is no strong perturbation of the disc. For massive planets the surface density near the planet is depleted and does not contribute significantly to the accretion. Therefore, this assumption is not expected to affect our results. We note that if we tried to reproduce gap shapes from hydrodynamic simulations accurately, these aspects would become important. The remaining integral in Eq. (12) is evaluated numerically. The corresponding mass is removed from the disc and added to the planet in every grid cell separately at each time step. Comparing Eq. (12) with the requirement ${\dot{M}}_{\text{gas}}=2\pi {\displaystyle {\int }_{r_{\text{p}}-R_{\text{f}}}^{r_{\text{p}}+R_{\text{f}}}r{S}_{\text{acc}}}$ dr gives the sink term in Eq. (1): $$S_{\text{acc}}\left(r,t\right)=\frac{1}{\sqrt{2\pi }}\frac{H}{r}{\rho }_0\left(r\right)C_{\text{B,H}}\text{erf}\left(\frac{\sqrt{R_{\text{f}}^2-{\left(r-r_{\text{p}}\right)}^2}}{\sqrt{2}H}\right){\upsilon }_{\text{rel}}.$$(13)

Equations (6) (DL08) and (11) (Zhu et al. 2012) make the implicit assumption that the planet’s radius is much smaller or larger, respectively, than the disc’s scale height. In the first case the cross section is $\pi R_{\text{f}}^2$; in the second it is 2R_f H. These expressions do not apply when the planet’s size is comparable to the scale height. We solve this problem by using the analytic expression for the vertical structure (Eq. (2)) in the calculation of ${\dot{M}}_{\text{gas}}$. Equation (12) is thus valid when the planet is smaller, larger, or comparable to the disc scale height.

In the Hill regime we studied two different configurations for the feeding zone radius: R_f = R_H and R_f = 3R_H. We found that R_f = R_H depletes the feeding zone too quickly, and that R_f = 3R_H reduces this effect (see Sect. 4.1). Using a different feeding zone radius means the C_H also needs to be different in order to obtain the same inverse growth timescale. The values we adopted for the different configurations studied are given in Table 1. In order to prevent unphysical effects due to a discontinuous jump in R_f, we apply a smooth transition between the Bondi and Hill regimes.

Fig. 1 Schematic of the accretion geometry. The feeding zone is centred at the location of the planet rp, at z = 0. The direction tangential to the planet’s motion is perpendicular to the page.

2.3 Migration

Our description of migration is motivated by the impulse approximation (Lin & Papaloizou 1979a,b, 1986). This ensures a two-way angular momentum exchange between disc and planet. Instead of using the classical (analytical) impulse approximation, we apply a more modern formalism. We investigated how the classical impulse approximation, as well as a modified version, behave in Appendix B. D'Angelo & Lubow (2010; 2010) perform 3D hydrodynamic simulations of locally isothermal discs to study the type I migration torque. We follow them in assuming that the torque distribution per unit disc mass has the form $$\frac{\text{d}\mathcal{T}}{\text{d}m}=ℱ\left(x,\beta ,\zeta \right)\Omega {\left(r_{\text{p}}\right)}^2{r}_{\text{p}}^2{\left(\frac{M_{\text{p}}}{M_{*}}\right)}^2{\left(\frac{r_{\text{p}}}{H\left(r_{\text{p}}\right)}\right)}^4,$$(14)

where ℱ is a dimensionless function describing the torque’s shape, and the radial scaling factor is $$x=\frac{r-r_{\text{p}}}{\max \left(H,R_{\text{H}}\right)}.$$(15)

The parameters β and ζ are the gradients of the surface density and temperature, respectively, such that Σ ∝ r^−β and T ∝ r−ζ. DL10 find an analytical fit for ℱ : $$ℱ\left(x,\beta ,\zeta \right)=\left[{\rho }_1{e}^{-{\left(x+p_2\right)}^2/p_3^2}+p_4{e}^{-{\left(x-p_5\right)}^2/p_6^2}\right]\tan h\left(p_7-p_{8}x\right).$$(16)

In Eq. (16) p_i ≡ p_i(ß, ζ), i ∈ {1…8} are parameters that DL10 find by fitting to their hydrodynamic simulations for a large grid of disc parameters. The values we use in this study are listed in Table 2. The first column for (β, ζ) = (0.5, 1) contains the values appropriate for our comparison with DL08 (see Sect. 3). The second column for (β,ζ) = (1, 0.5) is obtained using a linear interpolation from the values for a combination of (β, ζ) of (1, 0), (0.5, 1), and (1.5, 1) since no parameters are available for (1, 0.5) in DL10. We used this second set of parameters for our comparison with Fletcher et al. (2019) (FNS19) in Sect. 4.8.

Using these fitting parameters seems to defeat our goal to find a prescription of migration independent of the slopes in Σ and T. However, in practice it is sufficient to choose the fitting parameters approximately adequate to the simulation performed since the torque density distributions do not depend very sensitively on the precise values of the slopes (see Sect. 6 in DL10). We also tested whether this is the case by varying either β or ζ by ±50 or 25%, respectively and found that the difference in semi-major axis migrated is at most 12%. This sensitivity calculation can be found in Appendix C. Therefore, our prescription can also be used for discs with with (β, ζ) ≠ (1, 0), (0.5, 1), as long as the difference does not exceed a few tens of percent. This also depends on the required precision. As an example, Appendix C.1 discusses the impact of changing β or ζ by 0.25 in either direction. We note that in non-isothermal discs the dependency of migration on the temperature slope, for example, is no longer small. Unfortunately, grids of torque densities like those given in DL10 for locally isothermal discs are not available for non-isothermal discs. We plan to improve our prescription if such torque densities become available in the future.

Using torque densities in 1D models may result in premature gap-opening, as described in Hallam & Paardekooper (2017). The authors study gap formation in 1D and 2D models. They propose modifying the torque density, setting it to zero for |r − r₀| < 0.95 R_H, with a sharp transition for 0.95 < |r − r₀| < 1.05 R_H in order to prevent the formation of gaps that are too deep. We followed this approach, but made the region where the torque density is truncated larger, choosing 1.8 R_H instead of R_H. This is necessary because leaving it at 1 R_H still caused a gap to open too soon compared to the hydrodynamic simulations (Sect. 4.1). We call this modified torque density ‘torque mod’. The key here is that we only used the modified torque density in the torque acting on the disc. The migration rate is still calculated based on the unmodified torque. This is important since using the modified torque also for the calculation of the migration rate would again result in too slow migration¹. The advantage of this approach is that it can be applied to the entire range of torque densities given in DL10. There is a drawback: using different planet-disc and disc-planet torques violates the conservation of angular momentum. For this reason, we also investigated a different approach for the torque density. This torque model, which we call ‘high mass torque’ is based on Fig. 15 in DL10. In this figure, the authors show torque density distributions for a planet of increasing mass from their 3D hydrodynamic simulations. The distributions are scaled by ${\Omega }^2{r}_{\text{p}}^2$(M_p/M_*)²(r_p/H)⁴ and demonstrate the influence of the planet’s mass. The figure reveals that the scaled torque density changes and decreases in magnitude with increasing planet mass, in particular close to the planet. Starting at approximately 1 M (1 Jupiter mass), an inversion appears: the torque density takes on negative values inside the planet location and then reaches positive values outside the planet location. These features are also seen in Fig. F.1, where we plot the scaled high mass torque for various masses. As pointed out by DL10, these subtleties are not important once a deep gap has been opened in the disc. However, the inversion is key during the gap formation process, as we discuss in Sect. 4.1. The torque densities given in Fig. 15 of DL10 are for (β, ζ) = (0.5, 1) We use ℱ(x, 0.5, 1) for a planet of zero mass and interpolate linearly in mass using the (digitised) torque densities from the figure. The numerical data we used for the interpolation is available in tabulated form at the CDS. Figure 2 shows the scaled torque densities for the three cases we discussed. The most massive planet for which a torque density is given in DL10‘s Fig. 15 has a mass of 2$M_{0_{+}}$. We use this torque density also for more massive planets (no extrapolation). The inversion, as it develops with increasing planet mass, is additionally shown in Fig. F.1. Figure F.1 also shows the LP86 formula and the A02 formula (scaled for visibility) that are discussed in Appendix B.

As an example, we show the time evolution of the surface density of an accreting, migrating planet in Fig. 3. It illustrates, how the planet interacts with the disc over time. The “High mass torque” is used in this example, and the initial conditions are the ones from our baseline case that we will discuss in Sect. 4.1. The evolution of the surface density in the other cases we studied is presented in Appendix A. In this study we do not discuss the precise shape of the gap since this data is not available from DL08.

Fig. 2 Torque densities used in this work: ℱ(0.5, 1) (DL10), a modified version (torque mod) where ℱ is truncated near the planet when applied to the disc, and a variant model (high mass torque) based on an interpolation of torque densities obtained from a hydrodynamical simulation. The last case is for a 1.5 M planet. The torque densities are scaled by ${\Omega }^2{r}_{\text{p}}^2$(Mp/M*)2(rp/H)4. Additionally, the LP86 formula (thin black solid line) and the A02 formula (thin black dashed line) are shown (see Appendix B). These are divided by 3 and 15, respectively for better visibility.

Fig. 3 Surface density around an accreting, migrating planet at different times. At the beginning, the planet has a low mass (5 M⊕) and does not perturb the surface density appreciably. As the planet grows, it changes the surface density in its vicinity both by accreting gas and by pushing disc material away through tidal interaction. After 1000 orbits, the planet has migrated to ≈4.2 au, reached a mass of 1 $M_{0_{+}}$ and opened a gap in the disc. The planet’s initial location (5.2 au) is indicated with a thin black vertical line. This simulation is described in more detail in Sect. 4.1.

2.4 Autogravitation

In the present study the disc mass is low compared to the stellar mass. In this scenario the disc’s self-gravity (autogravitation) can be neglected. If the disc is comparable in mass to its host star, the disc’s scale height and angular frequency will change (Hueso & Guillot 2005). We already studied this effect in detail in Schib et al. (2021). Our prescriptions presented here can in principle still be applied. However, some of the equations given in Sect. 2 need to be modified. In particular, H needs to be replaced by c_s/Ω in Eq. (14). The expressions for accretion also change slightly due to the different vertical structure. The details are given in Appendix D.

3 Investigated parameter space and initial conditions

In this section we describe the set-up of our simulations for the different comparisons we performed.

For the comparison with DL08 (Sects. 4.1–4.5) we chose initial conditions similar to theirs. Our simulations start with a surface density profile of the form $$\text{Σ}\left(r\right)={\text{Σ}}_0{\left(\frac{r}{5.2\text{\hspace{0.17em}au}}\right)}^{{\beta }_{\text{Σ}}}\exp \left[{\left(\frac{5.2au}{r_{\max }}\right)}^{2+{\beta }_{\text{Σ}}}\right]\exp \left[-{\left(\frac{r}{r_{\max }}\right)}^{2+{\beta }_{\text{Σ}}}\right],$$(17)

where β_Σ = −0.5 is the surface density radial slope (the same as in DL08), and r_max is the radius of the outer exponential cutoff, for which we chose 30 au. The inner disc is truncated at 0.05 au. We note that for numerical reasons, DL08 study only a part of the disc’s radial extent (2.08−13 au for the cases relevant here). Due to the different nature of their 3D simulation, their boundary conditions are also different (see their Sect. 2.1.3). Since the viscosity is low in the majority of the simulations we performed, the effect of the different boundary conditions is not important. This changes in the case with higher viscosity, where the boundary conditions play a significant role, as we discuss in Sect. 4.5. We checked that increasing or decreasing our r_max by one-third has a very small influence on our results. The initial value of the surface density at 5.2 au, Σ₀ is chosen between 100 and 500 g cm⁻² according to the different cases in DL08. Using Eq. (4) we get a temperature profile proportional to r⁻¹ : $$T\left(r\right)={\left(\frac{H}{r}\right)}^2\frac{\mu \text{u}G{M}_{*}}{k_{\text{B}}}{r}^{-1}.$$(18)

We set (β, ζ) = (0.5, 1) in Eq. (16), as suggested by the slopes of Σ and T in Eqs. (17) and (18). A constant kinematic viscosity of the same magnitude as DL08 was also used (see Sect. 4).

For our comparison with FNS19 in Sect. 4.8, we used the initial surface density given in their Fig. 1 (digitised), and the temperature profile given in their Eq. (1): $$T\left(r\right)=200K{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-1}.$$(19)

This temperature profile agrees well with the initial profile in FNS19. However, while our temperatures stay constant in time, their temperatures start to deviate inside of ~30 au due to artificial viscosity heating. We expect that this will slow down migration somewhat inside of ~40 au. Instead of using a constant viscosity as above, we use the alpha-viscosity $$v=\alpha \frac{c_{\text{s}}^2}{\text{Ω}}$$(20)

(Shakura & Sunyaev 1973) with α = 3 × 10⁻². This value approximately corresponds to the artificial viscosity of the codes used in FNS19 (their Sect. 3.3.1). Hydrodynamic codes such as the ones used in their study typically produce artificial viscosities comparable to this value.

4 Results

In this section we calibrate our model and compare our results to those inferred using hydrodynamic simulations, first to DL08 and then to FNS19. We use the baseline case (see Sect. 4.1) from DL08 to calibrate the parameters C_B, C_H, the width of the feeding zone R_f, and the truncation width for the torque mod (see Sect. 2). Once these parameters are calibrated they remain unchanged. For practical use we recommend our high mass torque model with the parameters listed in Table 1.

In the following sections we compare the results for a number of different sets of initial conditions. Table 3 gives an overview of the cases studied in the following sections.

4.1 Baseline case: Calibration of the 1D model

In this simulation a 5 M_⊕ planet is inserted at 5.2 au into a disc with an initial surface density of 100 g cm⁻² at 5.2 au with a constant kinematic viscosity of 10¹⁵ cm² s⁻¹ (corresponding to α = 0.004 at 5.2 au) and a constant aspect ratio of H/r = 0.05 (corresponding to an initial temperature of ≈118K at 5.2 au). These parameters are identical to those used by DL08.

The planet is then allowed to migrate and accrete gas. We note that the study by DL08 and our work both focus on the ‘disc-limited’ mode of gas accretion (Sect. 2.2): it is assumed that the planet can accrete all the gas available from the disc (also known as the rapid gas accretion phase). To what degree these accretion rates can be absorbed depends very much on the planetary properties. For example, in the core accretion scenario, gas only starts to be accreted rapidly once the critical core mass has been reached (Pollack et al. 1996; Emsenhuber et al. 2021). The critical core mass is of the order of 10M_⊕. The disc-limited regime is typically reached at total planet masses of 50 M_⊕ to 100 M_⊕ (Movshovitz et al. 2010; Mordasini et al. 2014).

The top left panel in Fig. 4 shows the inverse growth timescale ${\tau }_{\text{G}}^{-1}$ as a function of planet mass when applying the DL10 torque with a feeding zone radius of 1 R_H (orange solid line). The accretion agrees well with that from DL08 in the Bondi regime (we chose C_B accordingly). For the Hill regime we chose C_H such that the inverse growth timescale agrees with the analytical fit at early times (corresponding to a planet mass of ≈20M_⊕). However, the accretion rate drops too soon. This happens because a gap opens very rapidly. This effect is seen in the bottom right panel of Fig. 4. The figure shows Σ_B, the surface density near the planet, as a fraction of the initial Σ₀. The surface density is averaged over a radial band of width 2 H, centred on the planet. The value of Σ_B drops somewhat faster in our calculation (orange solid line) compared to the hydrodynamic calculation. In order to prevent the premature depletion of the feeding zone, we increased its radius by a factor of three (and reduced C_H accordingly) since this gives the best agreement with the hydrodynamical simulation. The result (green dashed line in Fig. 4) is now very similar to that from DL08. Consequently, the mass evolution is now also very similar (top right panel in the figure). The planet reaches 300 M_⊕ at the end of the hydrodynamic simulation after ≈1200 orbits.

We next investigate the behaviour of migration in our model. The bottom left panel of Fig. 4 shows the planet’s semi-major axis versus time. It reveals that the migration rate in the DL10 torque case agrees well for the first ~500 orbits, but slows down too much when approaching the type II regime (orange solid line). This slowdown is even stronger with the increased feeding zone radius (green dashed line). While the increased feeding zone gives a better agreement for the accretion, it has an opposite effect on the migration: a gap is opened even more quickly due to the faster increase in mass, and the migration is too slow. This is again seen in the bottom right panel of the figure, where Σ_B now drops even more quickly and reaches almost zero after 600 orbits. The premature gap opening in 1D models is a well-known result. It is studied in Hallam & Paardekooper (2017).

Now we consider the torque mod described in Sect. 2.3. The effect is again seen in Fig. 4 (red dotted line). The agreement with DL08 is very good both in terms of mass evolution (top right panel) and migration (bottom left panel). The decrease in surface density is slightly slower than in the hydrodynamic calculation. Finally, applying the high mass torque described in Sect. 2.3 gives a very similar result. The reduction of the torque density in this case, in combination with the inversion near the planet, produces fast migration (as in DL08) and angular momentum is conserved. The surface density near the planet is still depleted more quickly compared to DL08, but it does not inhibit mass growth or slow down migration significantly. The planet mass reached at the end of the hydrodynamic simulation is »1% higher (torque mod) and ≈10% lower (high mass torque) respectively, compared to DL08. The bottom left panel also shows predictions from saturated and unsaturated type I theory, and type II theory (thin black dashed lines). As expected, the planet initially follows the type I prediction well, but slows down as it grows in mass and starts perturbing the disc, approaching the type II expectation. The theoretical type I tracks were calculated by DL08 based on the initial value of the surface density.

We find that when we calibrate C_B, C_H, R_f and (for the torque mod prescription) the width of the truncation region, the agreement with the hydrodynamic simulation is good. In the following sections we investigate whether this agreement persists in the other cases studied by DL08.

Fig. 4 Time evolution of the baseline case with an initial surface density of 100 g cm−2 at 5.2 au. The solid black line represents the 3D result from DL08. The solid orange line corresponds to our 1D result when using the DL10 torque. The green dashed line gives again the result with the DL10 torque, but with the feeding zone radius increased threefold. The red dotted line shows a calculation where, additionally, the torque density has been modified. The purple dash-dotted line shows the variant model, where interpolated torque densities for planets of different masses from DL10 were applied. The torque densities are described in Sect. 2.3. Top left: inverse growth timescale in units of the orbital timescale at the initial location (≈12 yr). The thin black dashed line represents the fit given by Eqs. (7) and (8). Top right: mass evolution. Bottom left: orbital migration. The thin dashed black lines represent (from left to right) predictions for saturated and unsaturated type I, and type II theory. Bottom right: evolution of the mean surface density around the planet. Time is given in orbits at the planet’s initial separation (5.2 au).

Fig. 5 Same as Fig. 4, but for increased surface density (300 g cm−2 at 5.2 au).

4.2 Higher surface density

Here we consider, as in DL08, a case where the initial surface density is increased by a factor of three relative to the baseline case (300 g cm⁻² at 5.2 au), while all the other conditions were left the same. The results are shown in Fig. 5. For clarity, we omit the DL10 torque and R_f = 3 R_H cases from now on.

We find a similar behaviour to the baseline case. The increased feeding zone radius, combined with the modified torque density, results in a migration rate (bottom left panel) and an inverse growth timescale (top left panel) very similar to that found by DL08. This is true for both variant torque models we considered and is true despite the fact that, in this calculation, the planet is growing much faster, reaching 300 M_⊕ already after 220 orbits. The masses we find are ≈5 (high mass torque) and ≈10% (torque mod) higher than what DL08 found after 350 orbits.

The decrease in surface density is again somewhat slower in the torque mod case and somewhat faster in the high mass torque case, compared to the hydrodynamic simulation (bottom right panel). We conclude that the calibrations done with the baseline case lead to good agreement of 1D and 3D simulations also when at a higher initial surface density.

Fig. 6 Same as Fig. 4, but with five times higher initial surface density (500 g cm−2 at 5.2 au).

4.3 Very high surface density

In this section we present the results when the disc’s initial surface density is increased by a factor of five relative to the baseline case (500 g cm⁻² at 5.2 au). The results are presented in Fig. 6.

For this simulation DL08 do not provide the mass evolution of the planet as a function of time, so it is missing in Fig. 6. They indicate it reaches 1 M after 130 orbits, which is shown with the black cross in the top right panel of the figure. We find that the mass evolution agrees well with this value and with the analytic fitting formula (top left panel). The mass reached after 130 orbits is 2% (high mass torque) and 8% (torque mod) lower compared to DL08. In the migration track there is some disagreement. The torque mod model shows almost no departure from type I migration during the first 125 orbits. It only does so after reaching approximately 0.65 times the initial separation (not seen in the figure). The high mass torque model agrees qualitatively with the hydrodynamic simulation, though slowing down slightly faster. The behaviour of the surface density around the planet (bottom right panel) is very similar to the first two cases we discussed.

4.4 Lower temperature

In this section we look at a variant calculation with a lower disc temperature (H/r = 0.04, this corresponds to an initial temperature of ≈76 K at 5.2 au). The results are presented in Fig. 7.

The mass growth is faster here in comparison to the baseline case, due to the higher Bondi accretion rate, but it slows down quickly as the Hill regime is reached. Consequently, inward migration is also faster until a gap starts to open after ~300 orbits. The agreement of our model with the calculation from DL08 is good for both of our torque models shown. The planet migrates in marginally further, but then transitions to type II migration at a very similar rate (bottom left panel). The evolution of Σ_B is again slightly different between the torque mod and high mass torque simulations. The latter matches more closely the result from the hydrodynamic calculation (bottom right panel). The mass evolution is also similar (top panels), though the mass reached at the end of the hydrodynamical simulation after about 850 orbits is 15% (high mass torque) and 20% (torque mod) lower compared to DL08.

Fig. 7 Same as Fig. 4, but for reduced temperature (H/r = 0.04).

4.5 Higher viscosity

In the simulation presented in this section, the disc viscosity is higher by an order of magnitude (10¹⁶ cm² s⁻¹, corresponding to α = 0.04 at 5.2 au) compared to the other cases. A precise agreement with the hydrodynamic simulations from DL08 or their analytic fit is not expected here, as we discuss below. The result are shown in Fig. 8.

The top left panel of Fig. 8 shows that our Hill accretion rate does not reach the value predicted by the analytic fit. Conversely, the hydrodynamic simulation reaches a value higher by more than a factor two compared with the fit. This leads to a planetary mass that is higher by a factor of five in comparison to our model (top right panel).

The migration tracks are not very different (bottom left panel), though in DL08 the planet becomes slower than expected from type II theory after ≈500 orbits, likely due to the fast increase in mass. In our simulation the migration continues according to the type II slope.

The difference in mass accretion is also seen very clearly in the evolution of Σ_B (bottom right panel). The disc remains almost unperturbed for 400 orbits for DL08, while it starts declining steadily much earlier in our model.

These different results are caused by the global disc evolution due to the high viscosity in combination with differences in boundary conditions. The viscous timescale (r²/ν) is only ~19 kyr, or less than 1600 orbits at the planet’s initial location. The consequence is a rapid decline in the disc mass in our simulation, unlike in the other cases we consider. This has a substantial influence on the planetary evolution: mass growth is reduced due to the flow of disc material across the inner and outer boundary. This explains the difference to the result from DL08: they only model a part of the disc in radial direction. Their boundary conditions at the outer edge of the grid does not allow outflow (their Sect. 2.1.3)². In this case our results are likely to be more realistic than those presented in DL08 because we model the entire disc and our boundary conditions allow for outflow of disc material. In Sect. 4.8 we study a case of a similar viscosity and show that our model performs reasonably well. Our results for the high-viscosity case are very similar for both of the torque models discussed.

4.6 Impulse approximation

We also investigated how the classical impulse approximation (Lin & Papaloizou 1979a,b, 1986) and an improved version (Armitage et al. 2002) compare to the simulations fromDL08 in the context of our model. The impulse approximation explicitly neglects co-rotation torques from material inside the horseshoe region. Therefore, it should only be applied when a gap is already present. However, since it is widely used in the literature, it is still interesting to study how it behaves in our test cases. For conciseness, these results are discussed in Appendix B.

It is found that both prescriptions give worse agreement for accretion, migration, and gap depth than torque mod and high mass torque. The modern approach with torque densities derived from hydrodynamic simulations should thus clearly be preferred.

Fig. 8 Same as Fig. 4, but for ten times higher viscosity (10 cm2 s−1).

4.7 Type I migration

Recently, a new interpretation of type II migration has been proposed. Kanagawa et al. (2018) argue that it is actually the same as type I migration, but using the surface density in the gap. We cannot directly investigate to what degree this approach would allow type II migration in our model since Kanagawa et al. (2018) do not consider accretion. Here we studied a similar approach instead. For this test we turned off the tidal interaction between disc and planet. The migration of the planet is then calculated as a type I torque based on the surface density at the planet location. This analysis can also be found in Appendix B. In summary, we find that this approach gives good agreement with the results from DL08 early in the planet’s evolution and when gap formation starts. Migration is slowed down significantly through the reduction of Σ by mass accretion alone. However, migration typically does not slow down enough because no deep gap forms due to the lack of a type II torque acting on the disc. Therefore, this approach does not seem adequate to model the long-term evolution of planets in the framework of our model. Nevertheless our investigation shows the importance of the accretion process on the planet’s migration and gives surprisingly good agreement with DL08 for accretion and for migration. We note that this approach does not conserve angular momentum.

4.8 Large separation

So far our analysis has concentrated on initially low-mass planets that were released at 5.2 au. FNS19 perform a code comparison project where they compare the performance of a number of hydrodynamical codes. Their focus is on massive planets in the outer disc, such as would be expected from disc fragmentation (see Kratter & Lodato 2016 for a review). However, they only study discs that are not self-gravitating (Toomre-Q-parameter >2, Toomre 1964). It is found that the seven hydrodynamic codes they compared agree qualitatively on the outcome, but differ quantitatively. The authors also include a comparison with existing population synthesis studies and find that the agreement among them is worse than among the hydrodynamic codes.

In order to assess whether our method is also usable for large initial separations, we repeat one of their two test cases that include accretion. A planet of 2 M is inserted in the disc at 120 au and left to accrete and migrate. The results are shown in Fig. 9. FNS19 also perform comparison runs with three different disc instability population synthesis codes: Forgan & Rice (2013), Müller et al. (2018), and Nayakshin & Fletcher (2015). These results are also shown in the figure (digitised from their Fig. 9).

We do not expect a precise agreement with the simulations from FNS19. The goal of this study is to compare the performance of different hydrodynamic codes. To this end, FNS19 choose highly idealised initial conditions. In particular, high-mass planets are inserted into an unperturbed disc. In a more realistic case, and if the planet had formed through core accretion, the disc would be significantly perturbed, and possibly already have a gap in the presence of such a massive planet. Furthermore, accretion is modelled with a sink particle prescription.

Nevertheless, we assessed whether we can find a qualitative agreement with their work. For this comparison, we used the same initial surface density and temperature profiles as in FNS19 (digitised from their Fig. 1). For the torque mod torque density we used different p_i parameters for Eq. (16), appropriate for the initial slopes in FNS19: (β, ζ) = (1, 0.5). The values are given in the rightmost column of Table 2.

The results are shown in Fig. 9. The green shaded region represents the range of numerical results found in FNS19. As is immediately evident from the figure, our high mass torque prescription does not work well in this case. Migration proceeds much slower, despite the reduced or inverted torque density (right panel). Mass growth does not stop due to the high mass reservoir available at large separations (left panel). The torque mod prescription fares much better. The mass evolution in this case is similar to that found by FNS19, though accretion stops later and more abruptly. Consequently the final mass reached in our simulation is approximately 15% higher than the highest value reached in FNS19 after 10kyr (data from FNS19 is only available until 10 kyr). Migration proceeds in line with the slowest code (GIZMO-MFM) applied in FNS19 during the first ~6 kyr. It slows down later, leading to a somewhat smaller separation after 10 kyr. This result is likely influenced by the different temperatures at small separations, as well as the higher planet mass in our simulation (see Sect. 3).

The difference in the migration tracks for our two torque models in this large separation case is noteworthy. In the other comparisons we performed the migration tracks were very similar. The reason for the discrepancy is that the surface density is unperturbed at the beginning of the simulation. In this scenario, the torque mod gives higher migration rates than the high mass torque since the amplitude of the latter is reduced (see Sect. 2). When starting with a low initial planet mass, the surface density near the planet is gradually reduced, and by the time the high mass torque is reduced in amplitude significantly, the corresponding region of the disc is already partly depleted and contributes little to the overall torque. We therefore expect the torque densities in the simulations of FNS19 to differ from our high mass torque. This cannot be confirmed since the torque densities are not provided in FNS19. Clearly, the high mass torque should not be used with massive planets in unperturbed discs. We note, however, that it may still be applicable in fragmenting discs if the fragment mass is removed from the disc (which is not done in FNS19). In this case the surface density would be reduced automatically in the relevant region.

The population synthesis models also shown in Fig. 9 exhibit a much larger spread in the final semi-major axis reached than what is seen in the hydrodynamic simulations. In the model of Forgan & Rice (2013) the planet opens a gap immediately, thus staying on a wide orbit. Conversely, the model from Müller et al. (2018) predicts very fast migration initially, and a gap is opened rather abruptly, when the planet reaches the inner disc. In order for a gap to open, the authors demand that the time to cross the gap is 100 times (solid line) or 1000 times (dotted line in Fig. 9) longer than the time it takes to open the gap. We note that Müller et al. (2018) use the migration timescale from Baruteau et al. (2011), which is valid in self-gravitating discs. The initially much faster migration rates are therefore expected. The result from the population synthesis code of Nayakshin & Fletcher (2015) is shown as the grey hatched region in Fig. 9. According to Sect. 4 of FNS19, the authors assumed that the viscosity α-parameter is a random variable taking values between 0.005 and 0.05. FNS19 show the two extreme values in their Fig. 9. Our Fig. 9 shows the entire range. The lower end of the grey shaded region in the right panel corresponds to α = 0.05. Since the viscosity is close to this value in FNS19, we expect the model from Nayakshin & Fletcher (2015) to agree quite well with the hydrodynamic calculations when using the same viscosity.

The mass accretion rate also differs substantially between the codes. It is typically much lower in the population synthesis codes compared to the hydrodynamic simulations. The exception is the model from Nayakshin & Fletcher (2015), when a high value for α is used (corresponding to the upper end of the hatched region in the left panel of Fig. 9). This result is in line with the results of FNS19. We note that accretion was neglected in Nayakshin & Fletcher (2015) (and in Forgan & Rice 2013), and is only added for the purpose of comparison in FNS19.

Fig. 9 Comparison of our model to results from the code comparison project in FNS19. Evolution of an initially 2 $M_{0_{+}}$ planet inserted at 120 au. Left: mass vs. time; right: separation vs. time. The green shaded region represents the region of parameter space covered by six different hydrodynamic codes. The results from our calculation are shown as thick dotted and dash-dotted lines (legend in the right panel). Results from different disc instability population synthesis codes are also shown (legend in the left panel). The grey shaded-hatched region corresponds to the region of parameter space covered by Nayakshin & Fletcher (2015). The orange dashed line shows the result from Forgan & Rice (2013). The result from Müller et al. (2018) is shown in blue. The dotted line (covered by the solid line in the left panel) represents a different gap opening criterion (see details in Sect. 4.8).

5 Comparison with published prescriptions

In Sects. 3 and 4 we describe our model and identify its advantages over existing approaches. We also demonstrate how different iterations of our prescription compare to 3D hydrodynamic simulations. In this section we compare our accretion model combined with the high mass torque to other existing prescriptions. We note that our model was specifically developed to be used in situations where existing prescriptions are inapplicable (Sect. 1). Nevertheless, it is interesting to compare different approaches in a regime where various published prescriptions are valid. For this comparison we adopted the model presented in Bitsch et al. (2015, BLJ15). The authors study accretion and migration of protoplanets in evolving discs in the pebble-based core accretion scenario. Since the accretion model presented here considers the disc-limited regime of gas accretion, we neglected the effects related to solid accretion in the comparison. For rapid gas accretion BLJ15 use the approach from Machida et al. (2010). Migration is described by the torque formula from Paardekooper et al. (2011), which includes non-isothermal and saturation effects. The discs in DL08 are locally isothermal. Therefore, in order to have a reasonable comparison, we applied the formula for the locally isothermal limit of Paardekooper et al. (2010). This means we are applying a modified version of the model in BLJ15, denoted BLJ15_li for locally isothermal. In order to account for the planetary growth, including the formation of a gap and the transition into the type II regime of migration, BLJ15 apply a reduction on the type I timescale based on the gap’s depth (Crida & Morbidelli 2007), a reduction of the type II timescale (Baruteau et al. 2014) and an additional smoothing between the two regimes. We use the same formulae and run the comparison for all the cases we studied in Sects. 4.1–4.5. The details can be found in Appendix E. We find that the model used by BLJ15 gives lower gas accretion rates because the underlying accretion model of Machida et al. (2010) already predicts lower accretion. In Fig. 10 we show the inverse growth timescale for the BLJ15_li model applied to the baseline case (Sect. 4.1). The predicted accretion rates are lower by up to a factor of two. At the end of the simulations, the resulting planetary masses are therefore lower than those from DL08. This is seen in the left panel of Fig. 11. It shows the time evolution of the planet’s mass in the baseline case discussed in Sect. 4.1. We find that gas accretion is significantly lower in the Bondi and the Hill regimes. The planet’s mass reaches ≈4 × 10⁻⁴ M_⊕, 60% less compared to DL08 after 1100 orbits. The lower planetary masses also lead to slower migration in the type I regime. As a result the planet remains further away from the star (right panel). Migration proceeds somewhat slower than seen in DL08, although the slope of the semi-major axis is similar (right panel) after migration slows down due to gap formation after around 800 orbits (i.e. fully in the type II regime). The main reason for the slower migration is the lower planetary mass at earlier times (the type I torque is proportional to M_p). We note that in the description of BLJ15, the surface density is not perturbed by the planet. For the comparison, we also mimicked this behaviour in our model: the accreted gas is not removed from the disc (no conservation of mass), and the disc does not feel the tidal interaction with the planet. Therefore, the planet continues accreting even if a gap is formed, and may eventually reach a higher mass than we see in the simulations, where accretion is modelled self-consistently.

Figure 12 shows the planet’s mass and semi-major axis as a function of time for the case with an increased initial surface density. Again a slower accretion leads to less migration compared to BLJ15_li, as seen in the left and right panels. Migration slows down due to gap formation at around 270 orbits. The comparison with the other two cases from Sect. 4 can be found in Appendix E. Given the differences seen among hydrodynamic simulations of gas accretion themselves (see e.g. Sect. 5 in Machida et al. 2010), the agreement between our model and that used in BLJ15 is still reasonable.

Fig. 10 Inverse growth timescale, as in the top left panel of Fig. 4. The brown dotted line shows the result from the locally isothermal version of the model in BLJ15.

6 Discussion

Our comparisons with DL08 show reasonable agreement in four out of five cases. The discrepancy in the high-viscosity case, where mass accretion is much lower in our model, is not surprising. The corresponding hydrodynamic simulation is dominated by global effects. Our result seems reasonable in this case as well. The agreement seen in the comparison with DL08 when using the high mass torque prescription indicates that applying torque densities obtained from hydrodynamic simulations in 1D models is a promising approach for modelling migration and accretion simultaneously, including gap formation.

This is attractive since angular momentum is conserved in this approach. Using the same torque densities for the high-mass wide-separation planets studied in FNS19 did not give a good agreement. The high mass torque should therefore not be used with initially massive planets in unperturbed discs. However, the torque mod prescription gives reasonable agreement even in this case.

A limitation of our model is that it was studied only in discs that are not self-gravitating. This is also true for all the hydrodynamic simulations we compared our results to. If discs are self-gravitating, additional effects will have to be taken into account. Migration may be much faster (Baruteau et al. 2011), and accretion may be significantly reduced (Shabram & Boley 2013). A way to model migration in a self-gravitating disc and still conserve angular momentum is suggested in Sect. 3 of Nayakshin (2015). It could be applied in our model as well.

We also did not investigate the long-term evolution of our systems, mainly because no data is available for comparison from DL08 and FNS19. Generally we expect a good agreement also on longer timescales both for accretion and migration. At the end of the simulations, accretion has almost ceased because the mass reservoir is depleted, and the slopes of our migration tracks agree well with those in DL08. Also, once a deep gap has opened, the exact shape of the torque density near the planet no longer has a strong influence on migration. There is one exception: once a deep gap has opened, mass flow through the gap is essentially stopped in 1D models. This is not seen in hydrodynamic models. On the contrary, significant mass transport is still observed even across deep gaps (Dürmann & Kley 2015, 2017; Zhu et al. 2011). A possible way to account for this is described in Lubow & D'Angelo (2006); D'Angelo & Bodenheimer (2016).

Our simulations also show that the surface density in the gap drops to extremely low values (≪c10⁻¹⁰g cm⁻²) at some point. This effect is not seen in hydrodynamic simulations. Forexample in Figs. 4, 5, B.1, and B.2, Σ_B/Σ₀ remains at a small (positive) value in DL08. This does not influence our results in terms of planet mass and semi-major axis, but this topic should be investigated further in the future.

Clearly, more work is needed to study the applicability of torque densities from complex simulations in 1D models. Future studies should focus on comparisons in a wide parameter space, realistic thermodynamics, and ideally also include the self-gravitating regime. Our work is a only first step in this direction.

Fig. 11 Same as Fig. 4, but for a comparison with BLJ15li. The brown dotted line represents the results when using their prescription for migration and accretion (see Appendix B.1). The results from our high mass torque model are shown for reference.

Fig. 12 Same as Fig. 11, but for increased surface density (300 g cm−2 at 5.2 au).

7 Summary, conclusions, and outlook

We developed a 1D prescription for migration and accretion of planets in discs and calibrated it with results from hydrodynamic simulations. Its performance was assessed with different initial surface densities, temperatures, viscosities and planet locations, comparing to the results of two different hydrodynamic simulations. The accretion model is based on the Bondi-Hill accretion given in DL08, refined and adapted to a 1D vertically integrated disc model. This model for accretion seems fairly robust, given the good agreement in our comparisons.

Migration is modelled by means of torque densities. We used the torque density formula from D'Angelo & Lubow (2010) and implemented a modification inspired by Hallam & Paardekooper (2017). In a second approach we corrected for the influence of the mass growth of the planet by directly applying torque densities for high-mass planets from D'Angelo & Lubow (2010). Our results from the two approaches are in good agreement with the hydrodynamic simulations of DL08. Our first modification also works well when applied to a 2 $M_{0_{+}}$ planet inserted at a large separation. We conclude the following:

• Our proposed model for gas accretion can be applied in 1D models for a range of different parameters.

• Our torque mod migration model gives reasonable predictions both for initially low-mass planets at small separation and initially massive planets on wide orbits. However, it does not conserve angular momentum, and the high mass torque model should thus be preferred.

• The high mass torque model we propose provides a description of migration that works well in the range of parameters we studied. Together with our accretion model, it describes accurately the mass and semi-major axis evolution of initially small planets that undergo runaway accretion and open a gap. It is self-consistent by conserving mass and angular momentum. We therefore recommend this combination for application in 1D models of planet formation.

Our model for gas accretion works reasonably well given the differences of at most 15% using the high mass torque model and 30% when applying the torque mod compared to DL08. The result of the torque mod model gives ~ 10–20% higher masses than the highest result in FNS19 for an initially massive planet on a wide orbit. Nevertheless, it is important to study this topic further. Gas accretion is a challenging topic to study, also in hydrodynamic simulations. It may be influenced by the formation of circumplanetary discs that require highresolution simulations to be described accurately (e.g. Zhu et al. 2012; Szulágyi et al. 2017; Oliva & Kuiper 2020). More such simulations, covering a wide parameter space, are necessary to improve our understanding of accretion. This is particularly important if the discs are self-gravitating. In this regime, accretion is thought to be inhibited significantly due to the truncation of the circumplanetary disc (Shabram & Boley 2013).

Migration is also a challenge for hydrodynamic models. A proper treatment of thermodynamics is key to predicting accurate migration rates (see e.g. Rowther & Meru 2020). Even so, only simulations applying simplified thermodynamics can run for long enough for a large number of initial conditions. Our results indicate that applying torque densities from hydrodynamic models in 1D codes is feasible, though more work is needed to test this in a wide parameter space. It is still unclear to what degree torque densities from hydrodynamic simulations that consider more elaborate thermodynamics can be applied in 1D models. In the future, torque densities should be measured in radiation-hydrodynamic simulations for different planet masses, separations, and viscosities. Their applicability in 1D models should be reassessed.

Future observations may give essential clues about the migration process. For instance, kinematical detections of exoplanets may help us understand the gap formation process further, since they can probe the gas dynamics in protoplanetary discs with high precision (Teague et al. 2019).

The prescription for migration based on the high mass torque and our accretion scheme can be easily implemented in 1D planet formation models. Details on the implementation are presented in Appendix F.

Our study demonstrated the complexity of the topic and it is clear that further work is needed. In the future, we plan to implement our prescriptions in population synthesis studies of planet formation. Our aim is to add a treatment of clump evolution to our population synthesis of discs (Schib et al. 2021). The prescriptions provided in this paper will be an important part of this. The model for migration and accretion derived here should be useful for 1D models in general. Being able to describe accretion, migration, and gap formation self-consistently, relying only on local disc parameters, can improve predictions from planet population synthesis models and our understating of planet formation.

Acknowledgements

We thank the anonymous referee for valuable comments. We also thank Gennaro D’Angelo and Clément Baruteau for the insightful discussions; and Bertram Bitsch, Aurélien Crida, Tom Hands and Lucio Mayer for helpful comments. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. O.S. and C.M. acknowledge the support from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. R.H. acknowledges support from SNSF grant 288828_188468.

Appendix A Evolution of the surface density

Here we present the evolution of the surface density for the cases not shown in Fig. 3. The top left panel shows the case of a higher initial surface density (Sect. 4.2, 300 g cm⁻² at 5.2 au). The gap opening proceeds qualitatively similarly to the baseline case, but much faster, as is expected from the faster increase of the planet’s mass. In the top right panel the very high surface density (Sect. 4.2, initially 500 g cm⁻² at 5.2 au) is shown. We note that the scaling of the y-axis is different from the other panels. Gap opening is even faster here due to the rapid increase in the planet mass. The bottom left panel shows the low-temperature case (Sect 4.4, H/r = 0.04). The evolution of the surface density is qualitatively very similar to the baseline case, though much faster due to the stronger torque at lower temperature. In the bottom right panel we finally show the high-viscosity case (Sect 4.5, 10¹⁶ cm² s⁻¹). Here the evolution of the disc is important, as discussed. The background surface density drops by more than a factor of two during the first 500 orbits. Due to the high viscosity, it takes a very long time for a gap to open. We note that the time interval shown here is much longer than in the simulation of DL08.

Fig. A.1 Evolution of the surface density for the cases not shown in Fig. 3. Top left: Higher surface density (initially 300 g cm−2 at 5.2 au). Top right: Very high surface density (initially 500 g cm−2 at 5.2 au). Bottom left: Lower temperature (H/r = 0.04). Bottom right: Higher viscosity (1016 cm2 s−1).

Appendix B Impulse approximation and type I migration

In this section we investigate how some other migration prescriptions used in the literature compare to our test cases. One of them is the classic impulse approximation (Lin & Papaloizou 1979a,b, 1986) $$\frac{d\mathcal{T}}{dm}=\text{sign}\left(r_p-r\right)\frac{f}{2}\text{Ω}{\left(r_p\right)}^2{r}_p^2{\left(\frac{M_p}{M_{*}}\right)}^2{\left(\frac{r}{\left|{\text{Δ}}_p\right|}\right)}^4,$$(B.1)

where the term |Δ_p| = max(H, |r − r_p|) excludes material closer than one scale height from the planet. The factor f is an order of unity constant taken to be 0.23 (Papaloizou & Lin 1984). The function sign(r_p − r) is replaced by (r_p − r)/H inside of one scale height to prevent a discontinuity at r = r_p (Lin & Papaloizou 1986). We call Eq. B.1 the LP86 formula in the following.

Another torque density commonly used is that given in Armitage et al. (2002): $$\frac{d\mathcal{T}}{dm}=\{\begin{array}{ll}-\frac{1}{2}\text{Ω}{\left(r_p\right)}^2{r}_p^2{\left(\frac{M_p}{M_{*}}\right)}^2{\left(\frac{r}{\left|{\text{Δ}}_p\right|}\right)}^4,\hfill & r<r_p\hfill \\ \frac{1}{2}\text{Ω}{\left(r_p\right)}^2{r}_p^2{\left(\frac{M_p}{M_{*}}\right)}^2{\left(\frac{r_p}{\left|{\text{Δ}}_p\right|}\right)}^4,\hfill & r>r_p.\hfill \end{array}$$(B.2)

It is similar to Eq. B.1, but it uses a modification to give a symmetric treatment for the disc inside and outside the planet’s orbit. We call it the A02 formula from now on. The LP86 formula and the A02 formula are shown in Fig. 2.

Recently, it has been proposed that type II migration is nothing other than type I migration that uses the reduced surface density in the gap (Kanagawa et al. 2018). The authors did not include gas accretion by the planet. We study a modification of this scenario, which we call ‘type I’, by applying a standard type I migration timescale τ_l (Tanaka et al. 2002) $${\tau }_I={\left(2.7+1.1\beta \right)}^{-1}\frac{M_{*}}{M_p}\frac{M_{*}}{{\text{Σ}}_p{r}_p^2}\left(\frac{c_{\text{s}}}{r_p{\text{Ω}}_p}\right){\text{Ω}}_p^{-1}$$(B.3)

without applying a torque to the disc. This violates the conservation of angular momentum, but it allows us to study an interesting effect. The surface density near the planet is still reduced due to accretion. The migration of the planet is slowed down exclusively through this reduction.

We do not discuss the case with increased viscosity again here. As noted in Sect. 4.5, this case is dominated by global effects. As for the earlier cases, the parameter C_H is chosen in such a way that the inverse growth timescale agrees with the analytic formula from DL08. In the case of the type I torque we also revert to using a feeding zone radius of 1R_H to prevent too high accretion rates. The numerical values chosen for these parameters are given in Table 1.

B.1 Baseline case

Figure B.1 shows the results for the three alternative migration treatments in the baseline case (Sect. 4.1). With the LP86 formula, the mass accretion (top left panel) is very similar to the result from DL08 and our earlier results. When we apply the A02 formula, the accretion rate drops earlier despite the slightly increased C_H (Table 1), leading to a lower final mass (top right panel). In the type I prescription the accretion rate drops later, which leads to a somewhat higher final mass. In this case we used a smaller R_f and a higher C_H in order to match the Hill accretion rate given in DL08. Clearly, in the absence of a tidal torque, less material is removed from the planet location. This is also seen in the bottom right panel: the surface density near the planet is depleted more slowly than seen in DL08. On the contrary, the depletion is much too fast when applying the LP86 formula, and even faster with the A02 formula. All three prescriptions lead to a mass accretion that is in reasonable agreement (within 25%) with that found by DL08. However, only the ‘type I’ method closely follows their migration track during the first ~ 500 orbits. The approximate agreement in the other two at later times may be a coincidence.

B.2 Higher surface density

Here we discuss the LP86 formula, the A02 formula, and our ‘type I’ approach in the case with an increased initial surface density. The results are given in Fig. B.2. The mass accretion (top panels) proceeds qualitatively similarly to the baseline case, with the A02 formula giving somewhat too little and ‘type I’ giving too high accretion. This time, the LP86 formula also gives higher accretion at late times than found in DL08. Migration (bottom left panel) is surprisingly similar to DL08 with the LP86 formula, though migration is too slow and gap formation too abrupt. The A02 formula shows a migration track that slows down too early and too abruptly. In this case, the ‘type I’ prescription shows almost no slowdown during the first ~ 250 orbits. A significant slowdown only occurs after ~ 400 orbits when the planet reaches less than 0.6 times its initial semi-major axis (not seen). The depletion of the surface density near the planet (bottom right panel) proceeds in an analogous fashion to the baseline case (B.1).

B.3 Very high surface density

In this section we discuss the alternative prescriptions in the case of an initial surface density of 500 g cm⁻¹. The results can be found in Fig. B.3. Mass accretion and depletion of the surface density near the planet once more proceed qualitatively similarly as above (top and bottom right panels). In this case none of the prescriptions give a migration track that agrees well with DL08. They either open a gap too far out or too far in, as seen in the bottom left panel.

B.4 Lower temperature

Here, we look at how the alternative torque models compare in the case of reduced temperature. The results can be found in Fig. B.4. Mass accretion proceeds again similarly (top panels). The difference between the torque models in this case are moderate. The surface density near the planet (bottom right panel) is affected differently again by the different torque densities, as described above. Due to the lower temperature, gap opening is much faster. This is most pronounced when using the A02 formula. In this case a gap starts to open immediately, leaving the planet too far out (bottom left panel). When using the ‘type I’ prescription the planet migrates in line with DL08, but slows down less, hence migrating much further in. The LP86 formula gives a reasonable agreement with DL08 for the migration track.

In summary, we find that none of the alternative prescriptions we studied in Appendix B gives an agreement to the results from DL08 similar to our torque mod and high mass torque prescriptions discussed in the main text. The LP86 formula and the A02 formula both open gaps too early in at least some of the cases, and deplete the surface density around the planet too much. The ‘type I’ recipe, on the other hand, gives good agreement for the semi-major axis at early times, but migration slows down much too late in some cases. This is related to the evolution of the surface density near the planet, which is depleted at a rate slower than seen in DL08. It is worth noting that this prescription still causes a reduction in the surface density near the planet by more than an order of magnitude through accretion alone, and consequently an eventual slowdown of migration.

Fig. B.1 Same as Fig. 4, but for alternative torque models. The solid brown line represents the results from the LP86 formula. The dashed pink line shows the results based on the A02 formula. The blue dotted line gives the results when applying only a type I torque using the surface density at the planet location.

Fig. B.2 Same as Fig. B.1, but for increased surface density (300 g cm−2 at 5.2 au).

Fig. B.3 Same as Fig. B.1, but for initial surface density of 500 g cm−2 at 5.2 au.

Fig. B.4 Same as Fig. B.1, but for reduced temperature (H/r = 0.04).

Appendix C Sensitivity of migration on surface density and temperature slopes

In Section 2.3 we argued that using torque densities based on the work of DL10 is a reasonable approach even if the slopes of surface density and temperature, β and ζ, vary slightly. Here we present a sensitivity calculation to support this statement. For this test we used the torque mod, since it can be applied to an arbitrary torque density covered by the grid from D'Angelo & Lubow (2010). Figure C.1 shows the time evolution of the planet’s semi-major axis when different functions ℱ (x,ß, ζ) are used. Either β or ζ is varied by 0.25 in both directions. The deviation from the nominal run with parameters (0.5, 1) is very small when the temperature slope ζ is varied. The change in semi-major axis (5.2 au minus r_p after 1300 orbits) is at most 4% different. The difference is larger, but still moderate when changing β (at most 12%). The p_i parameters used for the calculations shown in Fig. C.1 are listed in Table C.1 The p_i for these combinations of β and ζ are not directly available from DL10. Therefore, we interpolated the values for (0.75,1), (0.25,1), (0.5,1.25), and 0.5,0.75) by using the values for [(0.5,1), (1.5,1)], [(0.5,1), (0,1)], [(0.5,1), (0,2), (1,2)], and [(0,0), (1,0), (0.5,1)], respectively, given in their Table 1.

Fig. C.1 Evolution of the semi-major axis for different values of the slopes of surface density and temperature (see legend and the text in Appendix C).

Appendix D Accretion with autogravitation

Here we derive the expression for ${\dot{M}}_{\text{gas}}$ (see Sect. 2.2) when the disc’s autogravitation is taken into account. We follow the approach presented in Hueso & Guillot (2005). In this case equations 2 and 5 need to be replaced by $$\rho \left(r,z\right)={\rho }_0\left(r\right)\exp \left(-\left(\frac{\left|z\right|}{H_0}+{\left(\frac{z}{H_1}\right)}^2\right)\right),$$(D.1)

and $$H=\frac{H_1}{\sqrt{2}}\exp \left(\frac{H_1^2}{4H_0^2}\right)\left(1-\text{erf}\left(\frac{H_1}{2H_0}\right)\right),$$(D.2)

where $$H_0=\frac{c_{\text{s}}^2}{4\pi G\text{Σ}},$$(D.3) $$H_1=\frac{\sqrt{2}{c}_{\text{s}}}{\sqrt{\frac{G{M}_{*}}{r^3}}}.$$(D.4)

Plugging Eq. D.1 into Eq. 6 then yields $$\begin{array}{ll}{\dot{M}}_{\text{gas}}\hfill & =\sqrt{\pi }{C}_{\text{B,H}}{\displaystyle {\int }_{r_p-R_f}^{r_p+R_f}{\rho }_0\left(r\right)H_1\text{}\exp \left(\frac{H_1^2}{4H_0^2}\right)\times }\hfill \\ \hfill & \left(\text{erf}\left[\frac{\sqrt{R_f^2-{\left(r-r_p\right)}^2}}{H_1}+\frac{H_1}{2H_0}\right]-\text{erf}\left[\frac{H_1}{2H_0}\right]\right){\upsilon }_{\text{rel}}\text{d}r.\hfill \end{array}$$(D.5)

The expression for the angular frequency also changes when autogravitation is considered. For a more detailed description, see Section 2.2 in Schib et al. (2021).

Appendix E Comparison

In this section we present the results of our comparison with the migration-accretion prescription presented in BLJ15. The authors use the torque formula from Paardekooper et al. (2011) for type I migration. This formula was derived for non-isothermal type I migration. Since the simulations in DL08, on which our work is based, are locally isothermal, we used the torque formula for the locally isothermal isothermal limit given in Paardekooper et al. (2010) instead as it is more appropriate for our comparison: $${\tau }_{\text{lociso}}={\left(1.7+2\beta +1.8\zeta \right)}^{-1}\frac{M_{*}}{M_p}\frac{M_{*}}{{\Sigma }_p{r}_p^2}\left(\frac{c_{\text{s}}}{r_p{\text{Ω}}_p}\right){\text{Ω}}_p^{-1}.$$(E.1)

In the type II regime, migration is assumed to proceed on a viscous timescale ${\tau }_{\text{visc}}=r_{\text{p}}^2/v$ with a slowdown for massive planets (Alibert et al. 2005; Baruteau et al. 2014). This leads to a type II timescale of $${\tau }_{\text{II}}={\tau }_{\text{visc}}\times \max \left(1,\frac{M_p}{4\pi \text{Σ}{r}_p^2}\right),$$(E.2)

where the (unperturbed) surface density is evaluated at the planet’s location. In order to account for the slowdown in migration related to the formation of a gap, the migration rate is multiplied by the function f (Ƥ): $$f\left(p\right)=\{\begin{array}{ll}\frac{p-0.541}{4}\hfill & \text{if\hspace{0.17em}}p<2.4646\hfill \\ 1-\exp \left(-\frac{p^{3/4}}{3}\right)\hfill & \text{otherwise}.\hfill \end{array}$$(E.3)

Here the function Ƥ describes the gap opening and is defined as (Crida et al. 2006) $$p=\frac{3}{4}\frac{H}{R_{\text{H}}}+\frac{50}{qℛ},$$(E.4)

where q = M_p/M_* and $ℛ=r_{\text{p}}^2{\text{Ω}}_p/v$ the Reynolds number. A gap is assumed to be fully open (surface density reduced by a factor of ten) when Ƥ ≈ 1. In this case the type I migration timescale calculated as described here may still differ from τ_II. As a result, BLJ15 interpolate between the type I and type II timescales to ensure a smooth transition. They use a linear interpolation in P starting at Ƥ = 2.4646 (fully type I) going to Ƥ = 1 (fully type II, private communication with B. Bitsch). Rapid gas accretion is modelled using the fitting formula presented in Machida et al. (2010) that gives the accretion rate for two regimes, $${\dot{M}}_{\text{gas,low}}=0.83{\text{Ω}}_p\text{Σ}{H}^2{\left(\frac{R_{\text{H}}}{H}\right)}^{9/2}$$(E.5)

and $${\dot{M}}_{\text{gas,high}}=0.14{\text{Ω}}_p\text{Σ}{H}^2$$(E.6)

for low-mass and high-mass planets, respectively. The lower of the two is used as the effective accretion rate. In Equations E.5 and E.6, Σ and H are evaluated at the planet’s location. No gas is removed from the disc, so Σ is not perturbed. We implemented this model into our code and applied it to the systems described in Sect. 4.1 to 4.5.

The results for the comparison with the baseline case and with the case with higher surface density are given in Sect. 5. Here we discuss two additional cases. Figure. E.1 shows the planetary mass and semi-major axis for the case of the fivefold increase in initial surface density. Again, gas accretion is much lower than seen in the hydrodynamic simulation. Consequently migration is also slower. In this case we see no slowdown due to gap formation during the first 150 orbits.

Finally, figure E.2 shows the comparison with the low-temperature case, where gas accretion is also much lower. The behaviour of migration is somewhat different from the other cases. While it is also slower early due to the lower planetary mass, migration does not slow down significantly due to the (partial) opening of the gap. Therefore, the planet actually migrates further in during the first 900 orbits than it does in DL08. In summary, we find that our Bondi-Hill accretion model combined with the high mass torque gives a better agreement with the simulations of DL08. This is not surprising since we calibrated our accretion model with one of their simulations. We also applied the accretion rates from Eq. E.5 and E.6 to the test case from the code comparison project (Sect. 4.8). In this case it also leads to low accretion rates compared to all other models. The accretion rate is more than a factor of two lower than the lowest-mass case shown in the left panel of Fig. 9. We note that Machida et al. (2010) did not consider planets more massive than 1$M_{0_{+}}$, and therefore their fitting formula may not be applicable here. Nevertheless it is unclear whether the accretion rates given by Eqs. E.5 and E.6 or the values we use in our model are more realistic, as the models involved are very different (locally isothermal in DL08, beta-cooling in FNS19, and radiation-hydrodynamic in Machida et al. (2010)).

Fig. E.1 Same as Fig. 11, but with a five times higher initial surface density (500 g cm−2 at 5.2 au).

Fig. E.2 Same as Fig. 11, but for reduced temperature (H/r = 0.04).

Appendix F Practical use

In this section we describe how our prescription can be implemented in 1D models. We recommend using the high mass torque in combination with the Bondi-Hill accretion scheme, as given in Table 1: with a feeding zone radius of 3 RH, C_B = 10.0, and C_H = 0.19. The torque density for the high mass torque is obtained by an interpolation of the digitised torque densities given in D'Angelo & Lubow (2010). Some examples are shown in Fig. F.1 and the data can be downloaded from the CDS. The interpolation at the lowest masses is done by using ℱ(x, 0.5, 1) as zero mass. No extrapolation is perfomed above the highest mass. The high mass torque is reliable as long as the gradients of the surface density and temperature do not deviate from (β, ζ) = (0.5, 1) by more than a few 10%. The main limitation is related to the migration of low-mass planets (tens of M_⊕) in non-isothermal discs. The condition for non-isothermality is given by (Dittkrist et al. 2014) τ_cool ≫ τ_U–turn, where $${\tau }_{\text{cool}}=\frac{l_{\text{cool}}\rho C_{\text{V}}}{8\sigma T^3}\left(8\rho \kappa l_{\text{cool}}+\frac{1}{\rho \kappa l_{\text{cool}}}\right),$$(F.1) $${\tau }_{\text{U-turn}}=\frac{64x_{\text{s}}{h}^2}{9q{r}_p{\text{Ω}}_p}.$$(F.2)

In the above equations $$x_{\text{s}}=1.16r_p\sqrt{\frac{q}{h\sqrt{\gamma }}}$$(F.3)

is the width of the co-orbital region, q = M_p/M_*, h = H/r_p, l_cool = min(H, x_s) is a cooling length, $\rho =\sum /\left(\sqrt{2\pi }H\right)$ is the mid-plane density, γ is the adiabatic index, and C_V is the heat capacity at constant volume. In this regime non-isothermal effects may slow down or even invert migration, an effect that is not included in our prescription. As long as no torque densities in this regime are available, we recommend applying the torque calculated from the torque formula in Paardekooper et al. (2011) to the planet (while still applying the high mass torque density to the disc) as long as the departure from isothermality is strong. This violates conservation of angular momentum, but the effect at low masses is moderate.

In order to prevent a discontinuity in the accretion rate when transitioning from the Bondi regime to the Hill regime, we apply a smoothing in M_p. We set the effective accretion rate ${\dot{M}}_{\text{eff},i}$ in each grid cell to $${\dot{M}}_{\text{eff},i}={\dot{M}}_i\times \{\begin{array}{cc}1-0.25\left(0.5\tanh \left\{\frac{\left(M_p-M_t\right)}{M_t/3}\right\}+0.5\right)& M_p<M_t\\ 0.25+0.75\left(0.5\tanh \left\{\frac{\left(M_p-M_t\right)}{M_t/3}\right\}+0.5\right)& M_p\ge M_t\end{array}$$(F.4)

where ${\dot{M}}_i$ is the accretion rate in the i-th grid cell obtained from Eq. 12 and M_t is given in Eq. 10. We note that the function given in Eq. F.4 is not continuous at M_p = M_t. This is necessary to obtain a continuous global inverse growth timescale.

Fig. F.1 Torque densities in units of ${\text{Ω}}^2{r}_p^2{\left(M_p/M_{*}\right)}^2{\left(r_p/H\right)}^4$ for different planet masses taken from D'Angelo & Lubow (2010). The inversion developing as the mass increases is clearly seen near the origin. The high mass torque is obtained by an interpolation (see text Appendix F).

In this paper we did not investigate the regime where the disc is self-gravitating. Applying the torque densities we describe in this work for self-gravitating discs would require further modifications. This is also true for accretion, as we discuss in Sect. 6.

References

All Tables

All Figures

	Fig. 1 Schematic of the accretion geometry. The feeding zone is centred at the location of the planet rp, at z = 0. The direction tangential to the planet’s motion is perpendicular to the page.
In the text

Fig. 2 Torque densities used in this work: ℱ(0.5, 1) (DL10), a modified version (torque mod) where ℱ is truncated near the planet when applied to the disc, and a variant model (high mass torque) based on an interpolation of torque densities obtained from a hydrodynamical simulation. The last case is for a 1.5 M planet. The torque densities are scaled by ${\Omega }^2{r}_{\text{p}}^2$(Mp/M*)2(rp/H)4. Additionally, the LP86 formula (thin black solid line) and the A02 formula (thin black dashed line) are shown (see Appendix B). These are divided by 3 and 15, respectively for better visibility.

In the text

Fig. 3 Surface density around an accreting, migrating planet at different times. At the beginning, the planet has a low mass (5 M⊕) and does not perturb the surface density appreciably. As the planet grows, it changes the surface density in its vicinity both by accreting gas and by pushing disc material away through tidal interaction. After 1000 orbits, the planet has migrated to ≈4.2 au, reached a mass of 1 $M_{0_{+}}$ and opened a gap in the disc. The planet’s initial location (5.2 au) is indicated with a thin black vertical line. This simulation is described in more detail in Sect. 4.1.

In the text

Fig. 4 Time evolution of the baseline case with an initial surface density of 100 g cm−2 at 5.2 au. The solid black line represents the 3D result from DL08. The solid orange line corresponds to our 1D result when using the DL10 torque. The green dashed line gives again the result with the DL10 torque, but with the feeding zone radius increased threefold. The red dotted line shows a calculation where, additionally, the torque density has been modified. The purple dash-dotted line shows the variant model, where interpolated torque densities for planets of different masses from DL10 were applied. The torque densities are described in Sect. 2.3. Top left: inverse growth timescale in units of the orbital timescale at the initial location (≈12 yr). The thin black dashed line represents the fit given by Eqs. (7) and (8). Top right: mass evolution. Bottom left: orbital migration. The thin dashed black lines represent (from left to right) predictions for saturated and unsaturated type I, and type II theory. Bottom right: evolution of the mean surface density around the planet. Time is given in orbits at the planet’s initial separation (5.2 au).

In the text

	Fig. 5 Same as Fig. 4, but for increased surface density (300 g cm−2 at 5.2 au).
In the text

	Fig. 6 Same as Fig. 4, but with five times higher initial surface density (500 g cm−2 at 5.2 au).
In the text

	Fig. 7 Same as Fig. 4, but for reduced temperature (H/r = 0.04).
In the text

	Fig. 8 Same as Fig. 4, but for ten times higher viscosity (10 cm2 s−1).
In the text

Fig. 9 Comparison of our model to results from the code comparison project in FNS19. Evolution of an initially 2 $M_{0_{+}}$ planet inserted at 120 au. Left: mass vs. time; right: separation vs. time. The green shaded region represents the region of parameter space covered by six different hydrodynamic codes. The results from our calculation are shown as thick dotted and dash-dotted lines (legend in the right panel). Results from different disc instability population synthesis codes are also shown (legend in the left panel). The grey shaded-hatched region corresponds to the region of parameter space covered by Nayakshin & Fletcher (2015). The orange dashed line shows the result from Forgan & Rice (2013). The result from Müller et al. (2018) is shown in blue. The dotted line (covered by the solid line in the left panel) represents a different gap opening criterion (see details in Sect. 4.8).

In the text

	Fig. 10 Inverse growth timescale, as in the top left panel of Fig. 4. The brown dotted line shows the result from the locally isothermal version of the model in BLJ15.
In the text

	Fig. 11 Same as Fig. 4, but for a comparison with BLJ15li. The brown dotted line represents the results when using their prescription for migration and accretion (see Appendix B.1). The results from our high mass torque model are shown for reference.
In the text

	Fig. 12 Same as Fig. 11, but for increased surface density (300 g cm−2 at 5.2 au).
In the text

	Fig. A.1 Evolution of the surface density for the cases not shown in Fig. 3. Top left: Higher surface density (initially 300 g cm−2 at 5.2 au). Top right: Very high surface density (initially 500 g cm−2 at 5.2 au). Bottom left: Lower temperature (H/r = 0.04). Bottom right: Higher viscosity (1016 cm2 s−1).
In the text

	Fig. B.1 Same as Fig. 4, but for alternative torque models. The solid brown line represents the results from the LP86 formula. The dashed pink line shows the results based on the A02 formula. The blue dotted line gives the results when applying only a type I torque using the surface density at the planet location.
In the text

	Fig. B.2 Same as Fig. B.1, but for increased surface density (300 g cm−2 at 5.2 au).
In the text

	Fig. B.3 Same as Fig. B.1, but for initial surface density of 500 g cm−2 at 5.2 au.
In the text

	Fig. B.4 Same as Fig. B.1, but for reduced temperature (H/r = 0.04).
In the text

	Fig. C.1 Evolution of the semi-major axis for different values of the slopes of surface density and temperature (see legend and the text in Appendix C).
In the text

	Fig. E.1 Same as Fig. 11, but with a five times higher initial surface density (500 g cm−2 at 5.2 au).
In the text

	Fig. E.2 Same as Fig. 11, but for reduced temperature (H/r = 0.04).
In the text

	Fig. F.1 Torque densities in units of ${\text{Ω}}^2{r}_p^2{\left(M_p/M_{*}\right)}^2{\left(r_p/H\right)}^4$ for different planet masses taken from D'Angelo & Lubow (2010). The inversion developing as the mass increases is clearly seen near the origin. The high mass torque is obtained by an interpolation (see text Appendix F).
In the text