© ESO 2019

1 Introduction

Gas giants form an important class of planets in the solar system as well as in extrasolar planetary systems. As the most massive planetary class, their gravity determines the evolution of planetary systems after their formation (Davies et al. 2014). Ever since the first discovery of gas giant exoplanets (Mayor & Queloz 1995), their occurrence rates and correlations with host star properties have been mapped observationally. It is clear now that gas giants orbit primarily around stars that are rich in heavy elements (Gonzalez 1997; Fischer & Valenti 2005; Mayor et al. 2011); this is particularly pronounced for gas giants on eccentric orbits that likely resulted from planet–planet scattering in systems of multiple gas giants (Dawson & Murray-Clay 2013; Buchhave et al. 2018). The significant increase in the formation rate of gas giants with metallicity supports the core accretion model for the formation of such planets (Pollack et al. 1996).

The formation of a planetary core and subsequent gas accretion must happen during the gaseous disc phase that lasts between 1 and 10 Myr. The following formation of one or more gas giants can have a dramatic influence on further planet formation. The growing gas giants will form gaps, first in the dust and, if the planet becomes sufficiently massive, also in the gas. These gaps in dust (Paardekooper & Mellema 2004; Lambrechts et al. 2014; Bitsch et al. 2018; Ataiee et al. 2018) and the gas distributions (Ward 1997; Lubow & D’Angelo 2006; Crida et al. 2006) will cut off or at least heavily modify the influx of material onto the inner planets. The edges of planetary gas gaps can destabilise, become turbulent, and thus influence further disc evolution by changing the strength of turbulence or regulating their own depth (Hallam & Paardekooper 2017). The gap structure can then have a feedback effect on the planet. Gaps slow planetary migration through the parent disc (Lin & Papaloizou 1986) and both gap depth and migration speed can affect accretion independently (Dürmann & Kley 2017; Kanagawa et al. 2018).

Migrating gas giants will also excite the orbits of other planets in the system (Raymond et al. 2016) and change their migration speed (Masset & Snellgrove 2001). They can further influence the evolution of systems through trailing swarms of pebbles that accumulate at their gap edge and form seed locations for the growth of new planets (Ronnet et al. 2018). Therefore, an important priority of planet formation theory is to clarify the gas accretion rates of young gas giants.

The early 1D gas accretion models of Mizuno (1980) determined that beyond a certain critical core mass, the gaseous envelope of a young planet cannot remain in hydrostatic equilibrium. Subsequent work has then often focused on finding this critical mass as a function of disc parameters with increasing model complexities (Bodenheimer & Pollack 1986; Wuchterl 1990; Pollack et al. 1996). Important ingredients that allow planets to pass the critical core mass were firstly the build-up of mass via accretion of solids, and secondly the cessation of their accretion, which allows the circumplanetary gas envelope to cool, contract, and go into runaway gas accretion. Those approaches have been refined over the years, for example with more detailed opacity laws and structure models (Piso & Youdin 2014; Piso et al. 2015).

However, effects of higher dimensionality have been shown to be important even before crossing critical mass. Horseshoe flows that penetrate the planetary Hill-sphere in the midplane and interfere with the contracting planetary envelope can only be accurately modelled in at least two dimensions (Ormel et al. 2015). Vertical outflows as well as cooling of the gas into space even necessitate full 3D models (Tanigawa et al. 2012; Lambrechts & Lega 2017).

As well as dimensionality, the treatment of the energy equation(s) for gas dynamics also plays a crucial role in calculating the release of initial accretion heat, allowing the accretion process to continue. The treatment of the energy equation has gone through several iterations, ranging from isothermal (e.g. Stone & Norman 1992; D’Angelo et al. 2003; Tanigawa et al. 2012; Fung et al. 2015) or adiabatic, or mixes between them (e.g. Machida et al. 2010), through Newtonian-cooling with density-dependent cooling times (Gressel et al. 2013; Zhu et al. 2016), to radiative transport simulations (e.g. Ayliffe & Bate 2009, 2012; D’Angelo & Bodenheimer 2013) which approximate the transitions from optically thick to thin regions.

Since 2D and 3D simulations are computationally expensive, there are limitations on the maximum size of the simulation domain around the planet and limitations on how close to the planet the simulation can correctly resolve the physics. This has led to the two following approaches.

Local shearing boxes accurately resolve the planetary Hill- sphere by the use of mesh-refinement techniques (e.g. Zhu et al. 2016), but can incur problems when representing gap-opening planets or maintaining the global energy balance. Global simulations have large radial and full azimuthal extent, and are sophisticated enough to resolve the planetary Hill-sphere (e.g. Fung et al. 2015; Lambrechts & Lega 2017) or even the planetary radius (e.g. D’Angelo et al. 2003; D’Angelo & Bodenheimer 2013; Szulágyi et al. 2014), but usually suffer from run-time restrictions that render the simulation of evolutionary timescales of planets impractical. A notable exception here is the work by Ayliffe & Bate (2012) who managed to simulate the runaway collapse of planetary envelopes in 3D radiative SPH-hydrodynamics globally with self-gravity.

Finally, accretion of the planet is also subject to different treatments. It is general practice to either use a global simulation with a certain fraction of mass taken out of the central 10–50% of the planetary Hill-radius (Kley 1989; Klahr & Kley 2006), or simulations with very high spatial resolution to let the cooling process proceed naturally, without taking out mass artificially (D’Angelo & Bodenheimer 2013).

In this work we aim to measure mass accretion rates onto giant planets while relaxing physical assumptions as much as possible. To this end, we present a grid-based, global, 3D radiation hydrodynamical two-temperature setting where neither the accretion rate nor the accretion luminosity are imposed, but are instead self-consistently calculated from the contraction of the gas, viscous heating, and heat advection. We resolve the planetary Hill-sphere down to 5% of the Hill-radius, where this again is maximally resolved by up to ≈8–20 grid cells in order to assess the convergence of the results with increasing resolution of the Hill-radius. We also investigate the influence of two numerical parameters, namely the resolution and gravitational smoothing length, as well as of the key physical parameter that regulates the cooling of the planetary envelope, that is the opacity law, on the measured accretion rates.

This is intended to be part feasibility study, part physical investigation. As the outcome of any accretion signal should be an unambiguous increase in mass, we searched the literature for helpful hints concerning the required parameters to maximise this increase. The studies of Ayliffe & Bate (2009) and Bodenheimer et al. (2013), where our setup mostly resembles the former, showed that a Saturn-mass planet should yield a maximum in accretion rate. This is caused on one hand by the fact that accretion rates increase with planetary masses, but also by the fact that massive, gap-opening planets deplete their own feeding zones through the action of gravitational torques, decreasing accretion rates as they grow even further. Therefore, we focus here on Saturn-mass protoplanets and leave the dependence of the results on the mass of the protoplanet to a future study.

As our planets are massive enough to form a gap in the gas distribution and are beyond their pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018), we do not set any additional heating sources from solids, which has been explored in a similar setting to ours by Lambrechts & Lega (2017) for low-mass planets. A similar study to theirs was also done by D’Angelo & Bodenheimer (2013) for low-mass planets between 5 and 15 Earth masses. These latter authors, however, did a more careful modelisation of the thermodynamics, particularly the dependence of the adiabatic index γ on temperature, and also assumed matter-photon energy equilibrium, which implies a one-temperature approach, an approximation which we can relax in order to perform a more realistic exploration of the cooling behaviour of high-mass planets.

This paper is organised as follows. In Sect. 2 we present a brief overview of the physical model used and the free parameters it depends on. Section 3 outlines the numerical methods used to solve the physical model and explains three individual steps we use to maintain a reasonable computation time. Section 4 discusses a limited set of measured accretion rates and their dependency on numerical resolution, which is shown to be the key parameter for this problem. This section also defines our concept of a well-resolved simulation. Section 5 follows up to investigate the accretion rate as a function of smoothing length, as well as opacities for well-resolved simulations. Section 6 discusses the envelope structure of a few selected simulations to highlight the influence of the simulation parameters and to connect those back to the accretion rate measurements. Finally, Sect. 7 discusses the implications of this study and the various findings from our simulation runs. In Appendices A–C we expand on the main text by justifying the flux-limited diffusion approximation and luminosity measurements, and derive the potential temperature, a tool we use in this paper.

2 Physical model

The modelling of gas accretion rates requires the treatment of density, momentum, and two energy conservation equations. In this section we present those, and discuss a few of their important properties.

2.1 The Navier–Stokes equations

We use the standard, time-dependent system of compressional mass conservation and Navier–Stokes equations for viscous fluids $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho u\right)=0,\\ \frac{\partial u}{\partial t}+u\cdot \nabla u=-\frac{1}{\rho }\nabla P+\nabla \cdot \Pi +g,\end{array}$$

where t is the time, ρ is the scalar volume mass density of the gas, u is the local vectorial velocity field, P is the pressure scalar field, and Π is the viscous stress tensor, defined as $$\Pi =2\rho \nu \left[D-\frac{1}{3}(\nabla \cdot u)I\right]\text{},$$(3)

which is a function of density, physical viscosity ν, shear tensor $D_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$, and velocity divergence times the unit matrix I. Finally, g is the local vectorial gravity field.The field of gravity g contains thecontributions from the star as well as the planetary core. The planetary gravity field is additionally smoothed out in the innermost grid cells. Here we use the gravitational smoothing prescription from Klahr & Kley (2006), which sets the planetary gravitational potential to $${\Phi }_{\text{p}}=-m_{\text{p}}G\left[\frac{r^3}{r_{\text{s}}^4}-2\frac{r^2}{r_{\text{s}}^3}+\frac{2}{r_{\text{s}}}\right],$$(4)

where g = −∇Φ(r), r is the planetocentric radius, m_p the planetary mass, G the gravitational constant, and r_s is the smoothing length. At this point we also introduce the dimensionless gravitational smoothing length ${\tilde{r}}_{\text{s}}=r_{\text{s}}/r_{\text{H}}$, where r_H is the planetary Hill-radius, $$r_{\text{H}}=r_{\text{p}}{\left(\frac{m_{\text{p}}}{3m_{\odot }}\right)}^{1/3},$$(5)

with r_p being the planet-to-star semi-major axis and m_⊙ the stellar mass. Gravitational smoothing is necessary to reach numerically convergent solutions when increasing the simulation resolution. For works that study migration (i.e. Kley et al. 2009), gravitational smoothing additionally takes care of avoiding singular values for gravity that could occur when the planet migrates very close to a cell midpoint. We neglect momentum transfer forces from the photons onto the gas, as their ratio with respect to planetary gravity is only ~ 10⁻⁵ in the worst case of the highest opacities used in our model.

Equations (1) and (2) are solved in FARGOCA (Lega et al. 2014) via a classic operator split method (Stone & Norman 1992). More information on this methodology can be found in Masset (2000). As a frame of reference we chose the planetary co-rotating frame, while the planetary radial position is fixed at r_p = 5.2 AU distance from the star. Inertial forces resulting from this transformation are the Coriolis and centrifugal force and are not written explicitly into Eq. (2).

2.2 Energy equations

Our simulations solve two energy equations, one for the internal energy of the gas, e_int, and one for the radiative energy of photons, e_rad. Technically, this corresponds also to two temperatures, one for the gas and one for the radiation. We nevertheless only name one of those by name, the gas temperature or simply the temperature T, which connects to the internal energy over the equation of state $P=\left(\gamma -1\right)e_{\text{int}}$, with the constantadiabatic index γ set to 1.4 in our simulations. While γ can vary in a fully realistic setting, we restrict our simulations to potential depths and thus temperatures which would vary γ only weakly(Popovas & Jørgensen 2016). Only our deepest potentials reach T ≈ 2000 K and only in the few grid cells closest to the planet. Using a variable γ would then change our results only negligibly (O. Gressel, priv. comm.).

The two energy equations are (Commercon et al. 2011) $$\begin{array}{l}\frac{\partial e_{\text{int}}}{\partial t}+\left(u\cdot \nabla \right)e_{\text{int}}=-e_{\text{int}}\nabla \cdot u-\rho {\kappa }_{\text{P}}\left(B-c{e}_{\text{rad}}\right)+Q_{\text{visc}},\\ \frac{\partial e_{\text{rad}}}{\partial t}+\nabla \cdot F=\rho {\kappa }_{\text{P}}\left(B-c{e}_{\text{rad}}\right)\text{},\end{array}$$

where κ_P is the Planck-mean opacity, B is the black-body energy radiation density with B = 4σT⁴, σ is the radiation constant, c is the speed of light, Q_visc = νρΠ² is the self-contraction of the stress tensor, and F is the radiative flux density. Mechanisms for local change of internal energy correspond to the processes of advection, compression, coupling to the radiation field, and viscous heating, from left to right in Eq. (6).

The radiative fluxes F are computed via the classical flux-limited diffusion (FLD), which is a formalism introduced by Levermore & Pomraning (1981) that presses the radiation energy fluxes in the optically thin limit in a form compatible with the optically thick limit, in order to allow for solutions in both cases with only one solver. Appendix A shows how the fluxes transition from the optically thick to the optically thin limit. In general, the FLD radiative fluxes have the form $$F=\frac{\lambda c}{\rho {\kappa }_{\text{R}}}\nabla e_{\text{rad}},$$(8)

with κ_R being the Rosseland mean opacity and the interpolation function λ being a function λ = f(ρκ_R∕L_e) of the local optical depth per grid cell divided by the local radiative energy gradient scale length L_e. Here, λ is calculated for each cell interface in order to decide between optically thick or thin fluxes. Specific choices for κ_P and κ_R are discussed in the following section. The implicit FLD solver that we use is based on the work of Commercon et al. (2011), also described there in great detail, and has already been successfully used for example in Bitsch et al. (2013) and Lega et al. (2015).

A useful reformulation of the internal energy equation into that of the specific entropy S is (see Vallis 2006, Chap. 1.6) $$\rho \frac{\text{d}S}{\text{d}t}=-\rho {\kappa }_{\text{P}}\left(B-c{e}_{\text{rad}}\right)+Q_{\text{visc}},$$(9)

where d∕dt = ∂∕∂t +v ⋅∇ is the Lagrangian time-differential operator. This equation shows that specific entropy can only be generated through viscous heating or is lost through coupling to radiation. The specific entropy is directly proportional to a diagnostic quantity often used in atmospheric physics, namely the potential temperature ϑ (see e.g. Vallis 2006) and it is S = c_P ln(ϑ) for a gas of constant specific heat capacity at constant pressure c_P. The potential temperature can be thought of as a temperature corrected for effects of compression or expansion and is therefore a quantity derived from the simulation data, not an additional temperature that we solve for. Just like the entropy, ϑ remains constant in regions of active convection, and has a positive gradient versus gravity in radiative regions.

