© The Authors 2023

1 Introduction

Observational surveys have demonstrated that the exoplanet population is dominated by super-Earths and mini-Neptunes, while giant planets make up a sub-population that accounts for only ~10% of the overall planet number (Batalha et al. 2013; Mayor et al. 2011). Efforts have been made to construct an unbiased sample of giant planets based on radial velocity discoveries, and these suggest that the distribution of orbital periods has two maxima, separated by a valley that sits in the range 10 ≤ P ≤ 100 days (Cumming et al. 2008; Santerne et al. 2016). A significant majority of giant planets orbit with periods longer than 100 days, and we refer to these as cold Jupiters. The masses of the giant planets discovered by radial velocities show a fairly flat distribution between 0.5 ≤ M_p sin i ≤ 2 M_Jup, which then declines down to ~4 M_Jup before flattening off at larger masses¹. Currently, we do not have a good theoretical understanding of the mass and orbital period distributions of giant exoplanets.

The giant planet sub-population can be considered to be distinct from the broader population of exoplanets for two physically motivated reasons. First, the core-accretion model shows that Jovian mass planets accrete the majority of their gas during a runaway phase of growth, whereas mini-Neptunes and Neptunelike planets accrete their gas envelopes quasi-statically (Pollack et al. 1996). This means that the rate at which gas is supplied to the planet by the protoplanetary disc is likely to be an important factor in determining the masses of giant planets. Second, gas giant planets are believed to open deep gaps in protoplanetary discs (Lin & Papaloizou 1979, 1986a; Goldreich & Tremaine 1980), and this can strongly influence both the gas accretion and the orbital migration rates of these bodies (Bryden et al. 1999; Kley 1999).

Giant planets are believed to start forming close to the snow line, where the condensation of water enhances the density of solid material and enables the rapid formation of the first planetesimals (Ida & Guillot 2016; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017). The location of the snow line in young discs is expected to be at 5–7 au (Lee et al. 2021; van’t Hoff et al. 2020), indicating that the cold Jupiter population migrated only a few astronomical units in the protoplanetary disc, which is at odds with the predictions of most studies of disc-driven giant planet migration. Giant planets are expected to form a gap in the disc and undergo type II migration. In the traditional picture, this mode of migration occurs with the planet migrating within the gap at the viscous evolution rate of the disc. Accretion rates onto young stars have typical values ṁ ~ 10⁻⁸ M_⊙ yr⁻¹ (Hartmann et al. 1998), which requires the viscosity parameter to have values α ≳ 10⁻³, giving rise to migration timescales ~0.1 Myr (Nelson et al. 2000). Given that discs have typical lifetimes of a few million years. (Haisch et al. 2001) this makes it difficult to explain the cold Jupiter population, and indeed models show that a giant planet needs to start undergoing type II migration at ~ 15–20 au from the star in order to halt at a distance of a few astronomical units after a few million years of migration (Coleman & Nelson 2014; Bitsch et al. 2015). Although the classical view of type II migration has been revisited recently (Dürmann & Kley 2015; Robert et al. 2018), such that the migration speed is not exactly equal to the unperturbed drift speed of the disc, it is still proportional to the viscosity when α ≥ 10⁻⁴, and a significant reduction in α is required to reduce the distance over which type II migration occurs over a disc’s lifetime (Lega et al. 2021).

Despite the presence of a gap, numerous studies of gas accretion onto Jovian mass planets show that a mass flux through the gap of ṁ ~ 10 M_Jup Myr⁻¹ is sustained for α ≥ 10⁻³, essentially corresponding to the unperturbed accretion flow through the disc (Bryden et al. 1999; Kley 1999; Lubow et al. 1999). Hence, if the planet can accrete a large fraction of the gas that is supplied to it, then it is difficult to explain the giant planet mass distribution, since the mass doubling time of a Jupiter-mass body is only ~10⁵ yr. Recent studies have shown that thermo-dynamic effects in the planet’s Hill sphere may be important for slowing down gas accretion under some circumstances (e.g. Szulágyi et al. 2016; Moldenhauer et al. 2021, 2022); however, the radiation-hydrodynamic simulations of Lambrechts et al. (2019) indicate that Jovian mass planets can accrete at a rate M_p ≥ 10 M_Jup Myr⁻¹. It is possible that other effects such as magnetic fields may influence the accretion rate, but the nonideal magnetohydrodynamical (MHD) study of gas accretion onto giant planets by Gressel et al. (2013) showed accretion rates in line with the values quoted above. Hence, within the context of viscous disc models, it is difficult to explain either the mass or the period distribution of the extrasolar giant planets.

The above discussion applies to planets embedded in viscous discs, where it was originally believed that the viscosity in discs might arise from the magneto-rotational instability (MRI, Balbus & Hawley 1991). The very low ionisation fraction in the main body of a protoplanet disc, however, quenches the MRI (Gammie 1996) and instead a combination of global magnetic fields, nonideal MHD effects and the ionisation of disc surface layers can launch magnetically driven disc winds that also drive a laminar accretion flow through the surface layers of the disc (Bai & Stone 2013; Gressel et al. 2015; Béthune et al. 2017). Hence, planet formation may occur in discs with very low levels of turbulence and where accretion towards the star occurs in narrow surface layers that sustain rapid radial gas flows.

Motivated by this change in our understanding of the internal dynamics of protoplanetary discs, we have recently presented simulations of giant planets migrating in very low viscosity discs. In Lega et al. (2021), we examined migration in discs without wind-driven accretion flows, and showed that initially a giant planet migrates inwards because a vortex forms at the outer gap edge, but once this vortex dissipates the ‘vortex-driven migration’ ceases and migration essentially halts. In Lega et al. (2022), we examined the migration of giant planets in discs with wind-driven laminar accretion flows, and we showed that the behaviour crucially depends on whether or not the torque from the planet is strong enough to block the accretion flow. When the flow is fast and occurs relatively unimpeded by the planet, we found that migration is slow and occurs at a speed ${\dot{r}}_{\text{p}}~3$ au Myr⁻¹. When flow is blocked by the planet then migration can be fast, occurring at a speed ${\dot{r}}_{\text{p}}~15$ au Myr⁻¹. Hence, for appropriate disc parameters, giant planets undergoing slow migration can explain the cold Jupiter population.

In this study we examine the rate of gas accretion onto giant planets that are held on fixed circular orbits, and which are embedded in layered disc models with laminar accretion flows that are very similar to those considered in Lega et al. (2022). Our aim is to determine the conditions under which the gas flow towards the planet can be significantly decreased by the planet’s tidal torques, such that the planet regulates the rate at which it can accrete gas, and to understand how the joint mass and orbital evolution of a gas giant planet proceeds in layered disc models.

The paper is organised as follows. In Sect. 2 we present the basic equations, physical model and numerical methods. In Sect. 3 we discuss our theoretical expectations, and in Sect. 4 we present the simulation results. In Sect. 5 we discuss the results and their implications for the formation and early evolution of giant plates, and in Sect. 6 we draw our conclusions.

2 Basic equations and numerical methods

2.1 Equations of motion

We use a variety of coordinate systems in this paper. Cartesian coordinates are denoted by (x, y, z), spherical polar coordinates by (r, θ, ϕ) and cylindrical polar coordinates by (R, ϕ, z). We solve the continuity, momentum and internal energy equations $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0\hfill \\ \frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)v=-\frac{1}{\rho }\nabla P-\nabla {\text{Φ+S}}_{\text{visc}}+f_{\text{wind}}\widehat{\varphi }\hfill \\ \frac{\partial e}{\partial t}+\nabla \cdot \left(ev\right)=-P\left(\nabla \cdot v\right)-\rho c\text{v}\frac{\left(T-T{}_{\text{i}}\right)}{{\tau }_{\text{cool}}},\hfill \end{array}$$(1)