As we show in Appendix C, the potential temperature can, for an arbitrary reference density ρ₀ and a constant adiabatic index of γ = 1.4 which we use, be computed as $$\vartheta =T{\left(\rho /{\rho }_0\right)}^{-2/5}\text{}.$$(10)

It is therefore not only a quantity that can be used to assess dynamical instability, but is also a useful one to compare envelope structures to one another, without the need to inspect both density and temperature.

2.3 Physical parameters – viscosity and opacity

In this section we outline our parameter choices for the viscosity, which determines the disc structure, and opacity, which controls the cooling of the planetary envelope and the disc structure. Another important heating source, the compressional work, does not depend on the choice of any free parameters and is therefore not discussed here.

The physical viscosity ν is kept throughout our work at ν = 10⁻⁵ in code units (c.u.), wherethe planetary distance r_p and the Keplerian speed of the planet v_K are both unity. This corresponds for a planet at 5.2 AU around a solar-mass star to a physical viscosity of $\nu ={10}^{15}{\text{cm}}^2{\text{s}}^{-1}$. With the standard relation between physical and α-viscosity of ν = αH²Ω, our simulation values of H∕r = 0.021 … 0.036, correspond to α = 3.3 … 1.1 × 10⁻² at a distance of the orbit of Jupiter. Our adopted nominal viscosity value is thus comparatively large. We plan to study lower values of the turbulent viscosity in a future publication. It is likely that the angular momentum in protoplanetary discs is transported by laminar stresses due to the large-scale magnetic field that penetrates the disc (e.g. Turner et al. 2014). Some works therefore adopt a two-alpha approach, with a disc α setting the angular momentum transport and a turbulent α of much lower magnitude setting the turbulent viscosity (e.g. Ida et al. 2018). However for reasons of simplicity we adopt a constant ν model here.

In order to obtain correct post-shock entropies, we employ a standard method for capturing shocks (Von Neumann & Richtmyer 1950) that uses an artificial viscosity to smear out any shock over $\sqrt{2}$ cells. In our work, this is mainly relevant for the spiral arm shocks triggered in the disc by the planet. Because our envelopes feature quite substantial pressure support radially as well as vertically, no accretion shock is established for the investigated parameter values.

The opacity, however, is the main physical parameter of interest for this study. From Eqs. (7) and (8) it is evident that the Rosseland-mean opacity κ_R controls the diffusive properties of radiative energy. It is this opacity that determines whether the transition from one cell to the next is optically thick or thin. The Planck mean opacity κ_P serves to couple internal and radiative energy. If κ_P > κ_R in an optically thin simulation region, then the sourcing term in Eq. (7) proportional to ρκ_P will be larger than ρκ_R and the gas will lose internal energy more readily. In an optically thick simulation region, κ_P ≠ κ_R will not matter, as the sourcing term is zero anyway in equilibrium.

The two different opacities κ_P and κ_R from Eqs. (7) and (8) are taken from one of the following sets:

• 1.

  Constant opacities κ_P = κ_R ∈{1.0, 0.1, 0.01} cm² g⁻¹. Constant opacities are probably not realized in nature, but provide a basis for a relatively controlled environment to study the envelope structure and accretion process.

• 2.

  The Bell and Lin opacities (Bell & Lin 1994). Here the opacities are ${\kappa }_{\text{P}}={\kappa }_{\text{R}}=\varepsilon {\tilde{\kappa }}_{\text{BL}}$ and ɛ ∈{1.0, 0.1, 0.01}. A cut through the function ${\tilde{\kappa }}_{\text{BL}}={\tilde{\kappa }}_{\text{BL}}(\rho ,T)$ is plotted in red in Fig. 1.

• 3.

  The Malygin opacities from Malygin et al. (2014). These include updated Planck gas opacities together with the dust opacitiesfrom Semenov et al. (2003). We use ${\kappa }_{\text{P},\text{SM}}=\varepsilon {\tilde{\kappa }}_{\text{P},\text{SM}}\ne {\kappa }_{\text{R},\text{SM}}=\varepsilon {\tilde{\kappa }}_{\text{R},\text{SM}}$ and ɛ = 0.01. A cut through the functions ${\tilde{\kappa }}_{\text{P},\text{SM}}(\rho ,T),{\tilde{\kappa }}_{\text{R},\text{SM}}(\rho ,T)$ is plotted in blue in Fig. 1.

For the latter two sets, where the opacity has a non-constant form, one can think of the parameter ɛ as a reduction factor relative to the ISM opacities, such that ɛ ∈{1.0, 0.1, 0.01} correspond to {1.0, 0.1, 0.01} ISM opacities and thus dust-to-gas ratios of {1.0, 0.1, 0.01}%. The order of magnitude of the constant opacity set and the opacities from Bell & Lin (1994) then coincide roughly at the water-sublimation iceline. In set (3), at higher temperatures one finds a significant deviation of κ_P from earlier estimates (see Fig. 1). Malygin et al. (2014) show that this is because κ_P is mostly sensitive to the spectral line wings of the underlying sampled linelist, and this was undersampled in earlier work. The main simulations in our work have been performed using the first two opacity sets, while the third set serves mainly to inform us about the simulation behaviour when κ_R≠κ_P.

Although in the third opacity set the difference between κ_P and κ_R becomes very large above T ~2000 K, in our work its influence is negligible, as our limitations in resolution prohibit gravitational potentials that are deep enough for those temperatures. In general however, one would expect there to be a visible effect when κ_P ≠ κ_R: the function of κ_P is to couple the evolution of e_int and e_rad, while the fluxes that cool the envelope are a function of κ_R. The same is true for the optical thin–thick transition, or the τ = 1-surface, defined by $-{\displaystyle {\int }_{\infty }^z\rho }{\kappa }_{\text{R}}\text{d}z=1$. If κ_P = κ_R, then at the same height where τ = 1, the coupling between radiation and matter also stops being efficient. However, when κ_P > κ_R, matter keeps transferring internal energy into photons even above the τ = 1-surface. This allows for more efficient cooling and therefore higher values of accretion rates are expected for set (3) compared to set (2).

The most prominent differences between the two non-constant sets can be seen in the potential temperature structure of the envelope (see Fig. 13). The multibump structure in the curve for the Malygin et al. (2014) opacities is due to the multitude of chemical ingredients used to compute the opacities, each having their own sublimation temperature. The opacity jump then results in differing temperature gradients, which impact the potential temperature.

Fig. 1 Opacity as a function of temperature. Constant opacities take the three different values 1.0, 0.1, and 0.01 cm2 g−1, where in the plot we indicate only 1.0 cm2 g−1. The second opacity set from Bell & Lin (1994) increases with temperature, shows a transition at the water-iceline at ≈ 200 K and a steep decline of the opacity once dust sublimation become relevant at temperatures of 1000 K. The third opacity set from Malygin et al. (2014) shows a much richer structure in terms of opacity transitions, which is due mainly to the various chemical constituents that Semenov et al. (2003) computed. The large difference between the Planck opacity and the Rosseland opacity mainly plays a role in optically thin regions. Midplane temperatures in our simulations at 5.2 AU without a planet range from 30 to 50 K.

3 Numerical procedures and parameters

Here we outline the most important properties of the code used and the individual simulation steps. We employ three such steps (Lambrechts & Lega 2017), which separate physical effects according to timescales and ascertain that our final results are not influenced by, for example, the process of forming gaps or spiral arms. We also present the simulation runs that we performed, characterised by their parameter values. The parameters that we varied are the opacity sets, the dimensionless gravitational smoothing length ${\tilde{r}}_{\text{s}}=r_{\text{s}}/r_{\text{H}}$, and the number of grid cells per Hill-radius N_c. We also use N_c as a measure for the numerical resolution in the simulation domain.

The steps we employ are connected via restarting and interpolating onto a new grid. Step A consists of a 2D radial–vertical run without a planet to find the radiative equilibrium structure of the disc, which gets its structure from viscous heating alone. Step B restarts from the end of Step A in a full-3D low-resolution grid of the same domain dimensions and plugs in the planet smoothly in order to generate the planetary gap and spiral waves. Once the gap depth satisfies a convergence criterion, we cut a global ring out of the disc and use this to run a 3D high-resolution grid on it to measure accretion rates, which is Step C.

3.1 Numerical parameters – grid, resolution, and smoothing length

We employ a grid in spherical coordinates centred on the star, [r, ϕ, θ], where r is the radial distance from the star, ϕ is the colatitudinal angle, and θ is the azimuthal angle, varying along the disc orbit. This results in quasi-Cartesian coordinates surrounding any single point of interest around the planetary orbit, when sufficiently high resolution around this point is used. For flows orbiting the star this is a natural coordinate frame, since it minimizes numerical diffusion. However, flows that orbit a planet in the disc, and thus in a Cartesian region, are then subject to strong numerical diffusion. To remedy this problem, Lambrechts & Lega (2017)introduced non-uniform grid spacing between cells along each spherical coordinate in order to resolve the Hill-radii of low-mass protoplanets, similar to the grid with quadratic grid-spacing used in Fung et al. (2015). In this study, we introduce a modified version of this grid in the azimuthal direction that follows a Gaussian grid-spacing, meaning that we achieve high resolution in the planetary Hill-sphere and uniform grid-spacing at the simulation boundaries.

The application of the quadratic non-uniform grid is implemented in radius and co-latitude with a functional dependence, $$\text{d}r\propto {(r-r_{\text{p}})}^2,$$(11)

and $$\text{d}\varphi \propto {\left(\varphi -\frac{\pi }{2}\right)}^2.$$(12)

In the azimuthal direction we use two different strategies: Equidistant spacing for moderate resolutions of up to 65 grid cells per Hill-radius to allow the application of the FARGO algorithm for the orbital advection around the star. For a number of runs that correspond to our highest achievable resolution, namely 100 and 150 grid cells per Hill-radius, we apply a non-uniform grid in azimuth. For the latter cases, we applied azimuthal grids that vary Gaussian with quasi-equidistant spacing d θ_qe far away from the planet. The Gaussian grid has the following spacing. $$\text{d}\theta \propto (\text{d}{\theta }_{\text{p}}-\text{d}{\theta }_{\text{qe}})\left\{\exp \left[-{(\theta -{\theta }_{\text{p}})}^2/W\right]-1\right\}+\text{d}{\theta }_{\text{p}},$$(13)

with θ_p being the azimuthal position of the planet, dθ_p the azimuthal resolution at the position of the planet, and W the Gaussian width that is found in accordance with requiring a certain number of total cells.

The simulation domain is global in azimuth for all runs, meaning that θ ∈ [−π;+π]. The extent of the simulation box in co-latitude is chosen such that it represents ≈ 3r_H in physical space which results in ϕ ∈ [90°;83°]. We do not extend the simulation box below the midplane, but use a half-disc approach instead that mirrors the hydrodynamic quantities in the cell above the midplane into the first ghost cell below it. Thus, we assume symmetric physics in all variables below and above the midplane in order to speed up our code.

Radially, we choose a ring of r ∈ [0.4;2.5] ⋅ 5.2 AU for the low-resolution Steps A and B and r ∈ [0.7;1.3] ⋅ 5.2 AU for the high-resolution Step C, and the planet is kept fixed at r_p = 5.2 AU. Numerical resolution – the number of grid cells between the box limits and their individual spacing – varies between the simulation steps and is therefore documented individually below.

The temperature near the planet in Step C (discussed in Sect. 3.4 below) where the accretion takes place is mainly regulated by the opacity and the gravitational smoothing length r_s (Klahr & Kley 2006). In our study, the latter takes values of r_s ∈ [0.5;0.05] r_H and we define a convenient quantity, the Hill-sphere normalized smoothing length ${\tilde{r}}_{\text{s}}=r_{\text{s}}/r_{\text{H}}$.

3.2 Step A – heating/cooling equilibrium for the disc alone in 2D

The first step places an initial mass distribution of $\Sigma =120\text{g}{\text{cm}}^{-2}{(r/5.2\text{AU})}^{-0.5}$ with Σ₀ = 120 g cm⁻² being the gas mass surface density at 5.2 AU and Gaussian vertical structure in the simulation domain. According to classical literature (Frank et al. 2002, Eq. (5.9)), that is, v_r = 3ν∕(Σr^1∕2)∂_r(Σr^1∕2), this is the density profile required for zero accretion, which was also used previously in Bitsch et al. (2013). The resolution is N_r × N_ϕ = 288 × 32 with equidistant spacing in both dimensions. No planet is present in this step. Thereafter, the density and temperature structure of the disc is self-consistently calculated via solving Eqs. (1)–(7) and evolved until a steady state is reached. This steady-statesolution is then taken as the initial condition for Step B.

3.3 Step B – gap opening and spiral arms in low-resolution 3D

With the initial conditions being the end-state from Step A, we extend the disc to 3D with the azimuth θ ∈ [−π;+π], N_θ = 768 and equidistant angular spacing, and plant a seed planet into it. The planetary potential is centred at the corners of the eight grid cells around the point (r, ϕ, θ) = (5.2 AU, 90°, 0°) and is ramped up slowly over ten orbits from zero to its final mass of 90 m_⊕ with ${\tilde{r}}_{\text{s}}=0.5$ in order to avoid unphysical shocks in the disc and the planetary envelope. Subsequently, we let the planet continue for 400 orbits until a steady-state gap has been opened. For our relatively large viscosity the 90 m_⊕ planet opens a partial gap of about 30% of the unperturbed gap depth (see Fig. 3).

We compare our gap to the semi-analytical predictions from Crida et al. (2006), where we have taken the value of the gap depth after half a viscous timet = t_ν∕2 and $t_{\nu }=r_{\text{p}}^2/\nu ={10}^{4.2}{\Omega }^{-1}$. Our viscosity corresponds to a Reynolds number of Re ≡ t_νΩ = 10^4.2 and with this we can compute the expected gap depth from the Crida-formula, resulting in a predicted 40% gap. It is interesting that our gaps are slightly deeper for the same h∕r, possibly hinting at some additional physical processes contributing to the gap formation process.

In Fig. 2 we show the aspect ratio of the disc as a function of the distance to the star for the constant opacity cases. As the locally isothermal midplance scale height H depends on the local temperature T through H = k_B T∕μg, it is evident that one finds colder temperatures with decreasing opacity in the optically thick disc midplane. For κ = 0.01 cm² g⁻¹, the planetary gap becomes optically thin and is then directly heated up by the gap walls. Thus, the material inside the gap becomes hotter for κ = 0.01 cm² g⁻¹ than for κ = 0.1 cm² g⁻¹. Surface densities follow the temperature rather than the opacities, which indicates the importance of radiative transfer for predicting correct gap depths.

We use the end state of this simulation step again as initial conditions for Step C. Therefore, it is important to check whether the planetary gaps are converged before we do so. We confirm this in Fig. 3.

Fig. 2 Profiles of the average aspect ratio in the disc midplane for runs with constant opacity after 400 orbits in Step B. The planet is kept fixed at 5.2 AU. The three different opacities that we use for later accretion rate measurements are used as running parameter. For κ = 0.01cm2 g−1 the gap becomes optically thin, and is then directly heated by the gap walls, which explains the higher gap temperature compared to κ = 0.1 cm2 g−1.

3.4 Step C – high-resolution step in 3D

The previous steps ascertain that the disc around the planet is in radiative equilibrium once we apply the non-uniform grid (see Sect. 3.1). To limit computational cost, we now shrink the simulation domain in radius to r ∈ [0.7;1.3] ⋅ 5.2 AU. We keep the information in the physical variables from Step B and apply boundary conditions where the boundaries are relaxed to the values obtained in Step B. This approach is valid if the new boundaries are far enough from the planet to be unaffected by the physics inside the Hill-sphere. The planet sits at r_p = 5.2 AU and has a Hill-radius of r_H = 0.25 AU, thus the simulation domain radially extends ~ 6 r_H to both sides, and about 3 r_H in colatitude.

We then interpolate the data from Step B onto the Step C non-uniform grid, and deepen the potential depth during a time of t_ramp = 2 Ω⁻¹ from r_s = 0.5 r_H to smaller values of r_s = 0.2, 0.1, 0.05 r_H.

Fig. 3 Gap depth vs. time and convergence ratio: for three different masses with all other numerical and physical parameters being constant, we compare the resulting gap depth as a function of orbits and the convergence measure (Σ(t) −Σ(t + dt))∕Σ(t). We see that after 400 orbits the gap depth has essentially reached its final value (up to 10–20%). The convergence rates for the gap depths are independent of opacity after a few orbits and thus their relative depth will not change. Different opacities behave identical.

3.5 Summary of simulation parameters

After explaining the simulation steps in detail, we briefly list the most important numerical parameters in Table 1. When Step C is reached, accretion rates can be measured; we explain this in more detail in the following section. Various simulation runs have been performed in order to ascertain numerical convergence and investigate parameter dependencies of those accretion rates. We list these dependencies in Table 2.

3.6 Measurement of mass accretion rates

The central aim of this study is the measurement of accretion rates in Step C and to obtain an understanding of the influence of resolution, gravitational smoothing length, and opacity on the results. The influence of those numerical parameters will be presented inthe results section, and we now turn to describing the two methods with which the mass accretion is calculated.

1. Direct measurement of mass per snapshot in time. Two successive snapshots can be subtracted to calculate a mass accretion rate: $$\dot{m}=\frac{m_{\text{H}}^{n+1}-m_{\text{H}}^n}{t^{n+1}-t^n},$$(14)

where $m_{\text{H}}^n$ is the mass inside the Hill-sphere at time-step n.

2. Compute a net mass accretion rate via average mass fluxes in and out of the Hill-sphere from only one snapshot. With this, we get the accretion rate $$\dot{m}=-{\displaystyle {\int }_{\partial \Omega }\text{d}}n\cdot u\rho ,$$(15)

where dn is the normal vector on the Hill-surface ∂Ω. We find that both these methods agree reasonably well within 10%. Therefore, to simplify the presentation of our results, we only show results given by Eq. (14). We now proceed to introduce general results from the measurements of mass accretion rates and their dependencies on numerical parameters.

Our disc accretion rate on the other hand is set to zero. This is implemented self-consistently by choosing a surface density profile of Σ ~ r^−0.5 in Step A that corresponds to an equilibrium disc.

4 Dependence of the mass accretion rate on the numerical resolution

This section gives an overview of the accretion rates that we measure in Step C, without going into considerable physical detail at first; we are initially interested in a clear picture of how numerical parameters influence our results. We then present the dependence of the accretion physics on the opacity and interpretations of the flow structure inside the Hill-sphere in the following two sections. Figure 4 shows an overview of the structure of density ρ, stream lines, and temperature T in the midplane. Typical features for gas flows around gas giants are visible: the gas gap along the azimuth, spiral arm shocks, and horseshoe flows feeding gas into the planetary envelope. Two flow separatrices define an inner and outer envelope region. One might expect there to be a simple relationship of the outer region feeding the inner one, but the 3D structure is more complicated than that, with some flows coming from the vertical direction. The temperature distribution is mostly determined by the compressional heating from shocks and accretion balanced by radiative cooling and weak heat advection effects.

To help interpret the mass accretion rates we measure the potential temperature ϑ, as defined in Eq. (10). The potential temperature shows a marked drop in the central regions of the Hill-sphere. This is a result of entropy reduction by radiative cooling and is clear evidence that the gas envelope is slowly accreted in our simulations.

When measuring the increase in mass in the planetary envelopes we first noted a strong dependence of the accreted mass on resolution. This dependency led us to perform a convergence study, which we describe first before discussing the physical state of the envelope further. Here, it became clear that underresolving the smoothing length, not the Hill-sphere, is the parameter responsible for this strong dependence. At low resolution of the smoothing length, one finds only weak entropy gradients in the envelope. The envelope then quickly finds a steady state and constant mass. In this steady state, the envelope still features a strong flow of horseshoe gas through itself, but none of this gas actually gets accreted.

Only when increasing the resolution does the entropy drop significantly, and at this point a part of the mass from the horseshoe gas starts to be accreted. Finally, at sufficiently high resolution we find that the mass in the envelope and the accretion rate converge, as one would expect from a physical result.

Interestingly, we have found that an underresolved simulation in steady state can be restarted at higher resolution, whereafter that simulation will be accreting again. Therefore, we conclude that (1) the steady states found must be a numerical artifact and (2) while the decrease in entropy can be used to describe the accretion process, the entropy loss is not guaranteed, as it is instead a change in heating and cooling physics at high resolution that allows the entropy to drop further, allowing accretion to continue.

Figure 5 shows that the mass in the envelope compared to the core mass (90 m_⊕) is small at all times. The initial value of the envelope mass is a result of the midplane density in the disc into which we relax the planetary core mass. Growth continues only for a short time of a few tens of orbits, as this is limited by our computation time. As m_env ≪ m_p our assumption to neglect self-gravity in our model is validated for all run times. The neglect of self-gravity in our model is thus self-consistent inside this model, but a realistic core of 90 m_⊕ would probably contribute a significant portion of this mass to the inner envelope, which we do not resolve. We discuss this further in Sect. 6.

A look at the accretion rates in Fig. 6 allows a more detailed study of the numerical convergence behaviour, as well as a comparison with necessary accretion rates to build the solar system gas giants. With increasing resolution, accretion rates quickly converge during approximately the first five orbits, but then rapidly decline if the smoothing length is underresolved. However, even for our best-resolved simulations there is a slight decrease in mass accretion rate, the non-convergence saturating at a level that converges. We discuss this decrease in more detail in Sect. 6. If taken at face value, the accretion rates that we see in Fig. 6 would be enough to form a Jupiter-mass planet in 1% of the nominal disc lifetime of 3 Myr. This is further discussed in Sect. 7.

In order to investigate whether the gas accretion into the envelope is irreversible, one needs to compute the evolution of the specific entropy of the gas. As introduced earlier, we investigate the related potential temperature instead, which is simple to compute directly from our simulation data. We plot the central potential temperature ϑ_c in Fig. 7, where we emphasise that $\dot{\vartheta }>0$ can only happen through the action of viscosity, and $\dot{\vartheta }<0$ indicates cooling (see our earlier mention of Eq. (9)), possibly with a small positive contribution from viscosity. We find ${\dot{\vartheta }}_{\text{central}}<0$, and therefore the mass accretion is irreversible and limited by cooling. Interestingly, the evolution of ϑ_c with increasing resolution indicates that there are still physical effects that are insufficiently resolved, although the accretion rates converge. The central temperatures converge as well, thus the behaviour of ϑ_c indicates that ρ_c increases with resolution, or in other words, the mass inside the envelope concentrates more and more in the centre as resolution increases.

Fig. 4 Overview and snapshot of typical midplane variables during gas accretion in Step C at the time of ten orbits for our nominal simulation run with constant opacity of κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$, and Nc= 100. The Hill surface of the planet is approximated as a sphere, which is marked with the red circle. Top left: gas mass density in the midplane and stream lines overplotted with their respective momentary velocities. Spiral arm shocks launched near the Hill-sphere of the planet are clearly visible. The horseshoe flows interact with the shocks and are partially deflected. Colours indicate the stream line velocities and highlight flow deceleration at the shocks. Judging from stream lines only it would be impossible to say whether there is ongoing gas accretion, because most of the horseshoe mass merely passes through the planetary envelope, while only a minor fraction is accreted. The planetary gap is approximately ~ 2rH in width (as wide as the plot) and attains its unperturbed disc value at ~ 3rH. For comparison, the simulation domain in Step C extends radially to ±6rH. Top right: temperature structure in the midplane. The cold gas in the gap and the surrounding disc is clearly visible as dark background colour, while the planetary envelope and beyond is heated by compressional heat from accretion and the spiral shocks. A small fraction of internal energy is advected along the horseshoe orbits, visible in separate moderately warm (~ 80 K) arm structures marked with arrows. Bottom left: potential temperature ϑ overlaid with iso-ϑ contours. Entropy generation at the spiral arm shocks is strongest where horseshoe orbits hit the envelope radially. The central region cools most efficiently in the vertical direction, and therefore develops a drop in ϑ, indicating a radiative region. The innermost few pixels are subject to further flattening of the ϑ-profile, as herethe smoothing length provides a region essentially free of gravity. Bottom right: cylindric average of density as background colour and temperature in contours, as well as the τ = 1 surface, denoting the optically thin-thick transition. As this run has a constant opacity, the τ = 1 surface follows essentially an iso-density surface. This surface is bent towards the planet at small r, due to heavy cooling into the disc atmosphere and the consequential contraction. This plot only extends out to 1 rH as otherwise averaging effects smear out the physical properties outside the Hill-sphere.

5 Dependence on smoothing length and opacity

Besides the numerical resolution, there is also the gravitational smoothing length r_s as à free numerical parameter. In this section, we investigate the influence of the value of r_s on the accretion rates. When varying the smoothing length it is necessary to clarify how much resolution is needed for convergent results. After all, it is possible that more resolution is needed than in the cases presented before, as the gravitational potential is deepened. After analysing the dependence of the measured gas accretion rates on the resolution of the smoothing length, we use the numerical parameters that satisfy convergence in order to investigate the accretion behaviour when changing the opacity.

Fig. 5 Total envelope mass within the Hill-radius as a function of time, immediately after the potential ramp-up time of two orbits, plotted against the numerical resolution of the envelope. We show the runs with cells per Hill-radius Nc ∈ {25, 35, 50, 65, 100} and smoothing length ${\tilde{r}}_{\text{s}}=0.2$ for a constant opacity of κ = 0.01cm2 g−1. Under-resolving the Hill-radius leads to a gradual decrease in the measured accretion rate, which results in the envelope as a whole achieving a dynamic steady state. Increasing the resolution step by step we find numerical convergence of accretion rates and a less pronounced drop in accretion rates over time. We note that our highest resolution (Nc = 100) drops slightly in accretion rate compared to the case with Nc = 65.

Fig. 6 Mass accretion rate as a function of time. We present results of runs with cells per Hill-radius Nc ∈ {25, 35, 50, 65, 100} and smoothing length ${\tilde{r}}_{\text{s}}=0.2$ for an opacity of κ = 0.01 cm2 g−1. These values are two orders of magnitude higher than the necessary accretion rate to form Jupiter within three million years, which is approximately Ṁ = 10−4 m⊕ yr−1.

Fig. 7 Central potential temperature for the same simulation parameters as in Figs. 5 and 6, κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$. The steady decrease in potential temperature confirms that the accretion is final and irreversible. Just to be sure, we checked andconfirmed that the entropy in the rest of the envelope decreases as well.

5.1 Convergence criterion

When performing simulations with smaller smoothing length (see Fig. 8), we note that in fact an even higher resolution is needed before accretion rates converge. This indicates that it is r_s and not r_H that is the relevant length scale to be resolved. In order to quantify the resolution required to achieve convergence, we plot the values of the accretion rates averaged over ten orbits with ten samples per orbit against the number of cells per smoothing length. The result can be seen in Fig. 8. There, the accretion rate can be seen to first rise almost linearly with increased resolution for all r_s except for ${\tilde{r}}_{\text{s}}=0.5$, before leveling off after r_s has been resolved with ≥10 cells. We consider the runs with ${\tilde{r}}_{\text{s}}=0.5$ as anomalous, as the force-free region inside r_s there connects directly with the horseshoe flows, giving rise to strange velocity shears that otherwise would not exist. The value towards which the simulations converge is slightly different for the individual values of ${\tilde{r}}_{\text{s}}$ (indicated in Fig. 8 by the dashed blue and dashed red lines). Although small, we deem this difference real, as there are physical differences in the envelopes that are introduced by deepening the potential. We discuss this further in Sect. 6.

From Fig. 8 we conclude that we need at least ten cells per smoothing length (Lambrechts & Lega 2017 already argued for >5) in order to obtain convergent accretion rates. Due to runtime considerations this limits our minimum achievable smoothing length to ${\tilde{r}}_{\text{s}}=0.05$ and with this parameter we unfortunately do not achieve full convergence, but we can already observe a trend that is similar to that seen in simulations of larger smoothing length. The physics that need to be resolved by the resolution of those ten grid cells are also discussed in Sect. 6.

5.2 Trends with smoothing length and opacity once converged

We now analyse the behaviour of simulations at the highest available resolutions as a function of the smoothing length and opacities. Those simulations are either fully resolved or very close to being so, and therefore we can ignore effects of resolution in this section. The first set of simulations we introduce here are the constant opacity runs. The second and third sets are those with non-constant Bell & Lin (1994) and Malygin et al. (2014) opacities. We show the measured mass accretion rates at the highest resolutions in Fig. 9.

Our simulation runs with values of ${\tilde{r}}_{\text{s}}=0.2,0.1,0.05$ demonstrate that there is a coherent decrease of accretion rates when deepening the gravitational potential for all constant opacities as well as for more complex opacity functions. The runs with ${\tilde{r}}_{\text{s}}=0.05$ are marginally well-resolved according to our criterion, and it is conceivable that the converged value for this parameter for ṁ is slightly higher. However, our simulation runs are currently limited to the maximal resolution of N_c = 150. Furthermore, those runs seem to agree with the trend set by ${\tilde{r}}_{\text{s}}=0.2,0.1$. We therefore see no reason to discard these data points as insufficiently converged and we use them for the extrapolation of the trend of ṁ versus r_s. The run with ${\tilde{r}}_{\text{s}}=0.05$ and κ = 1.0 cm² g⁻¹ is discussedseparately in Sect. 7.2. Finally, the runs with ${\tilde{r}}_{\text{s}}=0.05$ and complex opacities showed stronger irregularities and variability than in their constant-opacity counterparts. Strong oscillations prohibited us from executing those runs, and correct execution probably requires an even higher resolution than we can achieve here.