where v is the velocity, ρ is the density, P is the pressure, e is the internal energy per unit volume and T is the temperature (T_i being its initial value and τ_cool being the cooling timescale on which temperature fluctuations relax towards T_i). The viscous force per unit mass is denoted by S_visc, and $f_{\text{wind}}\widehat{\varphi }$ is an azimuthal acceleration applied to the disc designed to mimic the back reaction on the disc when a magnetised wind is centrifugally launched from the disc surface ($\widehat{\varphi }$ being the unit vector in the ϕ direction). In this work we adopt an ideal gas equation of state, P = (γ − 1)e with γ = 7/5, and we set τ_cool = 1, in units of the local orbit period, as this timescale is long enough for the vertical shear instability (VSI) to be suppressed (Nelson et al. 2013). Recent radiation-hydrodynamic simulations that adopt realistic opacities indicate that the VSI does not operate in the inner regions (≲10 au) of protoplanetary discs (Flock et al. 2017), hence our choice for τ_cool is justified as this is the region of interest in this study.

The gravitational potential, Φ, arises from the central star and the planet (which is maintained on a fixed circular orbit), and is given by $$\text{Φ}\left(r\right)=-\frac{G{M}_{\ast }}{r}+{\text{Φ}}_{\text{p}}+{\text{Φ}}_{\text{ind}},$$(2)