The overall trend of decreasing ṁ with ${\tilde{r}}_{\text{s}}$ can be puzzling at first: A smaller gravitational smoothing length means stronger gravitational forces are acting on the gas close to the planet, which one would associate with a higher accretion rate. However, a deepened gravitational potential means first and foremost that there is more gravitational energy available, which is used to heat the gas. When deepening the potential, gravitational forces quickly compress the gas adiabatically, rising temperatures and increasing pressure support.

Nevertheless, only a fraction of the additional injected internal energy is able to diffuse outward as radiative energy, which results in a disproportionately small increase in envelope luminosity. When going from ${\tilde{r}}_{\text{s}}=0.2$ to ${\tilde{r}}_{\text{s}}=0.1$ and thereafter to ${\tilde{r}}_{\text{s}}=0.05$, we measure a luminosity increase each time bya factor of ~1.5, while internal energy and central temperatures rise by a factor of ~2.0, proportional to $\text{GM}/{\tilde{r}}_{\text{s}}$. With the help of Eq. (16) one can then post-predict the accretion rates and see that ṁ will indeed decline.

Another interesting feature is the departure of the simulation run $({\tilde{r}}_{\text{s}},\kappa )=(0.05,1.0{\text{cm}}^2{\text{g}}^{-1})$ from the trend seen for the preceeding κ = 1.0 cm² g⁻¹ (black triangles in Fig. 9). It is conceivable that this might be a simple effect of under-resolving this particular simulation, or that high opacities simply need more resolution than low opacities. Upon testing this, we could however not find a difference in the resolution trend seen in Fig. 8 which was done for κ = 0.01cm² g⁻¹ and those performed at 1.0 cm² g⁻¹. We conclude that this is not an effect of underresolving that particular simulation run, but is indeed a physical effect. This is further discussed in Sect. 7.2.

The differences between the constant opacities and the more advanced Bell and Lin and Malygin opacities are not very dramatic. Neither the downwards trend with lowering ${\tilde{r}}_{\text{s}}$ nor the spacing when changing the opacity prefactors seems to make a huge difference in the general behaviour. We analyse the envelope structure for different choices of opacity law in connection with the accretion rates in more detail in the following section.

Fig. 8 Our complete set of all simulation runs for constant κ = 0.01 cm2 g−1 where we varied the smoothing lengths as well as the numerical resolution, plotted as accretion rate vs. smoothing length (left) and as accretion rate vs. cells per smoothing length (right). All accretion rates shown are averages of the firstten accreting orbits after increasing the gravitational potential in Step C; see Figs. 5–7. Left: gas accretion rate inside the Hill-sphere as a function of the smoothing length, for varying resolution numbers. Each dot is a simulation run, with the resolution of the Hill-sphere indicated as cell numbers. It is evident that as resolution increases, accretion rates cluster closer together, indicating numerical convergence. However, more resolution is needed for deeper potential wells. Right: gas accretion rate inside the Hill-sphere plotted against the number of cells per smoothing length. Different colours indicate different smoothing lengths. It is evident that we need to resolve the smoothing length with at least ten cells before convergent results can be achieved (vertical lines). Thus, we define a “well-resolved simulation” as one that hits convergent values of Ṁ at ${\tilde{r}}_{\text{s}}\times N_{\text{c}}>10$. Once converged, the accretion rates differ only slightly when deepening the gravitational potential, indicated with dashed lines.

Fig. 9 Dependency of the accretion rates on the smoothing length and opacities for simulations with the highest available resolution as depicted in Fig. 8. Left: averaged accretion rates from 5 to 15 orbits, for ${\tilde{r}}_{\text{s}}=0.5$ down to ${\tilde{r}}_{\text{s}}=0.05$ and our three different opacity values. Simulation runs with ${\tilde{r}}_{\text{s}}=0.5$ have such a shallow potential that the horseshoe flows, which penetrate down to approximately ${\tilde{r}}_{\text{s}}=0.5$ also interfere with the envelope cooling process. Right: accretion rates with Bell and Lin opacities, with 1 (=solar), 0.1, and 0.01% dust opacities, as well as a run with Malygin opacities and ɛ=0.01%. Simulation runs with ${\tilde{r}}_{\text{s}}=0.2,0.1$ have been shown in Fig. 8 to be sufficiently resolved. We therefore interpret the seen downwards trend of Ṁ with ${\tilde{r}}_{\text{s}}$ as real (red dashed line), even if runs with ${\tilde{r}}_{\text{s}}=0.05$ are only at the verge of fulfilling our convergence criterion.

6 Envelope structure

This section aims to link the envelope structure resulting from the simulation parameters N_c, ${\tilde{r}}_{\text{s}}$ and κ, with the measured accretion rates. We use the previously introduced theoretical as well as measured luminosities and explain those with the observed density and temperature structures.

6.1 Comparison to simple 1D models

After establishing numerical convergence, we are interested in understanding whether the solution to which the simulation converges is also the correct physical one.

A popular way to estimate mass accretion rates in 1D models of gas giant formation (see i.e. Mizuno 1980; Pollack et al. 1996)is based on the stellar structure literature and connects the energy loss from the envelope, the luminosity L, with the mass accretion rate via $$L=\frac{\text{GM}\dot{M}}{r}.$$(16)

This relation would be true for slow Kelvin–Helmholtz contraction, and we can check our simulation for consistency with this mechanism. Also, this is a slight improvement over the methodology of Lambrechts & Lega (2017), where the compressional heating was assumed to be identical to the luminosity. From our simulations, we can measure both M and L and check whether they obey relation Eq. (16), to answer the question of which physical process we are investigating. To this end, we estimate the 3D position of the τ = 1-surface, integrating all optical depth elements dzρκ vertically from infinity (the simulation boundary) towards the planetary envelope. An example of the computation of this surface can be seen in Fig. 4. Once this surface is found, the temperature T_{gas, τ=1} on all surface elements is taken, and it is then $L={\displaystyle \sum \Delta }A\sigma T_{\text{gas},\tau =1}^4$, with ΔA being the area of a surface element and σ the radiation constant. Throughout this work we use the unit of a Jupiter-luminosity L_J = 3 × 10²⁴ erg s⁻¹. The example surface in Fig. 4 has a luminosity of 7 × 10⁴ L_J and its evolution in time is shown for a resolved and unresolved simulation in Fig. 10, confirming that it is consistent with the mass accretion rates. This is because the flux-limited diffusion luminosity quickly becomes the free-streaming luminosity F_FS = cE_rad when passing through the strong density gradients of the planetary envelope into the disc atmosphere and the planetary gap. This quick transition is also shown in Appendix A.

This approach is only valid for the simulations for which the τ = 1-surface forms an enclosed surface on the planetary envelope, that is when the planetary gap is optically thin. In fact we find good agreement between the mass accretion rates and the luminosities in those cases. The downwards trend of luminosity also follows the downwards trend for the mass accretion rate, hinting at the increasingly radiatively inefficient envelope. We checked this downwards trend more thoroughly with a full 3D computation of the luminosity based on the actual fluxes in the simulation. This is shown in Appendix B, where we conclude that the actual luminosity is 20% less than the one by our estimate on the τ = 1-surface, but that the two luminosities have a relatively constant ratio.

This is interesting for the purpose of our analysis in that we can now link differences in luminosity to differences in accretion rate behaviour when varying opacity and smoothing length. The connection of luminosity to the envelope structure will be discussed further in connection with the averaged 2D profiles of the envelope.

Fig. 10 Evolution of measured luminosity vs. measured accretion rate for our κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$ simulation run. Sufficiently resolved simulations show behaviour such that the measured luminosity and the measured mass accretion rate correlate well. Unresolved simulations however run into a state of zero accretion rate, but finite luminosity, which remains higher than that of the final luminosity in the resolved simulation. This hints at numerical effects or other physics providing an additional energy source in unresolved simulations.

Fig. 11 Density and temperature structures for constant opacity runs in radial and vertical slices after ten orbits: data are shown for runs with ${\tilde{r}}_{\text{s}}=0.1$ and the nominal run with ${\tilde{r}}_{\text{s}}=0.2$ and Nc= 100 for comparison. A deeper potential raises the whole envelope temperature as long as no convective region is formed. Optically thick regions are hotter for higher opacity (see inset), while the optically thin regions approach isothermal structure and their temperature trend reverses according to the position of the τ = 1 surface (indicated with triangles). The density structure follows simple expectations: a hotter envelope has lower density in optically thick regions.

6.2 Density and temperature

We now present an overview in Fig. 11 of 1D slices of density and temperature profiles of planetary envelopes. Five different simulations are presented, of which three are identical except for the chosen opacity (the solid red, blue, and black curves), twoare identical except for the smoothing length (the solid red and dashed red curves), and another two are identical except for the resolution (the dashed green and dashed red curves). In this way we highlight the physical differences in the envelopes when varying parameters.

6.2.1 Change in resolution

At first, wecompare a well-resolved and under-resolved run (N_c = 100 cells per r_H red dashed with green dashed N_c = 35 cells per r_H) in Fig. 11. Those profiles are only snapshots in time. However, they are representative of the differing evolutionary paths that simulations can take, particularly seen in Fig. 10. The underresolved simulation has a low central temperature, but this has no effect on the accretion rates. It is the potential temperature (see Fig. 14) that helps us to properly interpret the difference between resolved and unresolved simulations: being closer to a full adiabatic profile, the unresolved simulation must be radiatively less efficient and is therefore accreting less. The reasons for this are discussed in detail in Sect. 6.6.

6.2.2 Change in smoothing length

Secondly, we compare a well-resolved deep potential with a well-resolved shallow potential (straight red with dashed red curve) in Fig. 11.

Both runs start with the same initial envelope mass, so naturally for the deeper potential the central density is higher and the outer envelope densities are lower. This is true for the radial structure of the envelopes. This is different in the vertical direction straight above the planet, where at z > 1.0 r_H the deep potential has a higher density than the shallow potential. The density structure at those altitudes, in the transition region from planetary envelope to disc atmosphere, is important because the optically thin–thick transition happens at densities of around ρ ~ 10⁻¹³g cm⁻³ and thus determines the luminosity.

Central temperatures scale linearly with the smoothing length as long as no complex opacity laws are used. This stems from the fact that the energy input into the envelope is ρGM_p∕r_s which is then translated into internal energy c_pρT. The radiative temperature gradient is then followed upwards in z until the optically thin–thick transition is encountered. At this point, photons can escape into space and the temperature slope changes. The remaining temperature gradient comes from the radiative coupling with the disc atmosphere, which decreases in efficiency as ρκ_P; see Eq. (7).

An important point to discuss is the difference in envelope temperatures and luminosities. The deeper potential heats the envelope more efficiently than the shallow one. The resulting luminosities are higher for the deep potential. This does not automatically impart a higher accretion rate however. One can understand the accretion rate according to Eq. (16) as a balance between the luminosity and the energy content of the envelope. For our simulations, the energy content of the envelopes rises as GM ∕r_s, but the resulting luminosities increase only by a factor of ~1.5 for a doubling of the potential depth. Thus, the accretion rates decrease with increasing potential depth.

Fig. 12 Potential temperatures for opacity set 1 with Nc = 100. Deeper potential and lower opacity lead to a flatter potential temperature profile. Radially, the simulation runs with ${\tilde{r}}_{\text{s}}=0.1$ have flat ϑ-profiles; this does not denote the onset of convection, but rather the action of the horseshoe flows that can freely enter and exit the gravitational well. In z, a steeper profile corresponds to stronger accretion flows from this direction. The much higher values of ϑ result from the low density in the disc atmosphere, as we have rH≅Hdisc throughout our work.

6.2.3 Change in opacity

Inside the deep envelope, where the optical depth is larger than one, a higher opacity causes higher temperatures. This is because the radiation diffusion times are proportional to κ, while all simulations of the same smoothing try to radiate away the same initial compressional work. Density is inversely related to temperature in the central region. Temperature gradients in r and z are proportional to 1∕κ while optically thick. However, when following the temperature profiles outwards in z, low opacities hit the optically thin–thick transition earlier and from there on the temperature remains quasi-isothermal with only a weak, opacity-independent temperature gradient. Radially, both ρ and T show signatures of the spiral arm discontinuity and jump by a factor of three to four, indicating a strong hydrodynamic shock.

6.3 Potential temperature and lessons from it

The potential temperature ϑ is straightforwardly computed from density and temperature profiles (Eq. (10)) and its gradients provide information about the energetic state of the solution. Positive gradient versus gravity indicates a radiative region while a null-gradient indicates convective instability. Strong gradients in potential temperature indicate heavy cooling, and must be followed up by contraction. Thus, we take potential temperature gradients to explain differences in accretion rates based on envelope structures.

As is visible in Fig. 12 using higher opacities, deepening the gravitational potential and under-resolving flattens the potential temperature profile, corresponding to a decreasing Ṁ. This explains why a deeper potential does not increase the mass accretion rate: The temperature in the central region cannot diffuse away efficiently enough, even for κ = 0.01 cm² g⁻¹, and reaches values such that the entropy gradients are decreased overall, which must decrease accretion rates.

Issues with high temperature in the central region of high-mass planets have been reported (Szulágyi et al. 2016). One might worry that those temperatures originate from underresolution of the gravitational potential. Here we show that underresolution of a potential actually lowers the central temperatures. However, the density relative to the temperatures is still lower, which results in an increased entropy. At any given resolution, this unfortunately does not help to decide whether the physically correct solution has been found, but it removes doubts pertaining to overly high temperatures being the cause of non-accretion.

Comparing now a set of low-opacity simulation runs for different smoothing lengths, we can understand, step by step, the role of deepening the gravitational well on the envelope structure (see Fig. 12, straight red and dashed red lines). Both simulations start out with the same envelope mass from the Step B run, but the deeper potential will have a higher temperature and thus a lower density. With Eq. (10) it is clear that log ϑ ~ logT − 0.4logρ, and therefore both those factors contribute to increasing the potential temperature. While for the constant opacities the relative envelope structure in ϑ remains very similar for all r_s, the effect of opacity transitions on entropy is clearly visible in the black and green curves in Fig. 13.

As the horseshoe region penetrates down into 0.5r_H, it is this radius which defines a “feeding region”. There, the local entropy gradient will decide whether a part of the horseshoe mass flows into the inner envelope or not. As we see, the entropy gradients at r = 0.5 r_H are not very different when switching from constant to complex opacities, when leaving rs constant. This could be a possible explanation for their similarities in accretion rates.

Therefore, although not dramatically different in terms of mass accretion rate, those simulation runs with different opacities would treat a tracer dust species differently, having interesting consequences on the transport of dust. We address this in a future work.

Fig. 13 Potential temperatures in simulations with constant (left, red curves), Bell and Lin (right, black curves), and Malygin opacities (right, green curve) for different values of ${\tilde{r}}_{\text{s}}$. We note that both plots share the two legends. The rs = 0.1 have Nc = 100 and correspond to the simulations shown in Fig. 12. A major difference between constant and Bell and Lin opacities is that close to the water-iceline at around $\tilde{r}=0.25r_{\text{H}}$ there are indications of convective instability. The notion that this is related to the vicinity to the smoothing region is weakened, as the trend of inverted ϑ-profile continues as we deepen the potential in a well-resolved manner. The comparison to the more realistic Malygin et al. (2014)-opacities is done only for ${\tilde{r}}_{\text{s}}=0.2$ (compare the green and the straight black line). There, it is evident that the gradients in ϑ are weakened by a multi-step transition in κ.

6.4 Temporal-2D-averages: τ = 1 surfaces, luminosities

The previous paragraph presents a comparison of 1D-(ρ, T, ϑ)-structures for a set of important simulation parameters. We now focus on the comparison of opacities, and show a 2D-averaged version of the full-3D dataof those simulation runs to explain the resulting gas dynamics.

Figure 14 shows the density and temperature in contrast with each other. This is important for understanding the gas dynamics: in regions where density and temperature surfaces are parallel to each other, the 3D simulation can be treated as a 1D problem. From the figures it is evident that the region of spherical symmetry reduces in radial extent as opacity reduces. This is because stronger cooling allows the initially pressure-supported hot gaseous blob to lose pressure, and thus it flattens into a shape that is more strongly supported by rotation. Another important non-spherically symmetric feature at low opacities is the efficient vertical cooling that causes upward bending of the temperature isocontours and downward bending of density contours, which is responsible for the “pancake” structure in the first place. As is evident from the plots, this region starts at the τ = 1 surface, from which photons escape freely.

The run with κ = 0.01 cm² g⁻¹ exhibits an enclosed optically thick surface due to the gap being optically thin. This makes the computation of the energy lost from the envelope particularly simple. We compute the luminosity from this run as L =∑ dAσT⁴, where the surface elements dA are computed from the data cell by cell and T is the temperature of each cell. As an example, according to Eq. (16), the measured luminosity of 2.0 × 10²⁹ erg s⁻¹ for a Saturn-mass planet with characteristic length being the smoothing length of 0.1 r_H corresponds to an accretion rate of 0.01 m_⊕ yr⁻¹, while our measured accretion rate is 0.02 m_⊕ yr⁻¹. Therefore, although there are complications in the 3D structure, one can use the measured luminosities to reasonably understand the differences in accretion rates, and follow the envelope evolution, which we demonstrate below.

6.5 Temporal 2D averages: the baroclinic vortex and gas flows through the Hill-sphere

An interesting phenomenon in our simulations is the appearance of a new circulation feature inside the planetary Hill-sphere, unlike what was observed for gas giants in Lambrechts et al. (2017). It is well known from classical fluid mechanics that any region with non-zero $\frac{1}{{\rho }^2}\nabla \rho \times \nabla P$ poses a local source for vorticity. Expressed differently, any region where iso-ρ and iso-P contours are inclined towards each other will produce vorticity. The same argument applies to iso-ρ and iso-T surfaces. We show those surfaces in Fig. 14. There, it becomes clear that the stronger the cooling, the sharper the transition from spherical to near-isothermal temperature contours and spherical to vertically stratified density contours. This poses a stronger vorticity source for stronger cooling, which destabilises the flow into a vortical feature for κ ≤ 0.1 cm² g⁻¹. The same feature appears for the Bell and Lin opacities already at ɛ ≤ 1.0%. The resulting vortical flows in the outer envelope, can be seen in Fig. 15. The usage of the opacity set provided by Malygin et al. (2014) enlarges the vortex in radial extent by a factor of approximately two compared to the vortex in the simulation with Bell & Lin (1994) opacities, and thus provides a stronger coupling between the disc and envelope entropies.

There are caveats however to the interpretation of this vortical flow. Firstly, Figs. 14 and 15 show time- and cylindrically averaged views of the flow as well as structure variables. In 2D cuts in the radial–azimuthal direction, the correlation between flow and potential temperature are clearer, but for the sake of simplicity we omit those plots. It turns out that the vortical feature is mostly a radial feature, that is it is strongest in the r–z plane, as opposed to the ϕ–z plane. In the r–z plane it also appears centred around r ≈ 0.2 r_H and z ≈ 0.3 r_H, whereas in the cylindrical average it is seen at r ≈ 0.5 r_H and z ≈ 0.3 r_H, so dynamical information extracted from cylindrical averages needs to be taken with caution. The effect on the 3D gas motion is the following: at z> 0, horseshoe flows that enter the planetary Hill-sphere encounter the baroclinic region. The vortical motion resulting from the baroclinicityis superimposed onto those flows, being weaker than the forward momentum in azimuth that they already possess. Thus, the vortex causes modified gravitational swing-by trajectories, and does not pose a closed circulation feature. Secondly, the midplane flows, which transport an important part of the mass flux, remain relatively weakly perturbed by this vortical flow. Therefore, the flows as seen in Fig. 4 remain essentially constant, even in the presence of the vortex. Thirdly, we show in Fig. 16 that the vortex leaves no trace in the overall energetics of the envelope, either in compressional or viscous work. This strengthens the notion that this vortical motion is only a passive tracer of the overall baroclinic structure of the envelope, and is not a strong feedback mechanism into the envelope energetics.

Fig. 14 Time and cylindrically averaged density (as background colour) and temperature (in contours) distributions as functions of r and z and the optically thin–thick transition surface measured from the simulation boundary. Parameters are well resolved (Nc = 100) runs of opacity Set 1 and gravitational smoothing ${\tilde{r}}_{\text{s}}=0.1$. The time-average ran over ten orbits with ten output samples per orbit. Near the planet, density and temperature are both radially distributed. Farther away from the planet, the density stratifies in z in order to connect to the gap. Moving outward radially in the disc midplane, the temperature becomes quasi-isothermal as r moves into the optically thin gap. In the z-direction however, heavy cooling bends the iso-T planes upwards. The luminosity required for Eq. (16) is computed for the κ = 0.01 cm2 g−1 cases in the approximation ∑dAσT4, which is possible as in those cases the planetary gap is optically thin. Luminosities and accretion rates vs. time are shown in Fig. 10.

Fig. 15 Temporally and cylindrically averaged flow stream lines with their colour-coded velocity relative to the planet and potential temperature ϑ in contours for well-resolved runs of opacity Set 1 and gravitational smoothing ${\tilde{r}}_{\text{s}}=0.1$. The plot isslightly zoomed in compared to Fig. 14 in order to better feature the flow structures. Correlation of flows and isentropic surfaces is imperfect due to the averaging process. The vortex feature that appears for κ = 0.1, 0.01 cm2 g−1 is a feature that exists in the x–z plane, not in the x–y plane. This feature stems from families of horseshoe trajectories that become distorted due to the bend in isocontours of ϑ. This again results from the baroclinicity of the envelope, which can be understood from Fig. 14 as the angle between iso-density and iso-temperature contours. Those results show that there is roughly a structure of three sections for accreting material: in the midplane, material is generally flowing into the Hill-sphere, excess material that cannot be cooled is ejected at mid-latitudes. Accretion from the top generally is seen in all cases, but amounts to only 10% of the net accretion.

Fig. 16 Cylindrically averaged compressional (left) and viscous work (right) in the envelope for ${\tilde{r}}_{\text{s}}=0.2$, κ = 0.01 cm2 g−1, which is almost identical to the runs with ${\tilde{r}}_{\text{s}}=0.1$ in Figs. 14 and 15. We chose the run with ${\tilde{r}}_{\text{s}}=0.2$ because the important physics is more visible. Most of the heating in the optically thin part of the envelope (for the thin–thick transition compare to Fig. 14) is generated at the spiral arm shocks (at r ≈ 0.8 rH), but is radiated away efficiently. The optically thick region of the envelope is affected by vertically infalling material, but is dominated by contraction of the inner, spherically symmetric region. White regions on the left indicate decompression, which plays only a minor role for the total luminosity in the envelope. The decompressing region within 0.4 < r∕rH < 0.7 is due to two different effects: At z ≤ 0.2 there is material from the horseshoe flows that decompresses as it falls into the inner envelope. For z > 0.2 the decompression is due to the global meridional flow. The vortical flow does not play a role for the compressional work. Weak decompression is seen above the planet at r ≈ 1.0 rH where strong cooling allows for gas to be in free-fall for a very short distance. The viscous work is volume-integrated a factor of approximately ten times less than the compressional work in the envelope outside of rs, and there is, as for the compressional work, no visible heating from the vortical flow. The flow inside of rs (indicated by the white arrow, also seen in Lambrechts & Lega 2017) is probably an artifact that is responsible for the behaviour seen in Fig. 10 and further discussed in the text. It is presumably the feature that we need to resolve properly with ten grid cells, as found before in Fig. 8.

6.6 The planetary boundary condition and the energy balance from gas dynamics

Another interesting point in our work is that we are able to quantify the influence of the 3D gas flows on the energetics of the envelope. Because the gas compresses and decompresses as it flows through the envelope and is subject to shear which generates viscous heating, one might ask whether this can significantly affect accretion rates. We investigate this by averaging the compressional and the viscous work in the envelope. Figure 16 presentsthe results and we conclude that compressional and shear heating in the outer envelope (at r > r_s) are insignificant for the overall energetics. What is dominant is the compressional work in the central region of the envelope, where density and temperature possess spherically symmetric distributions.

However, the amount of central heating might be overestimated. This is because the gravitational smoothing forces a slow circulation (indicated with a white arrow in Fig. 16) that has non-zero shear. This, on average, heats the envelope viscously as efficiently as it does via compressional work, in a pattern that correlates with the smoothing length and thus seems to be a numerical artifact. This is probably the feature that is necessary to resolve with ten grid cells, explaining our results from Fig. 8.

Figure 10 shows that in an underresolved simulation, the accretion rate continues to quickly decline, and the luminosity runs into a state of constant non-zero luminosity: if this was a star, we would expect a nuclear burning source to supply energy; however, there is no such thing implemented in our simulations. Instead the non-zero luminosity is presumably provided by the viscous work, which is a factor of approximately three stronger than the compressional work in the unresolved case, when integrated over the whole Hill-volume. The situation is more complex for the well-resolved cases: There, the ratio of both work integrals over the smoothing length region is approximately one, but when integrated over the whole Hill-volume this ratio can be as low as one-tenth, giving the false impression that viscous work is unimportant overall for the energy balance in the simulation.

That the viscous work is so important in the smoothing length region, even in the well-resolved simulations, could be an artifact of the non-existent gravity inside r_s: once gas accretes into this region, it still has kinetic energy that needs to be dissipated, as there is no planet to fall on. A region of high shear is generated, until an energy balance is found.

This problem seems to be avoidable by the usage of a lower viscosity. However it seems possible that the flow could also find a solution with even higher shear, amplifying the viscous heating again. A careful future investigation into this matter is needed, or viscous simulations be avoided altogether.

6.7 Rotational properties of the envelope and their evolution

All our envelopes rotate with sub-Keplerian speed within 0.5 r_H, as seen in Fig. 17. It is conceivable that our choice of large disc viscosity α is not realistic for protoplanetary discs, as Saturn probably had a moon-forming circumplanetary disc, and thus was able to exchange pressure for centrifugal support.

The lack of rotational support may be a consequence of accretion flow coming mainly from above with little or no angular momentum. However, there is still a weak radial inflow that may allow the envelope to accrete high-angular-momentum material. As seen in the cylindric flow plots of Fig. 15, there is a net mass flux coming in through the midplane, which can help to speed up the envelope azimuthally, but only if not expelled through the envelope. We now investigate the physical nature of this flow.

One has to distinguish the evolution of the angular momentum and angular velocity separately. The angular velocity can slow down purelybecause angular momentum is conserved as mass is accreted. Constant angular momentum with time in spite of accretion is only rarely observed and reflects the originating conditions for the accreted gas: gas coming from the midplane possesses high angular momentum, while gas from above the planet possesses almost none. As our planetary envelopes usually lose angular momentum over time, we conclude that the mass that is seen flowing in or out through the midplane is in fact only “passing by”, while a significant part of the net accreted mass could conceivably come from the top.

This is uncertain however, as the spiral arm shocks have a contribution in changing the angular momentum: as stream lines are deflected after the shock, the relative angular momentum towards the planet also changes. This issue is not settled in our simulation and a follow-up study will investigate the origin of accreted mass by following the stream lines.

Accreted stream lines can transport dust particles of various sizes. In Bell & Lin (1994) and Malygin et al. (2014) the dust particles correspond to μm -sized dust. If the same opacity were found to be provided by larger dust grains (which have a smaller opacity per unit mass), then the total accreted dust mass would need to be higher. At low opacities, this would require the dust grains to be larger than μm in size, but then it becomes easier for the pressure bump at the gap edge to capture those particles. Whether or not moon-forming particles can cross this gap as a result of simple 3D effects or other physics remains an open question.

Fig. 17 Rotation profiles in the envelopes for κ = 0.01 cm2 g−1; our envelopes barely reach Keplerian rotation. There are radial components in the velocity flow that originate from infalling material. The inner envelope is strongly pressure supported. The pressure support is dependent on opacity, as can be expected from higher internal temperatures in the opaque cases.

7 Discussion

We follow up the main part of our paper with a number of reflections on the results of the various simulation runs.

7.1 Observability

Our estimate of the location of the τ = 1 surface and the temperature there (see Fig. 14) give a first-order approximation of the luminosities and peak wavelength band in which to observe an accreting Saturn-mass planet. For example, the run with κ = 0.01 cm² g⁻¹, ${\tilde{r}}_{\text{s}}=0.1$ has a luminosity of 3 × 10²⁹ erg s⁻¹ at an effective temperature between 100 and 150 K, therefore having black-body emission maxima between 20 and 30 μm.

Those temperatures are however only true for ${\tilde{r}}_{\text{s}}=0.1$, and we have dedicated a significant part of this study to investigating the influence of this numerical parameter in order to estimate simulation behaviour at realistic smoothing lengths of r_s ~ r_Jup, corresponding to ${\tilde{r}}_{\text{s}}~0.002$. Based on the observed scaling behaviour of envelope variables with decreasing smoothing length, the τ = 1-surface moves only a small distance, but temperatures there would rise to about ~200–250 K.

Those moderate temperatures are partially due to the absence of an accretion shock onto our planet, which can increase luminosities significantly (Marleau et al. 2017, 2019). This again is due to the high viscosity chosen. The high viscosity imparts a pressure support upon our envelope through a temperature increase, which is not bound to vanish when going to realistic ${\tilde{r}}_{\text{s}}$, quite the contrary. Even deeper potentials will prevent the free-fall of matter from the disc atmosphere through the strengthened pressure support, a trend we have generally observed.