where Φ_ind is an indirect term that arises because we work in a non-inertial frame centred on the star, and which includes contributions from the planet and the disc. Φ_p is the potential due to the planet, and we use the prescription of Kley et al. (2009) in our implementation of this term: $${\text{Φ}}_{\text{p}}=\{\begin{array}{ll}-\frac{G{M}_{\text{p}}}{d}\hfill & d>ϵ\hfill \\ -\frac{G{M}_{\text{p}}}{d}f\left(\frac{d}{ϵ}\right)\hfill & d\le ϵ\hfill \end{array}$$(3)

where d = |r – r_p| and f (x) = x⁴ – 2x³ + 2x. We adopt a softening length ϵ = 0.4R_H, where R_H = r_p (M_p/3M_*)^1/3 is the Hill radius.

2.2 Numerical methods

The above set of equations are solved using the 3-dimensional hydrodynamical code NIRVANA (Ziegler & Yorke 1997) in spherical polar coordinates (r, θ, ϕ). This code has been used extensively in the study of protoplanetary discs, with and without embedded planets (e.g. Nelson & Papaloizou 2004; Fromang & Nelson 2006; Nelson et al. 2013). NIRVANA uses an algorithm similar to the ZEUS code (Stone & Norman 1992). The equations are divided into source and transport terms, and operator splitting is used to update the state variables according to the different terms in this formalism. The van Leer (1977) upwind, monotonic advection scheme is employed to update the transport terms. We have also utilised the FARGO3D code (Benítez-Llambay & Masset 2016) to conduct simulations to compare with and validate the NIRVANA runs, and have obtained very similar results using the two codes. A comparison of the two codes is presented in Appendix A.

As discussed below, we employ disc models that have initial temperature distributions that give rise to constant aspect ratios, h ≡ H/r, where H is the vertical scale height. The computational domain in the meridional direction extends from the disc surface, located at 4 scale heights above the midplane, down to the disc midplane at θ = π/2. A symmetry boundary condition is employed at the midplane, and standard outflow conditions are applied at the disc surface. The azimuthal domain lies in the range 0 ≤ ϕ ≤ 2π. The radial domain covers the interval 1.04 ≤ r ≤ 46.8, where the unit of length is the astronomical unit (au). Further details about the system of units employed in this paper are given below. The planet is inserted at a radius of r_p = 5.2 au in all simulations, and remains on a fixed circular orbit. We employ wave damping boundary conditions near both the inner and outer radial boundaries, which are implemented by damping the density, radial and vertical velocities towards their initial values on a timescale of 0.1 local orbital periods (for more details see de Val-Borro et al. 2006). In addition, a linear bulk viscosity is employed near the outer boundary to provide additional damping. The implementation of this term is the same as in the ZEUS code (Stone & Norman 1992). In our runs the linear viscosity coefficient, Q_lin, ramps up linearly from 0 to 1 between the inner edge of the damping zone and the outer edge of the computational domain. Wave-damping is employed in the region r ≤ 1.612 au, and both wave-damping and linear viscosity are applied in the region r ≥ 35.672 au. The numerical resolution we adopt in all simulations is (N_r, N_θ, N_ϕ) = (576, 48, 384), giving a uniform radial grid spacing Δr ≃ 0.08.

The simulations are conducted in a reference frame that corotates with the planet. We initialise the disc models according to the physical parameters described below, and allow a disc to relax for more than 200 orbits at the location of the planet before the latter is inserted, in order to allow a steady mass flow to be established. Upon insertion, the planet-star mass ratio q_p = M_p/M_* = 6 × 10⁻⁵. Accretion of gas onto the planet occurs according to $$\frac{\text{d}{M}_{\text{p}}}{\text{d}t}={\displaystyle \sum _{i}K{\rho }_i\text{Δ}{V}_i,}$$(4)

where K is a constant (with 0 < Kdt < 1) and ρ_i ΔV_i is the mass contained in a grid cell, and the summation occurs over all cells that are located within a distance equal to half the planet’s Hill sphere radius. The accreted mass is removed from the disc and is added to the planet’s mass until q_p = M_p/M_* = 10⁻³ (corresponding approximately to the Jovian mass), after which time the planet’s mass is held fixed while the rate of gas accretion continues to be monitored. This procedure for introducing and growing the planet is adopted to ensure the disc is not shocked by the sudden introduction of a giant planet, to conserve mass during the planetary growth and gap formation phase, and also as a means of crudely mimicking the runaway gas accretion that is expected to occur onto a planet with q_p = 6 × 10⁻⁵ as it grows to become a gas giant. We initiate the simulations using the value K = 6 in Eq. (4), and at late times this is increased to K = 30. This was done because in earlier runs, adopting K = 6 did not lead to a clear separation of the long-term accretion rates onto the planets for different disc models, and in particular did not reflect the rate at which mass is supplied to the planet by the background disc. This problem is ameliorated by increasing the value of K, and we found that increasing the value to K = 30 led to a clear separation of the accretion rates experienced by the planets in the different runs, and to values that are in line with expectations given the rate of mass supply to the planet. Hence, we did not experiment with further values of K to examine the effects of varying this somewhat arbitrary parameter. Understanding what value of K is ‘realistic’ remains an area of active research (see Paardekooper et al. 2022, for example), and a discussion about related issues is provided in Sect. 5.3.

We note that our simulations do not resolve the structure of the gas flow within the planet’s Hill sphere because R_H/Δr ≃ 4.6 when q_p = 10⁻³. Hence, the simulations are designed to test the ability of the planet’s tidal torques to retard the large scale accretion flow of the gas into the gap and then into the Hill sphere, and do not address the question of how the gas then accretes onto the actual planet (see Sect. 5.3 for a discussion about gas supply to the Hill sphere and the planet’s accretion rate), which at this stage of the evolution is expected to have contracted to a size that is significantly smaller than the Hill sphere radius (e.g. Papaloizou & Nelson 2005). We further note that because we only simulate the upper hemisphere of the disc models, when quoting accreted masses, accretion rates onto the planet, or when adding mass to the planet due to accretion, we multiply the accreted mass by a factor of two to account for the mass that would be accreted from the missing hemisphere.

2.3 Disc models with spatially constant accretion rates

We adopt disc models for which the surface density, defined by $$\text{Σ=}{\displaystyle {\int }_{-\infty }^{\infty }\rho \text{\hspace{0.17em}d}z,}$$(5)

is initially axisymmetric. The midplane density is given by a power law in cylindrical radius $$\rho \left(R,z=0\right)={\rho }_0{\left(\frac{R}{R_0}\right)}^p,$$(6)

and the temperature is given by $$T\left(R,z\right)=T_0{\left(\frac{R}{R_0}\right)}^q.$$(7)

Expressions for the equilibrium densities and azimuthal velocities for these power-law disc models are given in Nelson et al. (2013). The isothermal sound speed $c_{\text{s}}=\sqrt{RT/\mu }$, and the pressure scale height of the disc H = c_s/Ω_K, where Ω_K is the keplerian angular velocity. We adopt the values p = −3/2 and q = −1 for the power-law indices in Eqs. (6) and (7), and a value of T₀ that gives an aspect ratio h = 0.05 throughout the disc. The value of ρ₀ is chosen so that $\text{Σ}=\sqrt{2\pi }{\rho }_{z=0}H=220$ g cm⁻² at r = 5.2 au (note that this value includes both the top and bottom hemispheres of the disc; the column density from the surface to the midplane is 110 g cm⁻²). We further note the surface density at 5 au in the minimum mass solar nebula has a value ~150 g cm⁻² (Hayashi 1981).

Test simulations have shown that when completely inviscid discs are considered, convergence between simulations performed at different resolutions can be unsatisfactory. Adding a small viscosity significantly improves the situation (Lega et al. 2022), so we include an alpha viscosity with α = 10⁻⁵ in our layered disc models and in a calibration run which does not contain an accreting layer (the No Wind simulation).

2.3.1 Layered accretion due to an external torque

To obtain an expression for the azimuthal acceleration, f_wind (first introduced in Eq. (1)), that drives a constant mass flux in a laminar disc, we first note that the mass flux at each cylindrical radius, R, is given by $$\dot{m}=2\pi R{\text{Σ}}_{\text{A}}\left|{\upsilon }_R\right|,$$(8)

where Σ_A is the column density associated only with the actively accreting layer of the disc. When considering layered accretion, Σ_A < Σ holds except when the whole vertical column of the disc is active and accreting. The specific angular momentum in a keplerian disc $j=\sqrt{G{M}_{*}R}$, and differentiating with respect to time gives $${\upsilon }_R=2\sqrt{\frac{R}{G{M}_{\ast }}}\frac{\partial j}{\partial t}.$$(9)

Combining Eqs. (8) and (9) gives $${f^{\prime }}_{\text{wind}}\widehat{\varphi }=-\sqrt{\frac{G{M}_{\ast }}{R^5}}\frac{\dot{m}}{4\pi {\text{Σ}}_{\text{A}}},$$(10)

where $\partial j/\partial t=R{f^{\prime }}_{\text{wind}}\widehat{\varphi }$. Equation (10) can be used to set up a disc with a chosen spatially and temporally constant radial mass flux, where the mass flux is confined to the surface layers where Σ(z) ≤ Σ_A, where Σ(z) is the column density measured from some height in the disc, z, up to the disc surface $$\text{Σ}\left(z\right)={\displaystyle {\int }_z^{\infty }\rho \text{d}z.}$$(11)

We refer to the torque that arises from ${f^{\prime }}_{\text{wind}}$ as the ‘wind-driven torque’, in recognition of the fact that this is supposed to mimic the torque that arises from the launching of a magnetised wind.

Our numerical implementation of Eq. (10) involves calculating the column density, Σ, at each position in the disc every 100 time steps, and then applying the following acceleration, which employs the Fermi function to transition the applied acceleration across the boundary between the accreting and non-accreting layers of the disc: $$f_{\text{wind}}={f^{\prime }}_{\text{wind}}\cdot \left(\frac{1}{\exp \left(\frac{\text{Σ}\left(z\right)-{\text{Σ}}_{\text{A}}}{{\text{Σ}}_{\text{T}}}\right)+1}\right),$$(12)

where Σ_T ≪ Σ_A sets the width of the transition between the active and dead layers. A further consideration is how mass is to be supplied at the disc outer radius so that it can maintain a constant mass flux over time. The use of a buffer zone at the outer disc edge, in which the density is continuously relaxed towards its initial value on sub-orbital timescales, ensures that our disc models maintain constant accretion rates through the main body of the disc at all times. As a final comment concerning the implementation of the wind-driven torque, we note that the effects of mass loss from the disc through the magnetocentrifugally driven wind are not included. Generally, the mass loss rate through the wind is significantly smaller than the mass flow rate in the accretion flow, and the run times we consider are too short for the mass loss in the wind to affect the global disc structure.

Equation (10) shows that for a given radial mass flux through the disc, ṁ, the applied acceleration is inversely proportional to the column density of the active layer, Σ_A. In this work we consider a single value of ṁ = 10⁻⁸ M_⊙ yr⁻¹ (each hemisphere of the disc provides half of this value) and values of Σ_A = 0.1, 1 and 10 g cm⁻². We adopt Σ_T = 0.1 for the Σ_A = 1 and 10 g cm⁻² cases, and Σ_T = 0.01 for the Σ_A = 0.1 g cm⁻² case. The steady state velocity profiles that arise from the application of Eqs. (10) and (12) are illustrated by Fig. 1, and as expected the radial velocity scales inversely with Σ_A. For reference, we note that at r = 5.2 au the isothermal sound speed c_s ≃ 0.65 km s⁻¹ and the Mach number of the radial accretion flow has values $ℳ\simeq {10}^{-1}$, 10⁻² and 10⁻³ for Σ_A = 0.1, 1 and 10 g cm⁻², respectively.

Figure 2 shows the azimuthally integrated mass fluxes corresponding to the velocities shown in Fig. 1, and vertically integrated radial profiles of the mass fluxes for the layered and viscous models (described below) are shown in Fig. 3, where it can be seen that both hemispheres of the disc combined produce a radial mass flux of 10⁻⁸ M_⊙ yr⁻¹. Hence, mass will be supplied to embedded planets by the large scale background accretion flow in the unperturbed discs at a rate of 10⁻⁸ M_⊙ yr⁻¹, equivalent to 10 M_Jup Myr⁻¹. The mass doubling time for a Jupiter mass planet being fed by such a flow would then be 10⁵ yr, similar to the values obtained by numerous previous simulations of planets embedded in viscous discs (e.g. Bryden et al. 1999; Kley 1999; Lubow et al. 1999). We note that the oscillations in the mass flux for the Σ_A = 10 case in Figs. 2 and 3 arise because the width adopted for the transition between the active and dead zones is only marginally resolved in this run. Although this is unlikely to adversely impact the results of the simulations presented here, we have confirmed that widening the transition removes these oscillations.

2.3.2 Viscous models

In addition to the externally torqued laminar models discussed in the previous section, for comparison purposes we also compute a viscous disc model that has the same radial mass flux. From standard viscous thin disc theory, we have the following expression for the torque per unit mass acting at each radius in the disc $${\upsilon }_R\frac{\partial }{\partial R}\left(R^2\text{Ω}\right)=\frac{1}{R\text{Σ}}\frac{\partial }{\partial R}\left(R^{3}v\text{Σ}\frac{\text{dΩ}}{\text{d}R}\right),$$(13)

where v is the kinematic viscosity. Given ṁ = −2πRΣv_R, we can write $$-\frac{\dot{m}}{2\pi }={\left(\frac{\partial }{\partial R}\left(R^2\text{Ω}\right)\right)}^{-1}\frac{\partial }{\partial R}\left(R^{3}v\text{Σ}\frac{\text{dΩ}}{\text{d}R}\right).$$(14)

If we adopt power-laws for the surface density and kinematic viscosity profiles, Σ(R) = Σ₀R^β and v(R) = v₀R^δ, and consider a keplerian disc with $\text{Ω}=\sqrt{G{M}_{*}/R^3}$, then we obtain the result that m is independent of radius when δ = −β. Hence, for the disc models with Σ ∝ R^−1/2 considered in this paper, we require δ = 1/2. Equation (14) also leads to the expression $$\dot{m}=3\pi v_0{\text{Σ}}_0,$$(15)

such that the kinematic viscosity needed to produce the required mass accretion rate can be determined once the surface density has been specified. The steady radial mass flux obtained from the viscous disc model we compute is shown in Fig. 3, where the magenta line is seen to sit under the lines for the other models at a value of ṁ = 10⁻⁸ M_⊙ yr⁻¹, as required.

Finally, we note in passing that the requirement for δ = 1/2 in a disc with Σ = Σ₀R^−1/2 is satisfied by the α model for the kinematic viscosity, v = αH²Ω, when H/r is constant, as is the case in our models.

Fig. 1 Contour plots showing the radial velocity distributions in the various layered disc models prior to insertion of the planets. The black dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. We note that only the upper hemisphere of the disc is simulated.

Fig. 2 Contour plots showing time-averaged and azimuthally-integrated radial mass flux distributions in the various layered disc models. The cyan dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. Note that only the upper hemisphere of the disc is simulated.

Fig. 3 Radial profiles of the gas accretion rates in the viscous and externally torqued disc models.

2.4 Summary of the runs

We present the results of five simulations in this paper: a No Wind simulation that does not have a laminar accretion flow and where a very small viscosity (α = 10⁻⁵) is employed; a Viscous run in which the radial mass flux ṁ = 10⁻⁸ M_⊙ yr⁻¹ at all radii; three layered models that contain laminar accretion flows with ṁ = 10⁻⁸ M_⊙ yr⁻¹ at all radii, and with Σ_A = 0.1, 1 and 10 g cm⁻².

2.5 Units

When discussing simulation results the unit of length is 1 au, the unit of time is years and the unit of mass is the solar mass or Jupiter’s mass (depending on the context, and where M_Jup = 10⁻³M_⊙). Surface densities, volume densities and velocities are generally quoted in cgs units.

3 Theoretical expectations

The long term steady gas accretion rate onto a planet embedded in the disc models should be determined by the balance between the tidal torque exerted on the disc by the planet, and the viscous or wind-driven torque acting to drive the radial mass flux through the disc. The viscous torque per unit mass is given by Eq. (13) and combining this with Λ_p, the torque per unit mass due to the planet, leads to the following expression for the radial mass flux in a viscous disc containing an embedded planet $$-\frac{\dot{m}}{2\pi }={\left(\frac{\partial }{\partial R}\left(R^2\text{Ω}\right)\right)}^{-1}\left[\frac{\partial }{\partial R}\left(R^{3}v\text{Σ}\frac{\text{dΩ}}{\text{d}R}\right)+R{\text{ΣΛ}}_{\text{p}}\right].$$(16)

An embedded planet can open a gap in the disc when tidal torques locally exceed the viscous torques, and Eq. (16) shows that the viscous torque readjusts to the changing disc structure such that in principle viscous and planet torques can cancel each other, resulting in ṁ = 0, with gas accretion onto the planet being switched off. The torque from the planet is highly localised, however, so the mass flux far from the planet is relatively unaffected and hence the gap structure must evolve on longer timescales as gas flows towards the planet, such that torque balance cannot be maintained. For a planet on a fixed circular orbit, a steady state then corresponds to mass accreting through the gap towards the planet at the rate supplied through the disc, part of which will be accreted by the planet. Indeed, numerous hydrodynamical simulations of planets embedded in viscous discs have shown that accretion through gaps is maintained at essentially the rate at which gas flows through the unperturbed disc (Bryden et al. 1999; Kley 1999; Lubow et al. 1999). Hence, in our viscous model we expect to see steady accretion onto the planet at a rate of ~10⁻⁸ M_⊙ yr⁻¹ (corresponding to 10 M_Jup Myr⁻¹).

We now consider the competition between wind-driven and planet torques in our layered disc models for a planet on a fixed circular orbit. We presented a similar discussion in Lega et al. (2022) that differs in detail with the discussion presented below, but which reached very similar conclusions. For the layered models the torque balance only needs to be considered in the active layer, and the equivalent expression to Eq. (16) can be written as $$-\frac{\dot{m}}{2\pi }={\left(\frac{\partial }{\partial R}\left(R^2\text{Ω}\right)\right)}^{-1}\left({\text{Λ}}_{\text{p}}+{\text{Γ}}_{\text{wind}}\right)R{\text{Σ}}_{\text{A}},$$(17)

where Γ_wind is the wind-driven torque per unit mass corresponding to Eq. (10): $${\text{Γ}}_{\text{wind}}=-\sqrt{\frac{G{M}_{\ast }}{R^3}}\frac{\dot{m}}{4\pi {\text{Σ}}_{\text{A}}}.$$(18)

We expect the planet to form a gap in the dead zone, and the radial mass flux passing through the gap and approaching the planet will depend on how Λ_p compares with Γ_wind, noting that in our model the wind-driven torque applied within the active column is constant and does not change with time for a given value of Σ_A. For |Λ_p| > |Γ_wind| the mass flux should be blocked by the planet, and the accretion rate onto the planet should be much smaller than the mass flux through the unperturbed disc. Unlike in the viscous case, here we do not expect an increase in the mass flow into the gap as mass builds up near the gap edge. This is because we assume that Σ_A is constant, and the torque acting in the active layer is constant, so any gas that accumulates near the gap edge will just join the dead zone, and hence will not flow towards the planet.

In the limit (|Λ_p| ≪ |Γ_wind|), mass should flow relatively unimpeded towards the planet, which should be able to accrete at close to the mass flow rate through the unperturbed disc. The transition between these behaviours should occur when |Λ_p| ≲ |Γ_wind|, and the accretion rate should be smaller than the mass flow rate through the unperturbed disc, but the accretion flow should not be completely blocked.

An expression for Λ_p obtained from the impulse approximation (Lin & Papaloizou 1986b) is often used in studies of disc-planet interactions, and takes the following form when the planet Hill radius exceeds the pressure scale height $${\text{Λ}}_{\text{p}}=\{\begin{array}{ll}\text{sign}\left(R-R_{\text{p}}\right)\frac{8}{9\pi }\frac{q_{\text{p}}^{2}GM}{2R_{\text{p}}}{\left(\frac{R_{\text{p}}}{{\text{Δ}}_{\text{p}}}\right)}^4,\hfill & \text{if}\left|R-R_{\text{p}}\right|>R_{\text{H}}\hfill \\ \frac{\left(R-R_{\text{p}}\right)}{R_{\text{H}}}\frac{8}{9\pi }\frac{q_{\text{p}}^{2}GM}{2R_{\text{p}}}{\left(\frac{R_{\text{p}}}{{\text{Δ}}_{\text{p}}}\right)}^4,\hfill & \text{otherwise}\hfill \end{array}$$(19)

where Δ_p = Max (R_H, R – R_p), q_p = M_p/M_*, and R_H is the Hill radius. The second term in Eq. (19) is a modification that ensures Λ_p is continuous as it passes through the planet’s location. Using 3D hydrodynamical simulations of disc-planet interactions similar to those presented here, D’Angelo & Lubow (2010) have examined how Λ_p varies in different disc models and for different planet masses, and they provide the following expression based on fits to their simulation results (so the expression agrees much more closely with hydrodynamical simulations than does Eq. (19)): $${\text{Λ}}_{\text{p}}=\frac{1}{2}ℱ\left(x,p,q\right){\text{Ω}}_{\text{p}}^2{R}_{\text{p}}^2{q}_{\text{p}}^2{\left(\frac{R_{\text{p}}}{R_{\text{H}}}\right)}^4,$$(20)

where $$ℱ\left(x,p,q\right)=\left(p_1{e}^{\left(-\frac{{\left(x+p_2\right)}^2}{p_3^2}\right)}+p_4{e}^{\left(-\frac{{\left(x-p_5\right)}^2}{p_6^2}\right)}\right)\times \tan h\left(p_7-p_{8}x\right).$$(21)

We note that the values p and q in $ℱ\left(x,p,q\right)$ refer to the power-law indices that define the density and temperature profiles in Eqs. (6) and (7), and the factor of 1/2 in Eq. (20) should only be included to account for non-linear effects when gap opening planets are considered. Similarly in this limit we have x = (R – R_p)/R_H. The constants p₁, p₂, etc. in Eq. (21) depend on the disc structure, and for the power-law indices p = −1.5 and q = −1 we consider in this work the values of these constants are given in Table 1 of D’Angelo & Lubow (2010).

Figure 4 shows the values of Λ_p from Eqs. (19) and (20), and also the values of Γ_wind for each of the layered disc models. It is clear that Eq. (19) predicts torque values that are too large. Comparing the values from Eq. (20) with the wind-driven torques for the different disc models, we expect the accretion rate onto the planet to be close to the mass flow rate through the unperturbed disc for Σ_A = 0.1 g cm⁻². For Σ_A = 1 g cm⁻², we expect the accretion rate to be moderately impeded, and for Σ_A = 10 g cm⁻² we expect the accretion flow towards the planet to be blocked, and hence for there to be a very significant reduction in the accretion rate onto the planet compared to the unperturbed mass flux through the disc.

Fig. 4 Planet torque per unit mass calculated according to Eqs. (19) (dashed black line) and (20) (solid black line), and wind-driven torques per unit mass for different values of ΣA such that one hemisphere of the disc has a mass flux of 5 × 10−9 M⊙ yr−1 travelling through its active layer.

4 Results

4.1 Disc structure

Figure 5 shows the evolution of the midplane density in each of the models, and a common feature is the formation of a deep annular gap around the planet’s orbital location. As described in our earlier work (Lega et al. 2021, Lega et al. 2022), the low viscosity models all show the development of a strong vortex that grows via the Rossby Wave Instability at the density maximum associated with the edge of the gap (de Val-Borro et al. 2006; Lovelace et al. 1999). The vortex is a source of strong spiral waves that propagate outwards and dissipate, forming a secondary gap and a secondary density maximum that can be observed as a high density ring just exterior to 10 au in Fig. 5. The vortices migrate inwards and dissipate as they move away from the pressure bump, as can be seen in the rightmost panels of the figure, and the low viscosity employed in these models means that the ring and gap features generated by the vortex survive over long timescales. These models maintain long memories of earlier evolution. As expected, the Viscous model does not show the development of a long lived vortex or the secondary features that might arise from one, and instead the disc exterior to the planet has a smooth structure on which spiral density waves are superposed.

The azimuthally averaged surface density profiles for each of the models are shown in Fig. 6. Here we see significant variation between the models. The gap for the Viscous model remains relatively narrow, since the formation of the gap causes the viscous flow of material into the gap from both the inner and outer disc. The No Wind model shows the formation of a much deeper and wider gap than the Viscous model due to the very low viscosity adopted in this run, allowing the planet to more strongly repel material on both sides of its orbit. The inner-outer disc asymmetry observed in this run is likely due to the presence of a vortex at the outer edge of the gap during most of the simulation. The spiral waves emitted by this, and its eventual dissipation, would have provided a source of mass flow towards the planet at the outer gap edge.

The models with wind-driven accretion flows all show the same strongly asymmetrical gap structure, where the asymmetry this time is significantly larger than in the No Wind case. This arises because the wind-driven torque drives the disc gas towards the planet in the outer disc, opposite to the direction that the planet torque tries to drive material, and away from the planet and towards the star in the inner disc, where the planet torque acts in the same direction as the wind-driven torque. We note that keeping the planet on a fixed orbit rather than allowing it to migrate somewhat exaggerates the asymmetry of the gap, as shown by the results in Lega et al. (2022). Allowing the planet to accrete gas efficiently also contributes to this asymmetry, by blocking the flow from the outer to the inner disc even in models where the wind-driven flow is not strongly impeded by the tidal torque from the planet. The Σ_A = 0.1 g cm⁻² run has the shallowest gap, showing that that gas is able to flow into the gap region and is not strongly impeded by the planet. The Σ_A = 10 case, however, displays a very deep gap indicating that gas flow into the gap is strongly diminished by the planet. The Σ_A = 1 case is intermediate between these two, showing it lies in the transition between weakly and strongly impeded radial gas flow. Finally, it is worth noting that an interesting implication of these different disc structures, all arising from the same global accretion rate through the disc, is that accurately inferring planet masses from observations of gap structures in protoplanetary discs becomes much more difficult if accretion through the disc is driven by an unknown combination of stresses originating from turbulent viscosity and a magnetised wind.

4.2 Gas accretion rates

The accretion rates onto the planets in each model are shown in the left panel of Fig. 7, and the results are in line with the expectations described in Sect. 3. The Viscous model shows a steady accretion rate corresponding to Ṁ_p ≃ 12 M_Jup Myr⁻¹, such that a Jovian mass planet will double its mass in ≲10⁵ yr. Time-averaged radial profiles of the vertically and azimuthally integrated mass fluxes through the discs with embedded planets are shown in the right panel of Fig. 7. Exterior to the planet, the viscous disc has a mass flux towards the planet close to the unperturbed value. The disc interior to the planet shows an outward mass flow, again towards the planet. Hence, the planet at r_p = 5.2 au is being fed from both sides, and the jump in the mass flux rate at the planet location corresponds to the steady accretion rate recorded in the left panel.

The planet in the Σ_A = 0.1 model has a steady accretion rate of Ṁ_p ≃ 6 M_Jup Myr⁻¹, giving a mass doubling time of 2 × 10⁵ yr. The mass flux far from the planet in the right panel is again close to the unperturbed value, and decreases somewhat closer to the planet. The jump in the radial mass flux at the planet location again corresponds to the accretion rate recorded in the left panel, showing that the gas supplied to the planet by the wind-driven torque is largely accreted by it, although a fraction of this gas flows past the planet and into the inner disc.

The planet in the Σ_A = 1 model sustains an accretion rate of Ṁ_p ≃ 2 m_Jup Myr⁻¹, giving a mass doubling time of ≃5 × 10⁵ yr. The right panel of Fig. 7 shows a similar picture to the Σ_A = 0.1 model, with the mass flux at large radius being equal to the unperturbed value and decreasing somewhat at radii closer to the planet. The jump in the mass flux at the planet location again corresponds to the accretion rate displayed in the left panel.

The Σ_A = 10 runs shows the greatest influence of the planet torque on the gas accretion rate. Here Ṁ_p ≃ 10⁻¹ M_Jup Myr⁻¹, giving a mass doubling time of 10 Myr. In agreement with the discussion in Sect. 3, the wind-driven radial gas flow towards the planet is almost completely blocked by the planet torque, giving rise to the dramatically reduced accretion rate.

Finally, the No Wind model shows an averaged accretion rate of Ṁ_p ~ 10⁻² M_Jup Myr⁻¹. This small accretion rate indicates that the inclusion of a small viscosity in the wind-driven models does not significantly contribute to the gas flow towards the planets. The right panel of Fig. 7 shows the spiral waves launched by the planet into the outer disc cause a modest outward flow of mass in the No Wind model, and beyond 20 au the disc shows no net flow of mass as expected. Closer to the planet in the gap region, there is a very weak inflow towards the planet which arises because of the small imposed viscosity acting within the very deep gap.

Figure 8 shows time averages of the azimuthally integrated radial mass fluxes in the three simulations with wind-driven accretion flows. The time averages were computed using 100 snapshots that were output with an interval of one planet orbit, equivalent to 11.86 yr, starting at times slightly before the times that correspond to the surface density profiles displayed in Fig. 6. The panels show the accretion flow being largely confined to a narrow region in the vertical direction close to the transition between the active and dead zones, a feature which persists all the way into the gap. It is noteworthy that in all models the full column density, Σ(0), falls below Σ_A in the gap, such that the integrated accretion flow must fall below its unperturbed value there because the wind-driven torque per unit mass applied to the active column is constant in the wind model we have adopted. Figure 8 illustrates the fact that in 3D hydro-dynamical simulations the flow is considerably more complex than suggested by the simple 1D picture presented in Sect. 3, and the influence of resonances such as the 2:1 outer Lindblad resonance induce spatial and temporal variability in the flow. Nonetheless, the simple picture of the long term net accretion flow onto the planet being determined by the balance between the applied wind-driven torque and the planet tidal torque seems to hold, and explains the behaviour observed in the simulations as a function of Σ_A and the applied viscosity.

Fig. 5 Contours showing the evolution of the midplane density for each of the models.

Fig. 6 Azimuthally averaged surface density profiles from all disc models computed once each run had evolved for ~49 000 yr after the value of K was increased from 6 to 30 in Eq. (4). The planet positions are indicated by the filled circles. Note that the inner boundary condition influences the surface density profiles close to the inner boundary. The profiles shown for each run correspond to the right-hand-most midplane density contours shown in Fig. 5.

5 Discussion

5.1 Combined migration and gas accretion rates

In Lega et al. (2022) we presented simulations of migrating and non-accreting Jovian mass planets embedded in layered disc models in which the total mass flux towards the star is ṁ = 10⁻⁸ M_⊙ yr⁻¹. Here we have presented simulations of accreting and non-migrating Jovian mass planets in similar layered models. In each of these scenarios, the behaviour is determined by the ability or otherwise of the planetary torque to block the wind-induced mass flow occurring in the active layer. In both studies we observe a dramatic change in behaviour when considering the two extreme cases simulated, Σ_A = 0.1 g cm⁻² and Σ_A = 10 g cm⁻², with the transition in behaviour occurring between these values.

For migrating planets, Lega et al. (2022) show that once the vortex that forms in the layered disc models has dissipated, the long term migration behaviour depends on the mass flux through the gap. When the gas flows unimpeded through the gap then we can consider the disc to consist of two distinct parts: an inert dead zone that contains most of the mass, and a low mass active zone that sustains a constant mass flux at all radii. The planet forms a gap in the dead zone and sits in a location that minimises the net torque arising from its interaction with the inner and outer disc. The flow of gas through the gap maintains only a small amount of mass in the gap that barely influences the migration rate. Hence in this case we have very slow inward migration of the planet. The results in this paper show that this slow mode of migration will be accompanied by rapid gas accretion, since an unimpeded accretion flow reaches the planet. We note, however, that this conclusion depends on the efficiency with which the planet can accrete the gas that is supplied to it. As we discuss below, there is considerably uncertainty about what this efficiency is.

If we consider what would happen when a planet is undergoing slow migration and is able to accrete gas at close to the supply rate, then the flow from the outer to the inner disc would be interrupted. Over time the wind-driven torque acting on the inner disc would deplete it, and this would have the effect of modifying the balance of Lindblad torques acting on the planet through interaction with the inner and outer disc. This is one way in which allowing a planet to accrete gas while it migrates could in principle modify the migration behaviour compared to a scenario in which the planet does not accrete gas. We speculate that in this scenario the planet would migrate away from the outer gap edge until the torque it experiences from it becomes very small, at which point its migration would stall. Gas accretion, however, would continue to occur, and hence there seems to be little prospect of halting accretion and migration simultaneously. Simulations that explore this scenario will be presented in a forthcoming publication.

When the gas flow is blocked by the planet torque, then the migration picture changes. As the planet migrates inwards and away from the outer edge of the gap, the laminar accretion flow fills in the gap behind the planet. Hence, migration is sustained by the gas inflow and the migration rate is determined by the rate at which the gap is refilled, leading to the following estimate for the migration rate $${\dot{r}}_{\text{p}}=\frac{\dot{M}}{2\pi r_{\text{p}}\text{Σ}},$$(22)

where Σ is the total surface density at the edge of the gap. For a planet at r_p = 5.2 au, the simulations in Lega et al. (2022) show the migration speed is ~15 au Myr⁻¹ when the accretion flow through the disc is 10⁻⁸ M_⊙ yr⁻¹ and Σ_A = 10 g cm⁻². The results in this paper show that the gas accretion rate will be very small because the rate at which gas is supplied to the planet is ṁ ~ 0.1 M_Jup yr⁻¹ while the planet is in this faster mode of migration.

To summarise, a Jovian mass planet that can efficiently accrete any gas that is supplied to it, embedded at r_p = 5.2 au in a disc with a wind-induced accretion flow of 10⁻⁸ M_⊙ yr⁻¹, will undergo slow migration at a speed ${\dot{r}}_{\text{p}}~3$ au Myr⁻¹, and rapid gas accretion at arate of Ṁ_p ~ 5 M_Jup Myr⁻¹ if Σ_A = 0.1 g cm⁻². If Σ_A = 10 g cm⁻², however, then the migration speed will be faster ${\dot{r}}_{\text{p}}~15$ au Myr⁻¹, and the accretion rate will reduce to Ṁ_p ≲ 0.1 M_Jup Myr⁻¹. Migration and gas accretion rates that are intermediate between these values will be obtained for values of Σ_A that lie between these limiting cases, as demonstrated here and in Lega et al. (2022) for the case of Σ_A = 1 g cm⁻².

Fig. 7 Gas accretion and radial mass flow rates. Left panel: gas accretion rates onto the planets versus time for all disc models. Moving from top to bottom, the horizontal black dotted lines correspond to accretion rates of 10, 1, 0.1, and 0.01 MJup Myr−1. The grey shaded region highlights the temporal domain within which the value of K in Eq. (4) is increased from 6 to 30. For the Viscous run and those with different values of ΣA, the increase coincides with the discontinuous changes to the accretion rates that can be observed. The value of K for the No Wind case was increased at time ~87 000 yr. Right panel: radial mass flux profiles for all disc models containing Jovian mass planets. The upper and lower horizontal black dotted lines correspond to radial mass fluxes of 0 and 10−8 M⊙ yr−1, respectively. The mass fluxes were computed at the times corresponding to the surface density profiles displayed in Fig. 6.

Fig. 8 Contour plots showing time-averaged and azimuthally-integrated radial mass flux distributions in the various layered disc models with embedded planets. The cyan dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. The vertical white dotted line denotes the location of the 2:1 outer Lindblad resonance. Note that only the upper hemisphere of the disc is simulated.

5.2 Behaviour as a function of orbital radius

We have considered gas accretion onto Jovian mass planets orbiting at r_p = 5.2 au in layered disc models. Here we consider how the evolution depends on orbital radius, assuming that Σ_A varies weakly with stellocentric distance. As discussed in Lega et al. (2022) and Sect. 3, the transition between rapid and slow migration (and hence between slow and rapid gas accretion) occurs when Λ_p ~ Γ_wind. Using Eqs. (20) and (18) we can determine the critical value of Σ_A that corresponds to this transition in behaviour. We note that the function F defined in Eq. (21) has a maximum value ${ℱ}_{\max }~0.02$. The critical value, ${\widetilde{\text{Σ}}}_{\text{A}}$, can be written $${\widetilde{\text{Σ}}}_{\text{A}}=0.285\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-2}\left(\frac{\dot{m}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right){\left(\frac{M_{\text{p}}}{M_{\text{jup}}}\right)}^{-2/3}{\left(\frac{r_{\text{p}}}{5.2\text{\hspace{0.17em}au}}\right)}^{-1/2},$$(23)

where ṁ corresponds to the total mass flux through the disc due to the wind-driven torque. When ${\text{Σ}}_{\text{A}}>{\widetilde{\text{Σ}}}_{\text{A}}$ the accretion flow will be blocked and we expect fast migration, and when ${\text{Σ}}_{\text{A}}<{\widetilde{\text{Σ}}}_{\text{A}}$ the accretion flow will not be blocked and migration will be slow.

As remarked upon in Lega et al. (2022), the $r_{\text{p}}^{-1/2}$ dependence of ${\widetilde{\text{Σ}}}_{\text{A}}$ means that a migrating planet can transition between the fast and slow modes of migration as it migrates inwards. Figure 9 shows how ${\widetilde{\text{Σ}}}_{\text{A}}$ varies with radius in a disc with ṁ = 10⁻⁸ M_⊙ yr⁻¹ containing a Jovian mass planet, and shows the transition radius (where the lines cross in the figure) between fast and slow migration for different values of Σ_A. An important issue is determining an appropriate value for Σ_A. Although X-rays and cosmic rays are important sources of ionisation in protoplanetary discs, and penetrate into columns of depth ~10 g cm⁻² and 100 g cm⁻², respectively, the strong magnetic coupling in the surface layers of protoplanetary discs that allows a magnetised wind to be launched is thought to depend on the ionisation of sulphur and carbon atoms by UV photons from the star (Perez-Becker & Chiang 2011). The column density of this layer is such that we expect Σ_A < 1 g cm⁻². As shown by Fig. 9, values of 0.1 ≤ Σ_A ≤ 1 g cm⁻² give rise to locations where the planet transitions between fast and slow migration (and vice versa for gas accretion) that cover a very wide range of radii in protoplanetary discs, from r ~ 50 au when Σ_A = 0.1 to r ~ 0.4 au when Σ_A = 1.

If we suppose the disc has Σ_A = 0.5 g cm⁻² at all disc radii, and a Jovian mass planet starts migrating at r_p = 5.2 au, then ${\text{Σ}}_{\text{A}}>{\widetilde{\text{Σ}}}_{\text{A}}\left(r=5.2\text{\hspace{0.17em}au}\right)$ and migration should be in the fast regime with gas accretion being in the slow regime. As the planet migrates inwards ${\widetilde{\text{Σ}}}_{\text{A}}$ increases, and at r_p ~ 1.6 au we have ${\text{Σ}}_{\text{A}}={\widetilde{\text{Σ}}}_{\text{A}}$, so that migration should start to enter the slow regime and gas accretion should transition to the rapid regime. This evolutionary scenario nicely explains the large number of giant planets observed with orbital periods P ≳ 100 days, but also suggests that giant planet masses should on average show an inverse relationship with orbital period near the transition radius.

Figure 10 displays mass-period diagrams for giant planets discovered by radial velocity surveys (top panel is for all planets, bottom panel is for low eccentricity planets with e < 0.3 that are assumed to have not undergone strong gravitational scattering after removal of the gas disc)², and they do not show strong evidence for a significant increase in planet masses at orbital periods less than a few hundred days, independent of the eccentricity range being considered. A local minimum can be observed at ~400 days, but this does not seem to be statistically significant, and is not consistent with the idea that Jovian mass planets that form at ~5.2 au, migrate inwards rapidly, while slowly increasing their masses, and then slow down their migration and significantly increase their gas accretion rates as they reach orbital radii ~1 au.

Fig. 9 ${\widetilde{\text{Σ}}}_{\text{A}}$ versus radius for a Jovian mass planet located in a disc with ṁ = 10−8 M⊙ yr−1. Also shown are various values of ΣA, indicating where an embedded planet would transition between fast migration and slow accretion to slow migration and fast accretion.

Fig. 10 Mp sin i versus period for exoplanets discovered by radial velocity surveys. Top panel is for all planets, bottom panel is for low eccentricity planets where e < 0.3. Red points show the average masses of planets binned in radius, where the bin width is 0.4 in log10 P. Error bars are for one standard deviation of the planet mass distribution in each bin.

5.3 Gas supply rate versus accretion rate

The discussion in the previous subsection suggests that giant planets somehow avoid accreting gas at the rate supplied by the disc over significant timescales, but instead either accrete for only relatively short time periods because they form towards the end of the disc lifetime, or accrete the supplied gas at a significantly lower rate than assumed in simulations such as the ones presented here.

There have been a number of recent high resolution 3D studies of the gas flows within the Hill spheres and surrounding regions for planets of various masses embedded at different radii in protoplanetary discs (e.g. Fung et al. 2019; Lambrechts et al. 2019; Szulágyi et al. 2016; Moldenhauer et al. 2021, 2022). These studies have often utilised radiation-hydrodynamics codes, and hence are able to examine the cooling of the gas and the rate at which it accretes onto the planet. Moldenhauer et al. (2021, 2022) examined gas accretion onto cores ranging in mass between 1 and 10 M_⊕, orbiting close to the star at 0.1 au, and showed that the cooling and contraction of the planet envelope can be switched off completely by the continuous advection of high entropy gas into the Hill sphere from the surrounding disc. The advection occurs at a rate that scales with the local orbital period, and so while this is likely to be an important effect for planets close to the star, it is not clear that it will remain so at large stellocentric distances. This point is supported by the 3D radiation-hydrodynamic simulations presented by Lambrechts et al. (2019), who showed that the advection of gas into the Hill sphere is still an important phenomenon for planets orbiting at 5.2 au, but that Jovian mass planets orbiting there accrete gas at a rate that leads to the mass doubling time being ≲ 10⁵ yr. Hence, the conclusion of that study was that thermal effects occurring within the Hill sphere are unlikely to slow down the accretion of gas if it is supplied by the disc at a rate of ṁ ~ 10⁻⁸ M_⊙ yr⁻¹. Other effects might be important in slowing down the accretion of gas once it has been supplied to the vicinity of the planet, for example magnetic effects occurring in the Hill sphere, such as considered by Gressel et al. (2013), or an increase in the magnetic torque in the gap leading to a faster radial flow of gas past the planet, such that it accretes at slower rate. At the present time, however, we do not have evidence of the efficacy of these effects in slowing down accretion, and it is difficult to escape the conclusion that the final masses of most giant planets are determined by the lengths of time that they are present within their protoplanetary discs, rather than by a physical process that acts to slow the accretion of gas that is supplied by the protoplanetary disc, or a process that slows the supply to the planet’s location deep inside the Hill sphere.

6 Conclusions

We have presented the results of 3D hydrodynamical simulations of accreting giant planets embedded in protoplanetary discs that globally sustain accretion flows towards the star of ṁ = 10⁻⁸ M_⊙ yr⁻¹. We consider a classical α viscous disc model, and a sequence of very low viscosity models in which a laminar accretion flow is driven in the surface layers of the disc by an external torque, that is supposed to mimic the effects of a magnetocentrifugally driven wind. In these latter models, the column density of the actively accreting layer is varied in the range Σ_A = 0.1−10 g cm⁻², such that the radial speed, and the torque per unit mass acting on the accreting gas, increases as the column density of the actively accreting gas decreases, in order to maintain the given accretion rate.

The main result obtained is that the accretion rate onto the planet varies significantly between the models, in spite of the fact that the accretion rate supplied through the disc is always ṁ = 10⁻⁸ M_⊙ yr⁻¹. The viscous model produces an accretion rate that leads to a mass doubling time for a Jovian mass planet of ~10⁵ yr (0.1 Myr), in agreement with numerous previous studies. The wind-driven models, however, produce accretion rates such that the mass doubling times vary between 0.2 and 10 Myr for Σ_A = 0.1 and 10 g cm⁻², respectively.

Fundamentally, this result arises because viscous discs adjust their accretion flows in the presence of a planet-induced gap, such that a large accretion rate through the gap and onto the planet can be maintained. The model for the wind-induced torque we have adopted, on the other hand, assumes the torque per unit mass acting on the accreting gas does not change in the gap regions, and hence for larger values of Σ_A the tidal torque from the planet can overwhelm the torque acting in the accreting layer, and effectively block the accretion flow. For small values of Σ_A, the planet torque is smaller in magnitude than that acting in the disc surface layers, and hence the accretion flow is relatively unimpeded by the planet, such that a large accretion rate onto it can be maintained. For the parameters adopted in our models, the transition in behaviour between fast and slow accretion occurs between 0.1 < Σ_A < 10 g cm⁻².

In a recent study (Lega et al. 2022), we considered the migration of non-accreting Jovian mass planets in layered disc models, similar to those presented here. In that case, there is also a transition in behaviour such that migration is fast when Σ_A = 10 g cm⁻² and slow when Σ_A = 0.1 g cm⁻². Hence, we expect that planets will accrete slowly when local disc conditions allow for rapid inward migration, and will accrete rapidly when conditions cause migration to be slow (assuming no process operates that can impede accretion onto the planet of gas that is supplied by the accretion flow in the background disc). We will present the results of an ongoing study of migrating and accreting giant planets embedded in layered disc models in a forthcoming publication.

The model of the wind-driven accretion torque we have presented is highly simplified, and serves the purpose of providing a framework for understanding how giant planets evolve when embedded in non-viscous and accreting protoplanetary discs. In future work, we will improve on the simplifications of the present model by employing non-ideal MHD simulations, and examine how giant planets migrate and accrete in discs that include a more complete treatment of the relevant physics.

Acknowledgements

We thank the anonymous referee for a constructive report that helped improve the paper. R.P.N. acknowledges support from STFC through grants ST/P000592/1 and ST/T000341/1. This research utilised Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT (http://doi.org/18.5281/zenodo.438845). This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. E.L. and A.M. acknowledge support by DFG-ANR supported GEPARD project (ANR-18-CE92-0044 DFG: KL 650/31-1). We also acknowledge HPC resources from GENCI DARI n.A0120407233 and from “Mesocentre SIGAMM” hosted by Observatoire de la Côte d’Azur. A.M. acknowledges support from the ERC project HolyEarth - 101019380. E.L. wishes to thank Alain Miniussi for maintenance and re-factorisation of the code FARGOCA.

Appendix A Comparison between NIRVANA and FARGO3D

Here we present a comparison between the NIRVANA and FARGO3D codes when applied to a standard disc-planet interaction scenario. A run was performed using the disc model and setup described in Sect. 2. We adopted a model without an accreting layer and with viscous parameter α = 10⁻⁵, as in the No Wind case that we presented in the main paper. A 20 M_⊕ planet was inserted in the disc at 5.2 au on a fixed circular orbit and the radial mass flux in the disc was computed in the same way as was done for the mass fluxes presented in the right panel of Fig. 7. Figure A.1 shows the results obtained by the two codes, where the lines represent the azimuthally and meridionally integrated mass fluxes at each radius, averaged between the times corresponding to 400-480 planet orbits. As expected, the outward and inward propagating spiral waves lead to an angular momentum flux, and an associated mass flux, that is concentrated around the planet where the waves damp. Further from the planet, the mass flux essentially disappears because a wind torque is not applied in this case and the viscosity is very small. Good agreement is obtained between the two codes.

Fig. A.1 Time-averaged mass flux versus radius for runs performed using NIRVANA and FARGO3D.

References

All Figures

	Fig. 1 Contour plots showing the radial velocity distributions in the various layered disc models prior to insertion of the planets. The black dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. We note that only the upper hemisphere of the disc is simulated.
In the text

	Fig. 2 Contour plots showing time-averaged and azimuthally-integrated radial mass flux distributions in the various layered disc models. The cyan dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. Note that only the upper hemisphere of the disc is simulated.
In the text

	Fig. 3 Radial profiles of the gas accretion rates in the viscous and externally torqued disc models.
In the text

	Fig. 4 Planet torque per unit mass calculated according to Eqs. (19) (dashed black line) and (20) (solid black line), and wind-driven torques per unit mass for different values of ΣA such that one hemisphere of the disc has a mass flux of 5 × 10−9 M⊙ yr−1 travelling through its active layer.
In the text

	Fig. 5 Contours showing the evolution of the midplane density for each of the models.
In the text

Fig. 6 Azimuthally averaged surface density profiles from all disc models computed once each run had evolved for ~49 000 yr after the value of K was increased from 6 to 30 in Eq. (4). The planet positions are indicated by the filled circles. Note that the inner boundary condition influences the surface density profiles close to the inner boundary. The profiles shown for each run correspond to the right-hand-most midplane density contours shown in Fig. 5.

In the text

Fig. 7 Gas accretion and radial mass flow rates. Left panel: gas accretion rates onto the planets versus time for all disc models. Moving from top to bottom, the horizontal black dotted lines correspond to accretion rates of 10, 1, 0.1, and 0.01 MJup Myr−1. The grey shaded region highlights the temporal domain within which the value of K in Eq. (4) is increased from 6 to 30. For the Viscous run and those with different values of ΣA, the increase coincides with the discontinuous changes to the accretion rates that can be observed. The value of K for the No Wind case was increased at time ~87 000 yr. Right panel: radial mass flux profiles for all disc models containing Jovian mass planets. The upper and lower horizontal black dotted lines correspond to radial mass fluxes of 0 and 10−8 M⊙ yr−1, respectively. The mass fluxes were computed at the times corresponding to the surface density profiles displayed in Fig. 6.

In the text

Fig. 8 Contour plots showing time-averaged and azimuthally-integrated radial mass flux distributions in the various layered disc models with embedded planets. The cyan dotted line indicates the location where the column density measured from the disc surface corresponds to ΣA. The vertical white dotted line denotes the location of the 2:1 outer Lindblad resonance. Note that only the upper hemisphere of the disc is simulated.

In the text

	Fig. 9 ${\widetilde{\text{Σ}}}_{\text{A}}$ versus radius for a Jovian mass planet located in a disc with ṁ = 10−8 M⊙ yr−1. Also shown are various values of ΣA, indicating where an embedded planet would transition between fast migration and slow accretion to slow migration and fast accretion.
In the text

	Fig. 10 Mp sin i versus period for exoplanets discovered by radial velocity surveys. Top panel is for all planets, bottom panel is for low eccentricity planets where e < 0.3. Red points show the average masses of planets binned in radius, where the bin width is 0.4 in log10 P. Error bars are for one standard deviation of the planet mass distribution in each bin.
In the text

	Fig. A.1 Time-averaged mass flux versus radius for runs performed using NIRVANA and FARGO3D.
In the text