We stress here however that all this is only true for the Saturn-mass planets. Future work must address the existence of accretion shocks in 3D and observability for a wide range of masses separately.

7.2 Adiabatic 3D effects

As can be seen in Fig. 18, there is the possibility that a detached adiabatic region forms in the envelope of a planet. In this case, we found a potential temperature inversion in the simulation run for nominal Bell and Lin opacities and a potential of ${\tilde{r}}_{\text{s}}=0.1$. This inversion is directly linked to the existence of an ice line when using this opacity model. The ice line is at r = 0.5r_H and thus at 0.5 r_H < z < 1.0 r_H where the Bell and Lin opacities are highest; this is a region that is inefficient in radiating energy. The building temperature gradient is sufficient to induce convection, seen in the erratic behaviour of the stream lines in the time-cylindric average. It is interesting to also compare this to the smooth flows in the constant opacity cases in the region 0.5 r_H < z < 1.0 r_H in Fig. 15.

Now checking with Fig. 9, we see that this simulation does not differ anomalously in terms of accretion rate from the others in the same opacity function family. The simulation run with constant κ = 1.0 cm² g⁻¹ and ${\tilde{r}}_{\text{s}}=0.05$ features a similar adiabatic feature at high z, but a much stronger drop in accretion rate. We speculate that this however might be caused by still underresolving this particular simulation.

A different behaviour of the accretion rates would be expected when resolving the radiative–convective boundary of a planetary envelope sufficiently: once resolved, the accretion rate should become almost constant with deeper and deeper potentials, as convective regions are able to transport high luminosities. Resolving the convective region of the planet would require a much deeper gravitational potential and consequently much higher resolution than we are able to achieve in this study. This is a topic for future research however.

Fig. 18 Potential temperature contours (in black) together with the average flows over ten orbits (left) and opacity map (right) for Bell and Lin opacities with ɛ = 1% and ${\tilde{r}}_{\text{s}}=0.1$. High optical depth in the vertical direction causes radiative transport to be inefficient, triggering a buoyant instability. This is visible in the irregular vortex structures that the time-average produces, compared to the simpler vortices in Fig. 15 and the potential temperature profile (contours) that is inverted in the z-direction within 0.5 rH < z < 1.0rH. The opacity map shows the colocation of potential temperature maxima with the opacity maxima corresponding to a pseudo-iceline according to Fig. 1. The radial location of the iceline is slightly distorted through the cylindrical averaging process. The net effect of this buoyant instability on the accretion rates in this particular run is weak, as only 10% of the net accretion rate falls in vertically, which is then blocked out.

7.3 Decrease of accretion rate with time and the problem of the onset of runaway gas accretion

Our simulation runs experience a slow decrease of accretion rates, as shown in Fig. 6, even when resolving the smoothing length well. Similar effects have been reported (semi-analytical 1D, Pollack et al. 1996; 2D hydrodynamics with Newtonian cooling, Papaloizou & Nelson 2005). In the former case, the effect was mostly balanced by increasing accretion luminosity through planetesimals, resulting in essentially constant accretion rates with time. The latter work did not include accretion luminosities due to planetesimals, and reports a similar effect as quasi-static contraction of the envelope. These latter authors show that gas accretion proceeds from a branch of quasi-static contraction to a branch of runaway accretion as envelope masses increase through ongoing accretion. Those branches are characterised by decreasing mass accretion rate and luminosity on the quasi-static branch (which we would refer to as Kelvin–Helmholtz contraction), and increasing mass accretion rate and luminosity in the runaway branch.

Later, Ayliffe & Bate (2009) presented a long integration run for a 33 m_⊕ core mass where they found a similar but non-identical decrease in mass that appears similar to the one we find in the present study. These latter authors found ṁ ~ t^−0.6, while our data decline as ṁ ~ exp(−c t) with a resolution-dependent constant c. Although they employed self-gravity, their envelope was not massive enough to go into runaway gas accretion. Thus, although the conditions for future runaway gas accretion were given, the envelope of the 33 m_⊕ core remained on the quasi-static contraction branch of gas accretion. They additionally introduced the notion that the decrease of ṁ and L with time that they were observing in their simulations was due to the increased optical depth of the envelope, as more and more material was accreted.

Here we report results that corroborate the findings of Ayliffe & Bate (2009). Furthermore, we find that the evolution of the luminosity can explain the decreasing mass accretion rates; see Fig. 10. For a well-resolved simulation, the luminosity decreases by the exact same amount as the mass accretion rate, and these parameters can be in fact translated into one another via Eq. (16). The post-prediction of accretion rates from the measured luminosity is innaccurate by a factor of ~ 1.5−2, but it does reproduce the decline of accretion rates well, and thus suggests that our planetary envelopes lie on the quasi-static contraction branch of gas accretion. Furthermore, the luminosity has been linked to the structure of the envelope in Sect. 6 and can be explained by the increasing mass in the envelope, which moves the τ = 1 surface to colder temperatures.

However, our planetary mass of 90 m_⊕ would be expected to be in the run-away accretion branch of gas accretion. On this accretion branch, the luminosity increases with increasing envelope mass and the accretion accelerates in time. However, as we actually insert a core of mass 90 m_⊕, and not a real massive envelope, this core is on the slow contraction branch where the gravitational potential is dominated by the core, that is m_env ≪ m_core. We nevertheless expect that the accretion rate that we measure will become independent of the core mass as the envelope mass is increased. Our envelope masses reach at most 20% of the core mass, but even such a moderate envelope mass should have an accretion rate that approaches the run-away branch with m_env > m_core; otherwise it would be very puzzling to find luminosities (10⁵ L_J = 10⁻⁴ L_⊙) and envelope accretion rates (ṁ = 10⁻² m_⊕ for κ = 0.01 cm² g⁻¹ and r_s = 0.1) consistent with runaway gas accretion, while being in quasi-hydrostatic contraction.

Ayliffe & Bate (2012) employed a very massive protoplanetary disc to force the envelope to accrete rapidly enough to capture the quasi-static contraction as well as the onset of run-away gas accretion. Their simulations behaved identically to those of Ayliffe & Bate (2009) at m_env ≪ m_core, and therefore thus demonstrated that the early decreasing mass accretion rates are indeed due to being on the hydrostatic branch of gas accretion. Lambrechts et al. (2017) present a more detailed discussion of the two branches of gas accretion and how 3D simulations on the slow branch can be used to probe the accretion rate on the fast branch as well.

7.4 Scaling of Saturn-mass accretion rates with constant opacity and disc-limited accretion rates

According to standard gas accretion theory, accretion rates should scale like ṁ ~ 1∕κ (Mizuno 1980; Ikoma et al. 2000). However, in 3D radiation hydrodynamics calculations this seems to only be true in the low-mass regime. At higher masses, this κ-dependency gradually weakens until disappearing at planetary masses much higher than that of Saturn (Ayliffe & Bate 2009; D’Angelo & Bodenheimer 2013). The disappearance of the opacity dependence happens when a gas planet reaches the mass accretion limit that can be supplied by the disc.

Our scaling of accretion rates with opacity is ṁ ~ κ^−1∕4. This is a similar scaling to that of Ayliffe & Bate (2009) of ṁ with opacity, as we have a similar setup and also take gap formation into account. Therefore, we have presumably not reached the limit of what the protoplanetary disc can supply. Also, we have a similar setup to that of Lambrechts et al. (2017), who state that the total mass flux through the protoplanetary envelopes can be one to two orders of magnitude higher than the accretion rates. Only when the accretion rates were to come close to those total fluxes, would we speak of disc-limited accretion.

In fact, the accretion rates in our simulations do not seem to be dictated by the surface densities or the volume densities at all. Neither the scaling of minimum gap surface densities nor of the disc surface densities match with the scaling of ṁ. Even though we measure strong accretion rates of ṁ ~ 10⁻² m_⊕ yr⁻¹, the disc is able to supply more than that amount of mass. This is also independent of the branch of gas accretion that we find ourselves on, as accretion rate is not related to the disc supply on either branch of gas accretion at this planetary mass. We also show in this paper that the accretion rates are a consequence of envelope properties alone, which makes sense if the disc can supply enough to satisfy any envelope contraction.

There is also the possibility that the sublinear scaling of accretion rate with inverse opacity is a result of the initial conditions, as a lower opacity implies a more dense midplane and thus the accretion flows can be more massive. This is a hypothesis that will be testable in future work, as we aim to resolve both lower- and higher-mass planets in the state of gas accretion.

7.5 Comparison to isothermal models

In a recent study, Ida et al. (2018) constructed a full model for the gap opening, growth, and migration of gas giants in the isothermal framework. While we cannot compare all of our work with an isothermal model, there are certain parts that we can post-predict from our simulation setup and compare results.

The work by Ida et al. (2018) does not give a full recipe for the gap shape, but instead provides a formula for the minimum gap depth. The gap depth of these latter authors is a strong function of H∕r and therefore the disc temperature. We can use the H∕r-values from our simulation (see Fig. 2) to compare our actual gap depths Σ(t_final)∕Σ₀ with the ones predicted through the gap depth equation from Ida et al. (2018), which is based on 2D isothermal simulations presented in Kanagawa et al. (2018).

For the simulation runs with κ = (1.0, 0.1, 0.01) cm² g⁻¹ we have unperturbed values of H∕r = (0.038, 0.029, 0.027) which would then give Σ∕Σ₀ = (15, 8, 6)% for the Kanagawa et al. (2018)-formula, while we find Σ∕Σ₀ = (38, 28, 28)% in our simulation runs. There is work by Bitsch et al. (2018) which points out that gaps in 3D isothermal runs are less deep than in 2D isothermal runs. These latter authors connect this difference to the meridional flows found by Morbidelli et al. (2014) presumably replenishing the gaps with mass. This effect might play a role in our 3D simulations to set our gap depths as well, even if the flow around Saturn-mass planets has different structure (see Figs. 15 and Fig. 18, also Lambrechts et al. 2017). Midplane flows towards the planet are not observed at all in isothermal simulations and are due to effects of heating and cooling. The impact on gap depth has not yet been studied, but it is clear that one cannot directly compare 3D isothermal runs with 3D radiative runs.

Using the gas accretion rate of Ida et al. (2018) on the runaway accretion branch, we find for our setup an expected ṁ ~ 10⁻¹ m_⊕ yr⁻¹, even considering the effect of gap depth on limiting the accretion rate. Thus, there is a clear bifurcation between our work and isothermal models: on one hand the isothermal gaps are deeper, but on the other, the isothermal accretion rates appear to be be higher.

The work of Lambrechts et al. (2017) indicates that the mass flow through the envelope is of the order of ṁ ~ 10⁻¹ m_⊕ yr⁻¹ for a Saturn-mass planet, while our actual accretion rates are ṁ_p ~ 1–3 × 10⁻² m_⊕ yr⁻¹. This is all the mass that the planet can cool sufficiently while matter flows through the envelope. Thus, we find, as in Lambrechts et al. (2017), that the flow into the Hill-sphere is a poor proxy for the actual gas accretion rate, even for a Saturn-mass planet. Instead, the accretion is mainly regulated by radiative entropy losses, which are not captured in isothermal simulations.

7.6 Viability of a sink-cell approach

In global simulations, it has been a common approach to have a poorly resolved Hill-sphere (since otherwise it is impossible to attain results over several hundreds and thousands of orbits) while at the same time implementing a mass sink. Usually the planet then has a comparatively flat gravitational potential of ${\tilde{r}}_{\text{s}}=0.5-0.8$ and a fraction of the mass inside of ${\tilde{r}}_{\text{s}}=0.5$ is simply removed per time step in order to mimic accretion. A typical recipe (see Klahr & Kley 2006) removes mass from the envelope in a smoother fashion according to a removal function ṁ (r) = f(r)m(r). Those approaches are usually generalised under the umbrella term of “sink cells”. As we do find gas accretion with deeper potentials of ${\tilde{r}}_{\text{s}}=0.05-0.1$, which decouple the horseshoe region from the planetary envelope, we can aim at justifying or refuting those approaches.

To this end, we define the quantity mass accretion inside fractional Hill-radius ${\dot{m}}_{{\tilde{r}}_H}$. This is simply the fraction of the total accretion rate that is accreted by a certain fraction of the Hill-radius. In our simulations this is a straightforward quantity to measure, and we see ${\dot{m}}_{{\tilde{r}}_{\text{H}}=0.2}=0.4\dot{m}$ and ${\dot{m}}_{{\tilde{r}}_{\text{H}}=0.1}=0.2\dot{m}$ for all constant opacities. This means that 40% of the total accreted mass “disappears” into r = 0.2 r_H, and 20% of the total accreted mass “disappears” inside r = 0.1 r_H. Those fractions remain remarkably constant for all constant opacity runs.

In the case of Bell and Lin opacities, we find ${\dot{m}}_{{\tilde{r}}_{\text{H}}=0.2}=0.65\dot{m}$ for the dusty cases of ɛ = 1.0, 0.1% and ${\dot{m}}_{{\tilde{r}}_{\text{H}}=0.2}=0.80\dot{m}$. This means that for a more realistic opacity law one finds more mass being accreted into the gravitational smoothing region, as opposed to outside of it.

This means that if we were to apply a sink cell with radius r = 0.1r_H, a mass equivalent of ${\dot{m}}_{{\tilde{r}}_{\text{H}}=0.1}=0.2\dot{m}$ must be removed from this region in order for the envelope structure to remain consistent with the remaining simulation domain. In an effort to partially resolve the physical processes of the Hill-sphere while using a small sink-cell region, one could produce inconsistent results if one instead applied the full accretion rate to the innermost regions. The full accretion rate is the one consistent with the luminosity according to Eq. (16). We therefore caution against using sink cells with small radii and full accretion rates.

Drawing conclusions from this for the formation of subdiscs and the mass reservoir available for satellite formation of gas giants, unfortunately our potential depths are insufficient to resolve the inner Hill-sphere where moons should form. As 1D simulations (see the work by Fujii et al. 2017, and others) show, circumplanetary discs often only begin at ${\tilde{r}}_{\text{H}}=0.1$. This is the innermost radius where we can draw comparative conclusions for all simulations. Thus we cannot distinguish between mass delivery towards the circumplanetary disc and mass delivery towards the planet.

8 Summary

In this paper, we take a detailed look at the numerical and physical effects that accompany gas accretion onto the envelope of a planetary core. We employed 3D radiative and hydrodynamical simulations to follow the detailed flow patterns from the protoplanetary disc into the Hill-sphere and the gaseous envelope.

The main aim of this paper was to relax as many approximations as possible, particularly by not using a sink-cell approach, while producing a global simulation and sufficiently resolving the planetary Hill-sphere. The Hill-sphere of the planet has an inner, central boundary condition that is a simple, smoothed gravitational potential characterised by the smoothing length r_s. Our main findings from the simulation framework that we have performed are as follows.

• We find that numerical resolution of the gravitational smoothing length plays a crucial role in facilitating the accretion process. Insufficient resolution leads to a quick cessation or even negative accretion, likely due to artificial entropy generation in the envelope at the lowest resolutions.

• Furthermore, simulation runs with ten or more cells per smoothing length resolution capture accretion physics accurately in the sense that the accretion rate is converged.

• For simulation runs with a constant opacity of κ = 0.01 cm² g⁻¹, the gas gaps are optically thin and the planetary envelopes are completely enclosed by their τ = 1 surface integrated from infinity towards the midplane. The luminosity that we estimate is lost from this τ = 1 surface corresponds well to the mass accretion in accordance with Eq. (16). We therefore use it as a tool to understand the simulation behaviour when deepening the gravitational smoothing length, while remaining numerically converged.

• We then observe that deepening the gravitational potential increases the internal energy content of the envelope more than it does the luminosity. As a consequence the accretion rate decreases for a deeper gravitational potential, consistent with Eq. (16). A signature of this increasingly energetic envelope is the flattening out of its potential temperature gradients which leads to less and less systematic inflow into the envelope.

• Extrapolating the smoothing length r_s to a more realistic size of ~ 1 r_J, we obtain an accretion rate of 10⁻³ m_⊕ yr⁻¹ for a low opacity κ = 0.01 cm² g⁻¹.

• Using constant opacities we find that for high opacities, planetary envelopes remain spherically symmetric and connect smoothly to the disc. At low opacities, the planetary envelope possesses a more complex structure determined by radiating away pressure support and replacing it with centrifugal rotation. The envelope then also possesses a completely enclosed optically thick surface, thus disconnecting it from the disc in terms of radiative energy diffusion.

• The scaling of the accretion rate with constant opacity is ṁ ~ κ^−1∕4. Considering the agreement of our results with the results from Ayliffe & Bate (2009), this is likely to be an intermediate result between ṁ ~ κ⁻¹ for low-mass planets and ṁ ~constant for Jupiter-mass planets. An investigation of accretion rates versus planetary mass was beyond the scope of this paper and is deferred to a future study.

• More complex opacity laws, namely the Bell & Lin (1994) and Malygin et al. (2014), differ from the constant opacity case mainly in their shape of the τ = 1 surface and the temperatures at that surface, while the luminosities are relatively unchanged, as are the flows inside the envelopes. Thus, for the accretion rates the influence of using complex opacity laws seems small in our work, but there are three important caveats: (1) in the spirit of constructing self-consistent theories one would rather use the non-constant opacities. (2) In reality line cooling might play an important role in radiating additional thermal energy from the envelope (Gustafsson et al. 2008) (3) Our simulations do not reach temperatures above 2000 K as reaching thehigher temperatures near the core would require prohibitively high resolution. The more realistic non-constant opacity lawswould only show major differences with the constant opacities at such high temperatures. It is nevertheless not clear that the opacity so deep in the envelope has any influence on the total radiative entropy loss and mass accretion rate of the envelope.

• The 3D gas dynamics through the envelope can be understood by combining pressure and temperature into potential temperature, which is a measure of the entropy that links the local conditions in the envelope to a disc temperature at the same entropy. The potential temperature structure appears to be sensitive to the used opacity function; particularly, the entropy reduction is more rapid at low opacity. The hydrodynamical flows in the Hill-sphere result in heating both by shear and compression, but those values are negligible compared to the compressional work of the deep planetary envelope. In the future, 3D simulations reaching down to the planetary core would be necessary to analyse this heating budget deep in the convective envelope.

• The role of viscosity and viscous work in our simulations is complex. In underresolved simulations, viscosity can out-competethe compressional work by a large factor and drive the envelope into a steady state with zero accretion but non-zero luminosity. In resolved simulations, it becomes clear that most of the viscous work amounts to unity compared to the compressional work and is generated by a shear flow inside the gravitational smoothing region. While it is possible to envision that this shear flow is necessary to dissipate the momentum of accreted gas inside the smoothing region, as it has nowhere else to go, it is not certain whether the amount of heat generated by this process is realistic.

• Our envelopes generally accrete less than the total mass flux through the Hill-sphere. Increasing the viscosity would therefore not help to accrete even more mass through feeding the planetary gap. Instead, increasing the viscosity would help to increase the accretion rate via heating the envelope and increasing luminosity, which allows further accretion.

• Planetary envelopes for our simulations generally lack rotational support, as pressure support dominates those envelopes interior of approximately half the Hill-sphere. However, we adopted a relatively high disc viscosity α ~0.01. Future simulations with a lower α would be needed to probe circumplanetary disc formation at viscosities that better reflect protoplanetary disc values.

Our investigations pose a challenge for the simplified view of core accretion, where reaching the critical core implies immediate gas accretion and formation of a gas-giant planet. This critical core mass is usually constructed as being a function of disc and envelope parameters (Mizuno 1980). Our simulations however are set in a regime where any critical core mass is far exceeded, yet we find ourselves in a state of Kelvin–Helmholtz contraction. Therefore, reaching the critical core mass is in itself not a guarantee that the core will accrete substantial amounts of gas. This is also true for low-mass cores, as explored in Lambrechts & Lega (2017), where gas accretion can be stalled by high opacities.

Our work demonstrates that the pre-existing concept of luminosity resulting from the temperature structure of a planetary envelope explains our found accretion rates well. This poses a connection point for future radiation hydrodynamics calculations in one dimension for structure calculations and full 3D radiation hydrodynamics, as ours and other similar setups can compute those luminosities self-consistently. Furthermore, future work exploring a broader physical parameter space than what was possible in our work should be able to profit from this.

Acknowledgements

We thank the anonymous referee for their help in improving the quality of the paper and making a number of constructive suggestions. The authors are indebted to illuminating discussions with Jeffrey Fung, Gabriel-Dominique Marleau and Michiel Lambrechts. The new complex opacity tables have been provided by Nick Malygin himself for which we are very grateful. M.S. and A.J. were supported by a project grant from the Swedish Research Council (grant number 2014-5775) and an ERC Starting Grant from the European Research Council (grant number 278675 – PEBBLE2PLANET). A.J. was further supported by the Knut and Alice Wallenberg Foundation (grant number 2012.0150) and the European Research Council (ERC Consolidator Grant 724687- PLANETESYS). B.B. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. All the simulations presented in this work were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc in Lund University, Sweden.

Appendix A Radiative fluxes in two-temperature flux-limited diffusion (FLD)

Here we present the behaviour of the radiative fluxes in optically thin regions in the flux-limited diffusion approximation with our two-temperature solver. Radiative fluxes should approach |F_rad| = cE_rad in the optically thin regions of the disc atmosphere. We briefly comment on how well this limit is approached in our code.

Figure A.1 shows how the reduced flux in component i f_i,red ≡|F_i|∕cE_rad evolves with height above the midplane. It is important to note that the correct limit f_i,red = 1 does not necessarily indicate the correct solution for the fluxes and energies themselves. This needs to be addressed in the future via the computed moments of a direct radiative transport solution.

Fig. A.1 Reduced vertical radiative flux as a function of the vertical coordinate. We show that in the far upper regions of the disc atmosphere our code reaches the correct flux values. The computed vertically integrated τ = 1 transitions at three different positions in the disc are marked with triangles. The end points of the curves correspond to the 7° angle of our simulations at those disc positions. Compared to the analytical |Frad|∕cErad ~ 0.55 at τ = 1 (Eq. (82.15) in Mihalas & Mihalas 1984) It is evident that our method locates the optically thin–thick transition further inward than the analytical value. This leads to higher temperatures on the τ = 1 surface and thus to an overestimate of L.

We probe three different locations in our disc as examples: the outer and inner rim of our disc simulation boundary, located each ± 6r_H from the planet, and the position directly above the planet. The curves for both rims show how the open boundary conditions for the radiation in z enforce the correct value of f_red. Just above the planet, density gradients are stronger than at the rims because of the gap formation and the accreting envelope. Thus, above the planet f_red achieves its free-streaming limit much earlier than the boundary, indicating correct cooling behaviour without the action of the boundary conditions.

Our computed optically thin–thick transition is somewhat lower than the analytic value would predict, which impacts our luminosity estimate from Sect. 6, which we refer to here as the “simple” estimate. When the actual optically thin–thick transition lies slightly higher, its temperature is lower, thus the actually radiated L_ex is less. This is what we see when we compare the simple with the exact luminosity in Appendix B.

Appendix B Validation of the luminosity estimate

This appendix is dedicated to the discussion of the quality of our approximate computation for the luminosity L in Sect. 6. We show a comparison of the simple estimate with the exact luminosities in Fig. B.1.

Fig. B.1 Different luminosities in comparison. The Kelvin–Helmholtz luminosity that would be required to balance the accretion Ṁ in red,with the accretion rate measured internal to two different shells. Our simple estimate of L from Sect. 6 is marked in green. The exact measurement of L (shown in blue) differs from the simple measurement by 20%. Thus, even with the misplaced τ = 1, the difference in L is relatively small and does not impact our main conclusions.

We attempt a validation of our simple method of estimating the envelope luminosity L_smpl via the approach from Sect. 6, which assumed $L_{\text{smpl}}={\displaystyle \sum \text{d}}A\sigma T_{\tau =1}^4$. We show that L_smpl is correctly estimated down to a systematic uncertainty of 20% and compare it to the luminosity that would be required by Eq. (16), that is the Kelvin–Helmholtz luminosity L_KH.

To calculate the full-3D luminosity L_3D, it is required that we transform the data from the spherical, heliocentric grid onto a new spherical, planetocentric grid, in order to compute the integral L_3D = d_t ∫ dV E_rad = ∮ dAF_{r, planet}. Here, d V is the Hill-volume or a fraction thereof, dA are the intersection surface elements of the planetocentric spherical grid with the heliocentric grid cuboids, and F_{r, planet} is the planetocentric radial flux that must be computed from the heliocentric flux components.

In this method, the orientation of the surface at the optically thin–thick transition is however non-unique in our implementation of the FLD approximation, because we use the same radiative diffusion coefficient in all three spatial directions. Therefore, we can only measure the luminosity on spherical shells of varying radii around the planet. The planetocentric surface elements d A are then computed semi-analytically from vector geometry and the grid data. Finally, the F_{r, planet} are resulting from the scalar product n ⋅F where n is the planetocentric normal vector on a sphere of arbitrary radius around the planet, and F is the heliocentric flux data coming taken from the simulation.

A remaining difference of about (10^29.5 − 10^29.4)∕10^29.5 ≈ 20% between L_smpl and L_3D is evident. This difference seems to be independent of the chosen shell for the exact flux, because nearly all the radiative flux is generated in the very hot region interior to r < r_s. L_smpl can only sensibly be computed on the τ = 1-surface, so no other comparison value exists.

This comparison was made for the simulation run with parameter values k = 0.01 cm² g⁻¹ and ${\tilde{r}}_{\text{s}}=0.1$. We picked this parameter value for the comparison instead of the fiducial k = 0.01 cm² g⁻¹ and ${\tilde{r}}_{\text{s}}=0.2$, as the latter was mainly a means of showcasing different numerical resolutions and finding the convergence criterion, while we intend to work with simulations of ${\tilde{r}}_{\text{s}}=0.1$ in the near future and it is thus reasonable to gauge those correctly.

We conclude that the match of L_smpl overestimates the radiative losses, and the fact that it is L_3D < L_KH indicates that some energy is lost possibly due to advection of low-density gas from the envelope (thus not negatively impacting the mass accretion rates; see also Fig. 4). The close fit of L_smpl with L_KH also occurs for the simulation run with κ = 0.01 cm² g⁻¹ and ${\tilde{r}}_{\text{s}}=0.05$, but not with ${\tilde{r}}_{\text{s}}=0.2$ as seen in Fig. 10. This indicates a larger contribution of radiative transport for the overall cooling versus advection for hotter envelopes.

The 20% discrepancy between L_smpl and L_3D does not influence our main conclusions of the paper, which were concerning the explanation and analysis of why accretion rates decrease with time and why accretion rates decrease with decreasing r_s.

Appendix C Entropy and potential temperature from hydrodynamic data

In this section we aim to derive the potential temperature and show its relation to other hydrodynamic and thermodynamic quantities.

We start with the first law of thermodynamics, expressed in terms of the Helmholtz potential H, for which it is d H = c_pdT. We choose this particular potential because it is simply a shortcut, as we do not want to go through a whole class of statistical physics: $$\text{d}H=T\text{d}S+V\text{d}P,$$(C.1)

and after demanding dS = 0 for an adiabatic process, we get $$c_{\text{p}}\text{d}T=V\text{d}P.$$(C.2)

Now with the ideal gas law PV = NRT and a fixed number of particles N we can change this into $$c_{\text{p}}\text{d}T=\text{NRTd}P/P,$$(C.3)

from which we directly obtain $$\frac{c_{\text{p}}}{\text{NR}}\text{dln}T=\text{dln}P,$$(C.4)

and thus, after integration from the arbitrary state (T, P) to the reference state (ϑ, P₀) under the adiabatic assumption we find a the temperature ϑ: $$\vartheta =T{\frac{P_0}{P}}^{\frac{\text{NR}}{c_{\text{p}}}},$$(C.5)

which is simply the temperature gas blob that we would have if we had moved it adiabatically from (T, P) to the pressure P₀. We can now take this as a new variable and compute it on the whole simulation domain. It is important to note that we know NR = c_p − c_v, and the adiabatic index γ = c_p∕c_v. Therefore, for γ = 1.4 we can rewrite $\frac{NR}{c_{\text{p}}}=\frac{c_{\text{p}}-c_{\text{v}}}{c_{\text{p}}}=1-1/\gamma =2/7\approx 0.29$.

In order to do this in post-processing, we rewrite the preceding definition ϑ = ϑ(P) into ϑ = ϑ(ρ, T) by usage of the ideal gas law. The result with P = ρkT∕μ is, after some algebra, $$\vartheta =T{\left(\frac{\rho }{{\rho }_0}\right)}^{-2/5}.$$(C.6)

The entropy can then be computed via S = k_B ln ϑ.

References

All Tables

All Figures

Fig. 1 Opacity as a function of temperature. Constant opacities take the three different values 1.0, 0.1, and 0.01 cm2 g−1, where in the plot we indicate only 1.0 cm2 g−1. The second opacity set from Bell & Lin (1994) increases with temperature, shows a transition at the water-iceline at ≈ 200 K and a steep decline of the opacity once dust sublimation become relevant at temperatures of 1000 K. The third opacity set from Malygin et al. (2014) shows a much richer structure in terms of opacity transitions, which is due mainly to the various chemical constituents that Semenov et al. (2003) computed. The large difference between the Planck opacity and the Rosseland opacity mainly plays a role in optically thin regions. Midplane temperatures in our simulations at 5.2 AU without a planet range from 30 to 50 K.

In the text

Fig. 2 Profiles of the average aspect ratio in the disc midplane for runs with constant opacity after 400 orbits in Step B. The planet is kept fixed at 5.2 AU. The three different opacities that we use for later accretion rate measurements are used as running parameter. For κ = 0.01cm2 g−1 the gap becomes optically thin, and is then directly heated by the gap walls, which explains the higher gap temperature compared to κ = 0.1 cm2 g−1.

In the text

Fig. 3 Gap depth vs. time and convergence ratio: for three different masses with all other numerical and physical parameters being constant, we compare the resulting gap depth as a function of orbits and the convergence measure (Σ(t) −Σ(t + dt))∕Σ(t). We see that after 400 orbits the gap depth has essentially reached its final value (up to 10–20%). The convergence rates for the gap depths are independent of opacity after a few orbits and thus their relative depth will not change. Different opacities behave identical.

In the text

Fig. 4 Overview and snapshot of typical midplane variables during gas accretion in Step C at the time of ten orbits for our nominal simulation run with constant opacity of κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$, and Nc= 100. The Hill surface of the planet is approximated as a sphere, which is marked with the red circle. Top left: gas mass density in the midplane and stream lines overplotted with their respective momentary velocities. Spiral arm shocks launched near the Hill-sphere of the planet are clearly visible. The horseshoe flows interact with the shocks and are partially deflected. Colours indicate the stream line velocities and highlight flow deceleration at the shocks. Judging from stream lines only it would be impossible to say whether there is ongoing gas accretion, because most of the horseshoe mass merely passes through the planetary envelope, while only a minor fraction is accreted. The planetary gap is approximately ~ 2rH in width (as wide as the plot) and attains its unperturbed disc value at ~ 3rH. For comparison, the simulation domain in Step C extends radially to ±6rH. Top right: temperature structure in the midplane. The cold gas in the gap and the surrounding disc is clearly visible as dark background colour, while the planetary envelope and beyond is heated by compressional heat from accretion and the spiral shocks. A small fraction of internal energy is advected along the horseshoe orbits, visible in separate moderately warm (~ 80 K) arm structures marked with arrows. Bottom left: potential temperature ϑ overlaid with iso-ϑ contours. Entropy generation at the spiral arm shocks is strongest where horseshoe orbits hit the envelope radially. The central region cools most efficiently in the vertical direction, and therefore develops a drop in ϑ, indicating a radiative region. The innermost few pixels are subject to further flattening of the ϑ-profile, as herethe smoothing length provides a region essentially free of gravity. Bottom right: cylindric average of density as background colour and temperature in contours, as well as the τ = 1 surface, denoting the optically thin-thick transition. As this run has a constant opacity, the τ = 1 surface follows essentially an iso-density surface. This surface is bent towards the planet at small r, due to heavy cooling into the disc atmosphere and the consequential contraction. This plot only extends out to 1 rH as otherwise averaging effects smear out the physical properties outside the Hill-sphere.

In the text

Fig. 5 Total envelope mass within the Hill-radius as a function of time, immediately after the potential ramp-up time of two orbits, plotted against the numerical resolution of the envelope. We show the runs with cells per Hill-radius Nc ∈ {25, 35, 50, 65, 100} and smoothing length ${\tilde{r}}_{\text{s}}=0.2$ for a constant opacity of κ = 0.01cm2 g−1. Under-resolving the Hill-radius leads to a gradual decrease in the measured accretion rate, which results in the envelope as a whole achieving a dynamic steady state. Increasing the resolution step by step we find numerical convergence of accretion rates and a less pronounced drop in accretion rates over time. We note that our highest resolution (Nc = 100) drops slightly in accretion rate compared to the case with Nc = 65.

In the text

	Fig. 6 Mass accretion rate as a function of time. We present results of runs with cells per Hill-radius Nc ∈ {25, 35, 50, 65, 100} and smoothing length ${\tilde{r}}_{\text{s}}=0.2$ for an opacity of κ = 0.01 cm2 g−1. These values are two orders of magnitude higher than the necessary accretion rate to form Jupiter within three million years, which is approximately Ṁ = 10−4 m⊕ yr−1.
In the text

	Fig. 7 Central potential temperature for the same simulation parameters as in Figs. 5 and 6, κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$. The steady decrease in potential temperature confirms that the accretion is final and irreversible. Just to be sure, we checked andconfirmed that the entropy in the rest of the envelope decreases as well.
In the text

Fig. 8 Our complete set of all simulation runs for constant κ = 0.01 cm2 g−1 where we varied the smoothing lengths as well as the numerical resolution, plotted as accretion rate vs. smoothing length (left) and as accretion rate vs. cells per smoothing length (right). All accretion rates shown are averages of the firstten accreting orbits after increasing the gravitational potential in Step C; see Figs. 5–7. Left: gas accretion rate inside the Hill-sphere as a function of the smoothing length, for varying resolution numbers. Each dot is a simulation run, with the resolution of the Hill-sphere indicated as cell numbers. It is evident that as resolution increases, accretion rates cluster closer together, indicating numerical convergence. However, more resolution is needed for deeper potential wells. Right: gas accretion rate inside the Hill-sphere plotted against the number of cells per smoothing length. Different colours indicate different smoothing lengths. It is evident that we need to resolve the smoothing length with at least ten cells before convergent results can be achieved (vertical lines). Thus, we define a “well-resolved simulation” as one that hits convergent values of Ṁ at ${\tilde{r}}_{\text{s}}\times N_{\text{c}}>10$. Once converged, the accretion rates differ only slightly when deepening the gravitational potential, indicated with dashed lines.

In the text

Fig. 9 Dependency of the accretion rates on the smoothing length and opacities for simulations with the highest available resolution as depicted in Fig. 8. Left: averaged accretion rates from 5 to 15 orbits, for ${\tilde{r}}_{\text{s}}=0.5$ down to ${\tilde{r}}_{\text{s}}=0.05$ and our three different opacity values. Simulation runs with ${\tilde{r}}_{\text{s}}=0.5$ have such a shallow potential that the horseshoe flows, which penetrate down to approximately ${\tilde{r}}_{\text{s}}=0.5$ also interfere with the envelope cooling process. Right: accretion rates with Bell and Lin opacities, with 1 (=solar), 0.1, and 0.01% dust opacities, as well as a run with Malygin opacities and ɛ=0.01%. Simulation runs with ${\tilde{r}}_{\text{s}}=0.2,0.1$ have been shown in Fig. 8 to be sufficiently resolved. We therefore interpret the seen downwards trend of Ṁ with ${\tilde{r}}_{\text{s}}$ as real (red dashed line), even if runs with ${\tilde{r}}_{\text{s}}=0.05$ are only at the verge of fulfilling our convergence criterion.

In the text

Fig. 10 Evolution of measured luminosity vs. measured accretion rate for our κ = 0.01 cm2 g−1, ${\tilde{r}}_{\text{s}}=0.2$ simulation run. Sufficiently resolved simulations show behaviour such that the measured luminosity and the measured mass accretion rate correlate well. Unresolved simulations however run into a state of zero accretion rate, but finite luminosity, which remains higher than that of the final luminosity in the resolved simulation. This hints at numerical effects or other physics providing an additional energy source in unresolved simulations.

In the text

Fig. 11 Density and temperature structures for constant opacity runs in radial and vertical slices after ten orbits: data are shown for runs with ${\tilde{r}}_{\text{s}}=0.1$ and the nominal run with ${\tilde{r}}_{\text{s}}=0.2$ and Nc= 100 for comparison. A deeper potential raises the whole envelope temperature as long as no convective region is formed. Optically thick regions are hotter for higher opacity (see inset), while the optically thin regions approach isothermal structure and their temperature trend reverses according to the position of the τ = 1 surface (indicated with triangles). The density structure follows simple expectations: a hotter envelope has lower density in optically thick regions.

In the text

Fig. 12 Potential temperatures for opacity set 1 with Nc = 100. Deeper potential and lower opacity lead to a flatter potential temperature profile. Radially, the simulation runs with ${\tilde{r}}_{\text{s}}=0.1$ have flat ϑ-profiles; this does not denote the onset of convection, but rather the action of the horseshoe flows that can freely enter and exit the gravitational well. In z, a steeper profile corresponds to stronger accretion flows from this direction. The much higher values of ϑ result from the low density in the disc atmosphere, as we have rH≅Hdisc throughout our work.

In the text

Fig. 13 Potential temperatures in simulations with constant (left, red curves), Bell and Lin (right, black curves), and Malygin opacities (right, green curve) for different values of ${\tilde{r}}_{\text{s}}$. We note that both plots share the two legends. The rs = 0.1 have Nc = 100 and correspond to the simulations shown in Fig. 12. A major difference between constant and Bell and Lin opacities is that close to the water-iceline at around $\tilde{r}=0.25r_{\text{H}}$ there are indications of convective instability. The notion that this is related to the vicinity to the smoothing region is weakened, as the trend of inverted ϑ-profile continues as we deepen the potential in a well-resolved manner. The comparison to the more realistic Malygin et al. (2014)-opacities is done only for ${\tilde{r}}_{\text{s}}=0.2$ (compare the green and the straight black line). There, it is evident that the gradients in ϑ are weakened by a multi-step transition in κ.

In the text

Fig. 14 Time and cylindrically averaged density (as background colour) and temperature (in contours) distributions as functions of r and z and the optically thin–thick transition surface measured from the simulation boundary. Parameters are well resolved (Nc = 100) runs of opacity Set 1 and gravitational smoothing ${\tilde{r}}_{\text{s}}=0.1$. The time-average ran over ten orbits with ten output samples per orbit. Near the planet, density and temperature are both radially distributed. Farther away from the planet, the density stratifies in z in order to connect to the gap. Moving outward radially in the disc midplane, the temperature becomes quasi-isothermal as r moves into the optically thin gap. In the z-direction however, heavy cooling bends the iso-T planes upwards. The luminosity required for Eq. (16) is computed for the κ = 0.01 cm2 g−1 cases in the approximation ∑dAσT4, which is possible as in those cases the planetary gap is optically thin. Luminosities and accretion rates vs. time are shown in Fig. 10.

In the text

Fig. 15 Temporally and cylindrically averaged flow stream lines with their colour-coded velocity relative to the planet and potential temperature ϑ in contours for well-resolved runs of opacity Set 1 and gravitational smoothing ${\tilde{r}}_{\text{s}}=0.1$. The plot isslightly zoomed in compared to Fig. 14 in order to better feature the flow structures. Correlation of flows and isentropic surfaces is imperfect due to the averaging process. The vortex feature that appears for κ = 0.1, 0.01 cm2 g−1 is a feature that exists in the x–z plane, not in the x–y plane. This feature stems from families of horseshoe trajectories that become distorted due to the bend in isocontours of ϑ. This again results from the baroclinicity of the envelope, which can be understood from Fig. 14 as the angle between iso-density and iso-temperature contours. Those results show that there is roughly a structure of three sections for accreting material: in the midplane, material is generally flowing into the Hill-sphere, excess material that cannot be cooled is ejected at mid-latitudes. Accretion from the top generally is seen in all cases, but amounts to only 10% of the net accretion.

In the text

Fig. 16 Cylindrically averaged compressional (left) and viscous work (right) in the envelope for ${\tilde{r}}_{\text{s}}=0.2$, κ = 0.01 cm2 g−1, which is almost identical to the runs with ${\tilde{r}}_{\text{s}}=0.1$ in Figs. 14 and 15. We chose the run with ${\tilde{r}}_{\text{s}}=0.2$ because the important physics is more visible. Most of the heating in the optically thin part of the envelope (for the thin–thick transition compare to Fig. 14) is generated at the spiral arm shocks (at r ≈ 0.8 rH), but is radiated away efficiently. The optically thick region of the envelope is affected by vertically infalling material, but is dominated by contraction of the inner, spherically symmetric region. White regions on the left indicate decompression, which plays only a minor role for the total luminosity in the envelope. The decompressing region within 0.4 < r∕rH < 0.7 is due to two different effects: At z ≤ 0.2 there is material from the horseshoe flows that decompresses as it falls into the inner envelope. For z > 0.2 the decompression is due to the global meridional flow. The vortical flow does not play a role for the compressional work. Weak decompression is seen above the planet at r ≈ 1.0 rH where strong cooling allows for gas to be in free-fall for a very short distance. The viscous work is volume-integrated a factor of approximately ten times less than the compressional work in the envelope outside of rs, and there is, as for the compressional work, no visible heating from the vortical flow. The flow inside of rs (indicated by the white arrow, also seen in Lambrechts & Lega 2017) is probably an artifact that is responsible for the behaviour seen in Fig. 10 and further discussed in the text. It is presumably the feature that we need to resolve properly with ten grid cells, as found before in Fig. 8.

In the text

	Fig. 17 Rotation profiles in the envelopes for κ = 0.01 cm2 g−1; our envelopes barely reach Keplerian rotation. There are radial components in the velocity flow that originate from infalling material. The inner envelope is strongly pressure supported. The pressure support is dependent on opacity, as can be expected from higher internal temperatures in the opaque cases.
In the text

Fig. 18 Potential temperature contours (in black) together with the average flows over ten orbits (left) and opacity map (right) for Bell and Lin opacities with ɛ = 1% and ${\tilde{r}}_{\text{s}}=0.1$. High optical depth in the vertical direction causes radiative transport to be inefficient, triggering a buoyant instability. This is visible in the irregular vortex structures that the time-average produces, compared to the simpler vortices in Fig. 15 and the potential temperature profile (contours) that is inverted in the z-direction within 0.5 rH < z < 1.0rH. The opacity map shows the colocation of potential temperature maxima with the opacity maxima corresponding to a pseudo-iceline according to Fig. 1. The radial location of the iceline is slightly distorted through the cylindrical averaging process. The net effect of this buoyant instability on the accretion rates in this particular run is weak, as only 10% of the net accretion rate falls in vertically, which is then blocked out.

In the text

Fig. A.1 Reduced vertical radiative flux as a function of the vertical coordinate. We show that in the far upper regions of the disc atmosphere our code reaches the correct flux values. The computed vertically integrated τ = 1 transitions at three different positions in the disc are marked with triangles. The end points of the curves correspond to the 7° angle of our simulations at those disc positions. Compared to the analytical |Frad|∕cErad ~ 0.55 at τ = 1 (Eq. (82.15) in Mihalas & Mihalas 1984) It is evident that our method locates the optically thin–thick transition further inward than the analytical value. This leads to higher temperatures on the τ = 1 surface and thus to an overestimate of L.

In the text

Fig. B.1 Different luminosities in comparison. The Kelvin–Helmholtz luminosity that would be required to balance the accretion Ṁ in red,with the accretion rate measured internal to two different shells. Our simple estimate of L from Sect. 6 is marked in green. The exact measurement of L (shown in blue) differs from the simple measurement by 20%. Thus, even with the misplaced τ = 1, the difference in L is relatively small and does not impact our main conclusions.

In the text