{'fraction_non_alphanumeric': 0.03172836403241056, 'fraction_numerical': 0.01587016624178665, 'mean_word_length': 4.463545412934022, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reduced gas accretion on super-Earths and ice giants

A large fraction of giant planets have gaseous envelopes that are limited to about 10 % of their total mass budget. Such planets are present in the Solar System (Uranus, Neptune) and are frequently observed in short periods around other stars (the so-called Super-Earths). In contrast to these observations, theoretical calculations based on the evolution of hydrostatic envelopes argue that such low mass envelopes cannot be maintained around cores exceeding five Earth masses. Instead, under nominal disc conditions, these planets would acquire massive envelopes through runaway gas accretion within the lifetime of the protoplanetary disc. In this work, we show that planetary envelopes are not in hydrostatic balance, which slows down envelope growth. A series of 3-dimensional, global, radiative hydrodynamical simulations reveal a steady state gas flow, which enters through the poles and exits in the disc midplane. Gas is pushed through the outer envelope in about 10 orbital timescales. In regions of the disc that are not significantly dust-depleted, envelope accretion onto cores of about five Earth masses can get stalled as the gas flow enters the deep interior. Accreted solids sublimate deep in the convective interior, but small opacity-providing grains are trapped in the flow and do not settle, which further prevents rapid envelope accretion. The transition to runaway gas accretion can however be reached when cores grow larger than typical Super-Earths, beyond 15 Earth masses, and preferably when disc opacities are below kappa=1 cm^2/g. These findings offer an explanation for the typical low-mass envelopes around the cores of Super-Earths.


Introduction
In the core accretion scenario, giant planets form by attracting a gaseous envelope onto a previously formed solid core (Pollack et al. 1996). These large cores of a few Earth masses (M E ) form from planetary embryos that are able to accrete larger than kmsized planetesimals (Rafikov 2004) and sweep up cm-sized pebbles with high efficiency (Lambrechts & Johansen 2012). As the core grows large, the envelope is accreted from the gaseous component of the protoplanetary disc, which only remains present for a few million years around young stars (Ribas et al. 2015). Core formation must therefore proceed at a rate of several Earth masses per 10 6 year.
The accretion of the gaseous atmosphere itself can be divided in two separate phases. Early on, the envelope slowly gains mass as it cools down, a process which can be further delayed by continued heat deposition from solid accretion. However, if the cooling process can proceed efficiently and for a sufficient amount of time, the envelope can grow as massive as the core. This process triggers a second phase of rapid gas accretion. The onset of a self-gravitating gas envelope facilitates continued envelope cooling and mass growth. Planets that enter this regime become gas giants with a high-mass envelope.
In this work we focus on planets in the process of attracting a gaseous envelope, before reaching the point of runaway gas accretion. This growth phase has been studied extensively with 1D numerical models that assume a spherical symmetric atmosphere and evolve the envelope assuming hydrostatic equilibrium at all times, a method with origin in stellar evolution calculations (Mizuno 1980;Ikoma et al. 2000;Piso & Youdin 2014). Such models show that the timescale to reach runaway gas accretion is very sensitive to the solid accretion rate, the core mass, the dust opacity in the outer envelope, and the amount of dissolved solids (Stevenson 1982). Recent work has argued that the dust grains do not contribute much to the opacity, because they coagulate and settle rapidly (Ormel 2014). This reduces the timescale for which a completed 10-Earth mass core starts runaway gas accretion to ∼ 1 Myr, which is shorter than the disc lifetime (Mordasini 2014). Similarly, the addition of sublimated solids raises the mean molecular weight of the envelope and even further reduces the time to runaway gas accretion (Hori & Ikoma 2011;Venturini et al. 2015Venturini et al. , 2016. Together these findings hint that all cores above a few Earth masses should have massive gas envelopes. These theoretical findings are however in conflict with the existence of the ice giants in the Solar System: Uranus and Neptune. Both planets are larger than 10 M E , but never accreted a massive gas envelope. It also does not agree with the characterized population of short-period exoplanets. Super-Earths, planets with sizes between 1-4 Earth radii (R E ), are found around about half of the sun-like stars (Winn & Fabrycky 2015). Also these planets, with masses up to about 25 M E , typically have hydrogen/helium envelopes that do not exceed 10% of the total mass of the planet (Hadden & Lithwick 2016).
New hydrodynamical simulations challenge the standard assumptions used in envelope growth models. Previous evolution studies (Pollack et al. 1996;Ikoma et al. 2000;Piso & Youdin 2014) assumed envelopes in hydrostatic balance that are closed systems. Interaction with the protoplanetary disc would only occur through an imposed outer boundary condition that matches the unperturbed protoplanetary disc. Instead, hydrodynamical simulations performed in three dimensions show that gas con-Article number, page 1 of 20 arXiv:1708.00767v1 [astro-ph.EP] 2 Aug 2017 A&A proofs: manuscript no. paper tinuously flows from the poles through the planetary envelope (Ayliffe & Bate 2012;Tanigawa et al. 2012;D'Angelo & Bodenheimer 2013). Moreover, recent simulations that were performed under the simplifying assumption of a constant temperature envelope and disc, quantified that gas is continuously advected through the envelope on orbital timescales, implying that almost no gas accretion would occur (Ormel et al. 2015;Fung et al. 2015). However, the simplifying assumption of an isothermal envelope corresponds to the limit case of a completely cooleddown planet, which cannot be obtained for planets heated by the process of formation. Therefore, these new results are difficult to interpret. Furthermore, hydrodynamical studies in 3D that have a more complete thermodynamical treatment of the envelope, claim results well in line with envelope growth rates from standard hydrostatic models (D'Angelo & Bodenheimer 2013).
Here, we aim to find find the quasi-static structure of a lowmass envelope around a massive core (5-15 M E ). We do not make an isothermal approximation, but instead use a radiative treatment similar to D' Angelo & Bodenheimer (2013). Additionally, we model the heat deposited by solid accretion. In order to correctly describe disc-planet interaction, we model a full annulus of the protoplanetary disc around the planet. Hydrodynamical simulations covering the complete evolution of envelope growth in time are numerically unfeasible. Therefore, we make temporal snapshots at different growth stages that show how envelopes are structured and grow.
The paper is structured as follows. A technical description of the methods can be found in Section 2. We first study the case of a planet with a luminosity that is dominated by the heat released by the accretion of solids (Sec. 3). This case is investigated in detail to establish the dynamical and energetic structure of planetary envelopes in their formation stage. Then, we proceed by studying planetary envelopes after solid accretion has come to a halt (Sec. 4). Our results depend on the opacity in the envelope and the mass of the solid core of the planet. High dust opacities tend to delay envelope growth (Sec. 5). Large cores facilitate the transition to runaway gas accretion (Sec. 6). In Section 7, we summarize the implications of our work on the growth of the envelope and the core. We place our results in the context of exoplanet observations that reveal a high occurence rate of Super-Earths, but argue for a low frequency of more massive gas giants (Sec. 8). Finally, a summary is presented (Sec. 9).

General description of the model
We model a full annulus encompassing the orbit of the planet, as illustrated in Fig. 1. This allows us to carefully analyze the exchange of gas between the envelope and the disc. The simulations are performed in 3D, as gas motion perpendicular to the midplane regulates the gas dynamics around the envelope. We use a non-uniform grid to be able to trace the flow of gas on small scales well inside the Hill sphere of the planet, following the approach by (Fung et al. 2015). Additionally we include the heat deposited by the accretion of solids, and model energy transfer by radiation, similar to D' Angelo & Bodenheimer (2013). In order to do so, we make use of the FARGOCA code (Lega et al. 2014). This radiation-hydrodynamics code is designed for global 3D studies of planet-disc interaction. It is based on the finite difference code FARGO (Masset 2000). Below, we report in more details the methods used to model the planetary envelope. Global view on the vertically-integrated density (Σ gas in g/cm 2 ) of the simulated annulus, taken from run RADL27k0.01. We can identify the spiral arms originating from the planet located at coordinates (1, 0) revealing a strong interaction between planet and disc.

Gas dynamics
We simulate a non-selfgravitating viscous gas on a fixed grid, corotating with the planet . The standard continuity and momentum equations, are solved using spherical coordinates (see Lega et al. (2014) for more details on the implementation in FARGOCA). Here, D/Dt is the Lagrangian time derivative, ρ is the density, P is the gas pressure, u is the velocity and T is the viscous stress tensor. The gravitational potential Φ is the sum of the contribution of the stellar potential, the potential of the planet and indirect terms that take into account the acceleration of the primary due to the gravity of the planet and the disc. The planetary potential is given by Φ p = −GM c /r outside a radius r sm form the planet. Here G is the gravitational constant and M c is the core mass. Within the smoothing length r sm , the central singularity is avoided by a cubic interpolation to Φ p = −2GM c /r sm (Klahr & Kley 2006). In appendix A, we describe this prescription for the planetary potential in more detail and explore its influence on our results.

Energy equation
We separately evolve the thermal energy density (E gas ) and radiative energy density (E rad ). The evolution equations are given Table 1. Main parameters for all the simulation that are listed in column 1. The following columns give the gas surface density Σ p at the position of the planet in units of M * /a 2 p , viscosity ν in units of a 2 p Ω p , the aspect ratio H/r at the position of the planet, the core mass M c in Earth mass, the accretion luminosity L acc in erg/s, the opacity κ in cm 2 /g and the smoothing length r smooth as fraction of the Hill radius. The second-to-last column gives the number of orbits at the position of the planet, that were performed at the corresponding resolution level listed in the last column (see Sec. 2.4 for a detailed description of the non-uniform grid that was used).  (Kley 1989;Commerçon et al. 2011;Bitsch et al. 2013): Here, F is the radiative flux. The fluid-radiation interaction term is described by I = ρκ(4σ SB T 4 −cE rad ), where σ SB is the Stefan-Boltzmann constant, κ is the Rosseland mean opacity, T is the temperature and c is the speed of light. In the optical thick limit, this term reduces to energy diffusion of the form D∇ 2 E gas . The diffusion coefficient, corresponds to the product of the velocity and the mean free path of the photons and the last factor in brackets weights the expression with the ratio between the specific heat of the photons with respect to the gas (c V ). The term P∇ · u in Eq. 4 represents compressional heating and Q + viscous heating (for more detail see Mihalas & Mihalas 1984). The final term, l, is a source term which represents the deposited heat per unit time by solid accretion. The implementation of this term is discussed in more details in Appendix B. The radiative flux is calculated in the framework of the fluxlimited diffusion approximation, where f λ is a flux-limiter (Levermore & Pomraning 1981). The implementation is described in more detail in Lega et al. (2014).

Opacity
Both gas and dust act as a source of opacity which absorbs and deflects radiation. Below temperatures of around 1700 K dust particles are the dominant opacity source (Bell & Lin 1994). In this work, we purposely choose a simple constant opacity prescription for our envelope and gas disc. We then explore this free parameter in Section 5.

Equation of state
Similar to our choice of a constant opacity, we have a employed a simple equation of state. We use an ideal gas, where the pressure is given by Here R is the gas constant. We used a fixed value for the mean molecular weight (µ = 2.35), corresponding to a solar H/He mixture. More realistic prescriptions would need to take into account the pollution of the envelope by heavy elements and the changes to the mean molecular weight and adiabatic index. Additionally, variables like the specific heat at constant volume c V are temperature sensitive and depend on the unknown ratio of ortho-to parahydrogen (D'Angelo & Bodenheimer 2013). These topics can be investigated in future work.

Coordinate system
We use spherical coordinates (r, θ, ϕ) where r is the radial distance from the origin, θ the polar angle measured from the z-axis and ϕ is the azimuthal coordinate starting from the x-axis. The midplane of the disc is at the equator θ = π 2 and the origin of the coordinates is centered on the star. We work in a coordinate system which rotates with angular velocity, The code operates in code units with the stellar mass M * , the semi-major axis of the planet a p and the Keplerian frequency set to unity.

Simulation domain and boundary conditions
Radially, the disc extends from a min /a p = 0.7 to a max /a p = 1.3. In the vertical direction the disc extends from the midplane (θ 90 • to 6 • above the midplane, i.e. θ 84 • ). We do not study inclined planets orbits therefore we do not need to extend the disc below the midplane.
We use periodic boundary conditions in the azimuthal direction, as we aim to simulate the full 2π-annulus. We use a reflective boundary on the midplane. This is appropriate for noninclined planets, where mirror symmetry can safely be assumed between the upper and lower hemisphere of the disc. The upper boundary condition is also reflective, to prevent mass outflow. In the radial direction we use evanescent boundaries with a wave suppression zone covering 10% of the inner and outer radius, to prevent the reflection of density waves (de Val-Borro et al. 2006).

Disc model
Before inserting the planet, we bring the disc in radiative equilibrium (Lega et al. 2015). This equilibrium is obtained for a 2D (r, z) axisymmetric disc with resolution (N r , N θ , N ϕ ) = (200, 52, 2) in about 100 Ω −1 p . Viscous stresses are the only heating source, as we ignore stellar heating. This is a good description for the inner disc, where viscous heating dominates. We have used a constant viscosity ν = 10 −5 a 2 p Ω p (code units). This low viscosity will be necessary to properly capture the 3-dimensional flow. We noticed higher values of viscosity (of the order ν = 10 −4 a 2 p Ω p ) effectively force 2-dimensional flow, in agreement with Fung et al. (2015). We also verified that extremely low values of ν = 10 −6 a 2 p Ω p do not alter our results significantly, but at times may produce unsteady fluctuations.
The gas surface density follows a power law profile, Σ gas = Σ p (a/a p ) −1/2 g/cm 2 , which gives an equilibium disc without radial gas motion. At the location of the planet, we record a surface density of Σ p = 1.3 × 10 −3 M * /a 2 p , corresponding to Σ p = 400 g/cm 2 at a distance of 5.2 AU. The aspect ratio is the result of the equilibrium between viscous heating and radiative cooling. Depending on the choice of the opacity, we find a gas scale aspect ratio of H/a p = 0.03 for κ = 0.01 cm 2 /g or H/a p = 0.05 for κ = 1 cm 2 /g at the position of the planet. We note however that our results on the small scale of planetary envelopes only weakly depend on the exact choice of the large scale disc structure.

Introduction of the planet
After the procedure to generate an equilibrium disc, we refine the grid azimuthally and work in three dimensions. Planets of 5 or 15 Earth masses are then embedded in the equilibrium disc and held on fixed circular orbits in the midplane, with code units r p = 1, ϕ p = 0, θ p = π/2. Our mesh is chosen so that the planet is at the intersection between cell's surfaces in the 3 direction, i.e. at the corner of 8 grid cells.
Since hydrodynamical 3D calculations of fully radiative discs are expensive in computational time, we split our simulations in 2 phases. In the first phase, we employ a uniform grid so that we can use the FARGO algorithm (Fast Advection in Rotating Gaseous Objects, Masset 2000) and benefit from large time steps. We recall that in a Keplerian disc the traditional Courant condition provides very small time steps due to the fast orbital motion at the inner boundary of the disc. In the FARGO algorithm, the time step is instead limited by the perturbed azimuthal velocity with respect to the Keplerian one. We run each simulation with the uniform grid until the disc relaxes to the presence of the planet. This occurs after approximately 30 Ω −1 p . Fig. 1 shows the gas surface density when the stationary state is reached. For this uniform grid, our nominal resolution corresponds to (N r , N θ , N ϕ ) = (200, 52, 2048) grid cells. This choice provides a resolution of 0.003 a p , corresponding to 10 gridcells in the Hill sphere of a 5 M E planet.

Increased resolution around the planet
In the second phase, we use a non-uniform grid with increasing resolution towards the planet to more accurately model the planetary envelope. The grid layout is based on the formulation by Fung et al. (2015), which gives small nearly uniform grid cells in the Hill sphere around the planet and large grid cell sizes further out.
Our medium resolution corresponds to a width of the inner grid cells of 0.00135 a p , which resolves the Hill sphere of a 5 M E planet by about 24 grid cells. In total the grid has (N r , N θ , N ϕ ) = (200, 52, 1512) number of cells in the radial, polar and azimuthal direction. At this resolution we have about 1 orbit in 8 hours of computation with 60 cores, using hybrid MPI+OpenMP parallelization. We typically run the code for at least 8 orbits to relax the system to equilibrium.
In high resolution runs, we use inner grid cells of width 4.27 × 10 −4 a p , corresponding to about 84 grid cells across a Hill sphere for a 5 M E planet. This finer grid has (N r , N θ , N ϕ ) = (300, 76, 2268) number of cells in the radial, polar and azimuthal direction. At this resolution we have about 1 orbit in 50 hours of computation with 60 cores. In this final step, we simulate 4 planetary orbits, which we found to be sufficient to capture the interior structure of the planet.
Finally, we note that the FARGO advection algorithm cannot be used on the nonuniform grid; instead the time step is set by the usual Courant condition (Stone & Norman 1992).

Simulations
Our main numerical experiments are described in detail in Sections 3, 4, 5 and 6. Table 1 gives an overview of the performed simulations and the parameter space covered. Additionally, we performed several performance tests of the code. In Appenix A, we present several convergence tests. We investigate the dependency on the resolution in the code, with a focus on properly resolving the planetary potential. Additionally, long-term integrations are performed to demonstrate steady-state stability.

Overview
We start with a presentation of our results on the envelope around a planetary core in the process of accreting solids (run RADL27k0.01). The structure of the envelope is discussed in detail. This allows us to compare, to this reference case, results from following sections where we explore the influence of key parameters. In particular, we discuss the role of interior heating by solid accretion in Section 4 (run RADL0k0.01). We then proceed to investigate the role of the opacity in Section 5, when accreting solids (run RADL27k1), or not (run RADL0k1). We leave exploring the role of the core mass itself to Section 6 (runs RADL27k1M15 and RADL0k1M15).

Luminosity from accreting solids
The structure of the envelope is set by the balance between the central gravitational attraction of the planetary core and outward pressure from the heated interior. We now consider a planet where the interior heat generation is dominated by energy deposition from solids sinking towards the core. For such a planet, in the process of forming its 5 M E -core, the accretion luminosity is Article number, page 4 of 20 M. Lambrechts and E. Lega: Low-mass gas envelopes whereṀ s is the solid accretion rate and r c and ρ c the core radius and density. Here, we have assumed most heat is ultimately deposited deep down in the interior of the envelope . The magnitude of the solid accretion rateṀ s is poorly constrained. It depends on the unknown amount of solids available in the protoplanetary disc, but must be of the order of 10 M E /Myr in order to grow the core to completion before disc dissipation.
In Appendix E, we motivate the solid accretion rate in more detail. . Relevant length scales for a low-mass envelope, expressed with respect to the disc scale height, as function of planetary mass. The yellow line represents the Hill radius (r H ), the blue line the Bondi radius (r B ), and the black curve the radius above which radiative energy transfer is possible (r rcb ), assuming κ = 0.01 cm 2 /g and H/a p = 0.03.

Envelope structure
We first note that gas is perturbed up to a scale height, away of the core. Here c s is assumed to be the unperturbed sound speed. Therefore, gas feels the presence of the core even outside the gravitational reach of the core, given by the Hill radius, where M is the stellar mass. Outside this radius, the stellar tidal field starts dominating over the gravity of the core. Since we consider low mass envelopes, we make the approximation that the core dominates the potential (M p ≈ M c ). The Bondi radius, corresponds to the radius where the pressure is perturbed by the gravitational force from the planet, assuming no temperature perturbation. This latter assumption makes this scale, often used in envelope studies, not practical in the context of our radiative hydrodynamical calculations. We therefore prefer scaling our results with respect to the Hill sphere. However, for completeness we have indicated the Bondi radius with a square symbol in Fig. 2.
Interior to these outer length scales, spherical hydrostatic models predict the existence of a radiative-convective boundary (Appendix D). This corresponds to the location in the envelope where heat transfer goes from being convection-dominated in the interior to radiation-dominated in the outer envelope. The depth of this radiative shell, r rcb , is given by where T mid , P mid are the midplane temperature and pressure and ∇ ad the adiabatic temperature gradient. We discuss this expression in more detail in Appendix D.
The radius of the core itself is located deep inside the envelope, at about 10 −3 Hill radii, and cannot be resolved.
In this paper, we consider low-mass planets embedded in the gas disc. Therefore the above length scales are typically ordered as r rcb r B r H H, as can be seen in Fig. 2. The ordering does depend on the planetary mass and solid accretion rate. Fig. 3 illustrates these various length scales as function of planetary mass. We indicate the depth of the radiative zone with the black curve, assuming a luminosity provided by the core accreting pebbles (see Fig. E.1, Appendix E). Cores with masses above about 20 M E are typically not longer embedded in discs with aspect ratios below H/a p = 0.05. Then the Bondi radius exceeds the Hill radius, which has become larger than the local gas scale height of the disc.
Taken together, the panels of Fig. 2 indicate a good qualitative agreement with the expected thermodynamical structure of a planet, based on 1 dimensional analytical calculations. Overall, we note little temperature increase in the part of the envelope we model (Panel A of Fig. 2). The outer region of the envelope remains nearly isothermal, while interior to approximately 0.2 r H the temperature slope steepens. This transition occurs around r rcb , which we calculated for the given accretion luminosity. This gives a first hint that heat transport in the interior of the planet does not occur radiatively. This result is thus in line with standard 1D models. Note that the turnover of the slope at the smallest scales towards 0.01 r H is largely an artefact of the necessary smoothing of the potential.
Panel B of Fig. 2 shows the azimuthally-averaged density structure in the envelope (black curve), as measured in the midplane. We recover the characteristic exponential increase in the nearly isothermal outer layer. This can be seen by comparison to the gray dotted line in panel B of Fig. 2, which gives the density profile for a purely isothermal envelope (an isotherm). Interior to approximately 0.2 r H , corresponding to the location of r rcb , we observe a turnover to a power law-like slope. This turnover qualitatively agrees with a polytropic interior, with adiabatic index γ = 1.4, which one would expect from an inner convective interior. This is illustrated by the purple dotted curve in Fig. 2, where we have taken into account the smoothing of the potential, which is responsible for the flat slope at the smallest radii.

Overview
A large fraction of the gas inside the envelope, both in mass and volume, is not bound to the core. Instead, gas parcels enter and leave the planetary atmosphere in a complex fashion on orbital timescales. We find these flows reach a (nearly) steady state pattern and do not exceed the sound speed. Fig. 4 illustrates the characteristic orbits of gas parcels that enter the Hill sphere of the planetary core. At high latitudes above the disc midplane, fluid elements are only weakly perturbed by the planetary core, and perform the expected u-turn associated with horseshoe orbits. However, at lower latitudes, streamlines get bend significantly when they enter the Hill sphere and exit the atmosphere closer to the midplane (blue and green curves in Fig.4). Closest to the core, we note that field lines (red curve in Fig.4) remain trapped and perform a circumplanetary motion with a weak retrograde tendency.

Force balance
We now inspect the force balance on the fluid elements that approach the core. Assuming the Coriolis term dominates over advection (small Rossby number approximation), the momentum Fig. 4. Sample of gas streamlines, arriving from the upstream direction with respect to the planet (run RADL27k0.01). The core is located at x = r − a p = 0, y = a p (φ − φ p ) = 0, z = 0, as indicated by the black circle, and surrounded by a transparent purple sphere with radius r = r H . Results are shown in the frame co-rotating with the planet. Blue and green streamlines make a perturbed horseshoe orbit. Gas parcels arrive well above the midplane (z > r H ), enter the Hill sphere and subsequently exit closer to the midplane (z < r H ). The red and orange curves are examples of streamlines that interact strongly with the core. The orange streamline circulates closely around the core (within r < 0.5 r H ), while the red streamline indicates a gas parcel that remains trapped for the duration of the integration. equation in the rotating frame takes the form: Here we ignored the viscosity term and assumed steady state. This expression is similar to geostrophic balance, but with the addition of an effective potential term, which is the sum of the gravity terms and the centrifugal force (∇Ψ ). By multiplying Eq. 14 by ρ and taking the curl, we find the form We can further rewrite this expression by expanding both terms By making use of the steady state continuity equation, we can reduce this expression to Here, the directional derivative on the left hand term expresses the vertical gradients in the mass flux. Outside the Hill radius, densities are nearly uniform within a gas scale height from the midplane. Therefore the above expression reduces to ∂u/∂z = 0, which corresponds to the Taylor-Proudman theorem. As u z = 0 is zero in the midplane, it remains so at higher altitudes as well. Distant horseshoe orbits do not develop vertical motion and appear approximately identical from the midplane up to higher altitudes.
However, the planetary potential perturbs densities and the above approximation breaks down near the planet. From Eq. 17, it is nevertheless clear that a vertical mass flux gradient (∂(ρu z )/∂z) requires special non-aligned density gradients with 0.0 0.5 1.0 1.5 2.0 2.5 3.0 -10  Fig. 4, with panels A to D decreasing in height above the midplane. Magenta curves show the kinetic energy contribution, which is negligible. The purple curve is the effective potential term (see main text for full description). The green curve gives the thermal energy along the streamline. The black curve is the sum of the energy contributions, the Bernoulli parameter. Circles mark the points where the field line enters and exits the Hill sphere. The square symbol indicates the point where the fluid parcel crosses the x = 0 plane, which corresponds to the orbit of the planet. When multiple crossings occur the point closest to the core is selected. respect to the gravity force term, which are perpendicular to axis of rotation. Such a situation occurs when horseshoe-like streamlines come close to the planet and hit the spiral density perturbation triggered by the planet. Indeed, we note a sudden drop in altitude towards the second quadrant (x < 0, y > 0) for the blue, green and orange streamlines of Fig. 4. As Eq. 17 is general in nature, similar vertical motions are present in the isothermal approximation (Fung et al. 2015;Ormel et al. 2015).
We can also see that streamlines that get very close to the core and approach the midplane, rotate in a retrograde fashion. This hints that the pressure term slightly dominates the gravity terms towards the midplane, as can be seen from Eq.14. Such resulting retrograde motion, as seen in our prograde rotating frame, is reminiscent of vortices with high pressure centres in protoplanetary disc.
Closer to the core we can expect hydrostatic balance to set in. By definition this corresponds to the left hand terms of Eq. 14, or equivalently the left hand term of Eq. 17, being equal to zero. Assuming again no vertical velocities in the midplane, Eq. 17 implies that there can be no significant vertical mass flux into an innermost hydrostatic region.

Energy balance
As fluid elements enter the envelope they convert part of the liberated potential energy to thermal energy and exchange energy with the atmosphere. The energy budget on a streamline is given by for steady flow in a non-time dependent potential, while ignoring heat exchange. Here, the term (Ω × r) 2 /2 comes from the centrifugal force in the rotating frame and e is the internal energy per unit mass. The last two term of the expression form together the enthalpy. If one assumes no heat exchange, and takes the inviscid limit, the quantity C Be is conserved along the field line (Bernoulli's constant). However, heat is exchanged between the disc and the planet, through compressional heating and radiative exchange. We carefully investigate the energy balance for the streamlines depicted in Fig. 4. This analysis is presented with one panel for each streamline of Fig. 5. We first can note that the kinetic energy contribution is negligible for all streamlines. This may look surprising given the short timescales, on the order of 10 orbital timescales, on which fluid elements cross the Hill sphere, as one can read from the horizontal time axis. However, from an energetic viewpoint, they can be safely ignored. The potential term depicted here is the sum of stellar potential, the potential of the planet and the centrifugal term Ψ = Ψ grav − Ω 2 r 2 /2. For clarity, in Fig. 5 we only show the deviation from the background state at the orbital radius of the planet (Ψ back = −Ωr − Ω 2 r 2 /2 = −1.5 in code units).
As fluid elements approach the core, the combination of compressional heating and irradiation increase the thermal energy (E th = e + P/ρ). The peak value occurs in the point closest to the planetary core. The panels of Fig. 5 shows that in general internal energy increases as the potential component diminishes. However, this balance is not perfect, as one can see in the change in the Bernoulli parameter during envelope crossing. Fluid parcels loose gravitational energy as they approach the point closest to the core, put pick up energy on the way out of the envelope. Overall, the tendency for streamlines is to produce a net deposit of energy, as one can see by comparing the 6. Streamlines and gas density in the midplane of the planet, shown in a co-rotating frame centered on a planet. Here the core is accreting solids at rate of about 10 M E /Myr. (RUNHIGHRESL27). The x-axis points radially, the y-axis follows the azimuthal direction. Significant advection of mass occurs through the Hill sphere (marked by the white dashed line). The inner (left) and outer (right) circulating streamlines are shown in gray. Orange streamlines indicate the horseshoe region. The white streamlines show infalling gas from the poles being deflected outwards. Inside this region gas is unsteady in nature but nearly bound. A slow retrograde rotation can be seen, with respect to the rotating frame. The midplane gas density is shown in the background (in g/cm 3 ).
Bernoulli parameter at the entry and exit point of the Hill sphere (dot symbols in Fig. 5).
So far, we can conclude the following. The inspection of the averaged density and temperature structure reveals a planet which appears to be in near hydrostatic equilibrium, confirming our standard picture of planetary envelopes. In contrast, the analysis of the streamlines paints a more complex picture. There is no spherical symmetry. Moreover, streamlines reveal a strong steady state flow through the planetary envelope, which explicitly relies on breaking hydrostatic equilibrium. Little gas is bound to the envelope. Additionally, the flow leaves an imprint on the energy transport in the envelope, which we will investigate in more depth in Section 3.7.

Midplane
For completeness, we also show the flow morphology in the midplane in Fig. 6, as is often done (Tanigawa et al. 2012;D'Angelo & Bodenheimer 2013;Fung et al. 2015;Ormel et al. 2015). Care has to be taken to not over interpret these midplane streamlines as being a complete representation of the flow, because of the strong vertical flow component that is present (Tanigawa et al. 2012). Additionally, because we forced mirror symmetry, we force u z = 0 in the inner turbulent region where this condi- Fig. 7. Azimuthally-averaged specific angular momentum, normalized by the Keplerian value, measured in the midplane. Full lines indicate retrograde motion, while dashed lines mark prograde motion. The red curve shows an isothermal simulation (run ISO) with prograde rotation inside the Hill sphere. The black and blue curve correspond to, respectively, the reference simulations RADL27k0.01 and RADL0k0.01. Gray dotted lines indicate two limit cases, gas purely in Keplerian rotation around the core (h z /h K = 1) and gas in purely Keplerian rotation around the star, as if no planet was present (Eq. 22). tion does not necessarily needs to hold. In Fig. 6, we recognize the horseshoe orbits, depicted by the orange curves. The horseshoe orbits are strongly asymmetric across the y-axis, compared to standard 2D results (Fung et al. 2015). Keplerian shear enters deeply within the Hill sphere, as indicated here by the gray curves. The streamlines closest to the planet are unsteady close to the planetary core. This is in line with the notion that the inner envelope is convective. Further from the planet we note the two channels in the midplane that are aligned with the density wave perturbations launched by the planet, through which flow can escape from the planet and the horseshoe region. Fig 6 can be compared to the velocity field in an isothermal midplane, as depicted in Fig. C.2. The latter is characterized by two broad arms with flow escaping from the planet, as depicted by the white curves, in line with Fung et al. (2015) and Ormel et al. (2015). In our non-isothermal calculation, the width of this two arms is reduced, but the overall flow morphology is similar.

Angular momentum
Isothermal simulations show the emergence of prograde discs around planetary cores (Machida et al. 2008;Tanigawa et al. 2012). In our isothermal simulations we indeed recover prograde rotation (Fig. C.2). Gas is orbiting close to Keplerian velocity, consistent with the results by Tanigawa et al. (2012) (their Figure 8). In contrast, our non-isothermal simulations show the atmosphere does not build up significant angular momentum. Only a weak retrograde motion is measured, with respect to our prograde rotating frame. This is illustrated in Fig. 7, which shows the azimuthally averaged specific angular momentum, normalized by the Keplerian value, for our reference simulations (runs RADL27k0.01, RADL0k0.01, ISO). Here, r, φ are radial and azimuthal coordinates centered on the planet. We have thus found that envelopes are supported by a strong pressure force compared to the gravitational force, as opposed to results in the isothermal approximation. This change in the radial force balance altered the rotation direction from prograde to retrograde, as can be understood from Eq. 14. The planet thus only attracts low angular momentum accretion through the planetary poles. We do not see the development of a flattened disc, as the thermal energy contained in the envelope is not negligible compared to the gravitational potential. In fact, we find a result more consistent with Keplerian shear motion persisting throughout the envelope. In such a case, the angular momentum distribution takes the form, by ignoring azimuthal variations in the density. We can further linearise the velocity field to u y = − 3 2 Ω K x, as in the shearing sheet approximation, to find Finally, by making use of the definition of the Hill radius (Eq. 11), we obtain This relation is shown in gray in Fig. 7. The scaling is approximately maintained inside the Hill sphere, although averaged velocities do become slower closer to the core.

Mass flux through the envelope
The advection of gas is significant. Gas flows on tens of orbital timescales through the Hill sphere and involves a significant mass fraction of the envelope. Fig. 8 shows how mass is distributed in the envelope, by displaying the cumulative mass as function of the planetary radius. The total mass inside the Hill sphere is about 0.5 M E , about 10 % of the core mass, as shown by the black curve in Fig. 8. In our simulations we find the bulk of the envelope mass lays outside a radius of ≈ 0.2 r H . The steep slope in the cumulative mass function at smaller radii ( 0.1 r H ) is influenced by potential smoothing, as we could also see in panel B of Fig. 2. We now introduce a vertical mass flux timescale, where M(r) is the mass within a radius r, and the dominator is the inwards directed vertical mass flux through the surface S of the sphere. This is a mass-weighted measure of how rapidly gas is transported through the envelope. It corresponds to the total vertical flow, which both includes gas pushed in and again out of a spherical shell. This quantity is displayed in Fig. 9, as function of envelope radius. Outside ≈ 0.2 r H , we note gas is pushed through the envelope at a rate of about 10 % of the envelope mass per orbit, up to a Hill radius from the core. These rates are consistent with previous measurements of isothermal simulations, which are of the order of 1 % of the envelope per Ω −1 p (Ormel et al. 2015) to 100 % per Ω −1 p (Fung et al. 2015). At larger radii, outside the Hill radius, the vertical flux diminishes as gas comes in from directions perpendicular to to rotation axis. Closer to the core, around r 0.2 r H , we note a peak in t flux . It indicates that at this radius a large amount of the vertical mass flux is deviated away from the core. This transition is illustrated in more detail by the azimuthally-averaged vector field displayed in Fig. 10. This location broadly corresponds to the location of the radiative-convective boundary.
Article number, page 9 of 20 A&A proofs: manuscript no. paper <�> ��� � �� ��� ��� ��� Fig. 10. Azimuthally-averaged velocity field around an accreting 5 M E core, as seen in the co-rotating frame (run RADL27k0.01). The background shows the azimuthally-averaged temperature (in K) around the core. A large fraction of the gas that falls in from high elevation through the poles gets compressed above the inner planetary envelope, which gives rise to a mushroom-shaped compressional heat map. The inset shows in more detail the azimuthally-averaged velocity field, where overturning motions are driven by the heat release from solid accretion.
At even closer radii to the core, we do not measure an increase in t flux , which one would expect from a clear transition to a static bound envelope. Instead, an opposite trend is seen, where the vertical envelope mass flow becomes more efficient. We associate this with the emergence of convective flows, as can be seen in Fig. 10. At close distances to the core (r 0.05 r H ), the analysis is complicated by poor resolution and future work will need to address the inner envelope in more detail.

Heating and cooling a 3-layer envelope
The continuous flow of gas through the envelope also impacts heat transport. The outer shells, where gas flow is dominant, cannot be seen as parts of a static envelope that partake in secular cooling. Instead, they are dynamic regions, where a complex interplay of advection, heat deposition by compression (vertical direction) and cooling (horizontal direction) take place.
In order to further investigate the modes of heat transport, we consider two timescales: an advection timescale, and similarly a timescale to diffuse heat through radiation over the characteristic length scale of the planetary envelope, Outside the Hill radius, at heights above the scale height of the disc, heat exchange by radiative diffusion is efficient. Inside the Hill Radius, around 0.5 r H , we find a balance between advection and diffusion, as timescale become comparable and are of the order of 10-100 Ω −1 p . However, closer to the core, around 0.3 r H , we note that the diffusion timescale is much shorter than the advection timescale, which corresponds to a radiative zone in the envelope. At even closer distances to the core, within 0.2 r H , convective-like overturning motions become the dominant mode of energy transport, as heat diffusion timescales reach timescales of 10 4 Ω −1 p . Blue curves show results for the same planet, after the phase of solid accretion (run RADL0k0.01). Advection and radiative timescales follow a similar evolution with depth in the envelope. In the interior however, timescales increase by more than an order of magnitude.
We find a 3-layer envelope structure. In the outer layer, between 0.5 to 1 r H , advection of heat occurs on a similar timescale as radiative diffusion, as can be seen in Fig. 11. Outside of the envelope, at r 1 r H , radiative diffusion starts to dominate as we move the heights above the disc midplane.
In a middle layer, between 0.2 to 0.5 r H , we find that radiative diffusion is the primary heat transport mechanism. We can also identify this thin radiative shell in Fig. 10, as region with little advection.
In the inner shell, within r 0.2 r H , advection is again clearly the most efficient channel for heat transport. This is consistent with an interior envelope, where convective motions transport heat outwards.

Summary
So far, we have found that a planet in the process of accreting solids is not in hydrostatic equilibrium, but does reach a steady state flow pattern. The envelope acquires a characteristic 3-layer structure, made up of a convective interior, a radiative shell, and an advective outer envelope, which we illustrated with a cartoon in Fig. 12. We observe that the gas flux remains substantial, even through the deep interior, where the mass flux is on the order of 10 % of the enclosed mass per orbit (Fig. 9). At larger radii, rapid advection of gas and radiative diffusion both help transport heat away from the planet.
The envelope around the solid accreting planet we described is stable and would maintain steady state, were it not that on longer timescales the core mass increases and solid accretion rates change. When solid accretion dominates the luminosity budget, the cool-down process of the planet can safely be ignored. Then planetary envelopes move from one pseudoequilibrium state, similar to the one described here, to the next. However, when solid accretion rates become too low, or come to a complete standstill, this balance comes to an end. In the next section, we describe an envelope in this latter stage.

Results: after solid accretion
During the formation of a planet there is no guarantee solid accretion onto the core continues during the whole life time of the gas disc. Cores can get isolated from planetesimals due to scattering (Tanaka & Ida 1999) and resonant trapping (Levison et al. 2010). Similarly, pebble accretion comes to a halt when cores grow sufficiently massive to perturb the disc enough to from a pressure bump that prevents inwards drift and accretion of pebbles . Or pebble accretion simply comes to an end because radial drift depleted the available solid reservoir before the onset of gas disc dissipation (Sato et al. 2016). Therefore, we now focus on the same planetary core, placed in the same disc environment, but without interior heat generation by solid accretion (run RADL0k0.01).

Luminosity from compression
When the accretion of solids starts winding down, the dominant form of heat generation becomes compressional heating of gas that enters the planetary potential. Fig. 13 shows the azimuthallyaveraged compressional heating near the planetary centre and a heated cloud at the poles of the envelope where gas gets compressed as it gets deflected from the more bound interior towards the midplane. The integrated net heating, as function of the planetary radius, is shown in Fig. 14. We find a luminosity of about L ≈ 10 27 erg/s, from compressional heating.
Note that we have obtained this new equilibrium envelope with corresponding luminosity profile, without modelling the remnant interior heat in the envelope left by the formation of <���(ϵ ���� )> ��� -���� -��� -��� -��� -��� Fig. 13. Azimuthally-averaged velocity field around a 5 M E -core, similar to Fig. 10, but without the addition of accretion luminosity (run RADL0k0.01). The background shows the azimuthally-averaged compressional heating around the core. The inset zooms in on the vector field close to the core. A small amount of convergent motion can be identified, which can be placed in contrast with the inset in Fig. 10, which shows convective overturn.  the core and the secular cooling of the envelope. We have thus modelled a cooled-down planet that has lost this heat contribution, just before selfgravity becomes important. This latter stage in the evolution of the envelope is not accessible in our simulations, as we de not model self-gravity.

���(ρ) ���
-���� -���� -��� -��� -��� -��� Fig. 15. Streamlines and gas density in the midplane of the planet, shown in a co-rotating frame centered on a planet, not accreting any solids (run RUNHIGHRESL0). The x-axis points radially, the y-axis follows the azimuthal direction. Significant advection of mass occurs through the Hill sphere (marked by the white dashed line). The inner (left) and outer (right) circulating streamlines are shown in gray. Orange streamlines indicate the horseshoe region. The white streamlines show infalling gas from the poles being deflected outwards. Inside this region gas is bound and rotates slowly retrograde.
Compressional heating is also relevant even when solid accretion occurs. Fig.14 shows compression of gas generates a similar amount of luminosity as solid accretion. This energy is mostly liberated in the outer envelope. Therefore, the location and total amount of energy deposition between the runs with and without solid accretion are similar, but envelope structures are not identical.

Envelope structure
The thermodynamical structure of the envelope appears similar to the planet in the process of accretion solids, as can be seen in Fig. 2. The radial temperature profile is nearly identical. However, the density profile shows a sharper increase towards the centre, compared to the solid-accreting case. This central density increase results in a slightly more massive inner envelope, as can be inspected in Fig. 8.
Inspection of the midplane streamlines, as displayed in Fig. 15, also illustrates a similar flow structure as in the solidaccretion phase (Fig. 6). Rotation is again weakly retrograde, as also shown in Fig. 7. However the interior streamlines, do not reveal unsteady motion, instead we see bound field lines in the midplane.
Indeed, the inner envelope shows little motion. The interior is in near-hydrostatic balance and only a small amount of radial flow can be identified. This can be seen in Fig. 11, the advection timescale becomes exceedingly long as we approach the core. Nevertheless, the advection timescale is still shorter that the radiative diffusion timescale in the inner optically thick envelope (r 2 r H ).

Advection of gas
We again measure a strong mass flux through the outer envelope, of the order of 10% of the envelope mass per orbit, similar to the case of the solid accreting core, as shown in Fig.9. Closer to the core, within r 2 r H , t flux , the less turbulent interior leads to a slower vertical mass flux timescale. Previous studies have argued that this mass flux through the envelope can be estimated as a modified Bondi flow, where u char is a characteristic velocity, often taken to be the Keplerian shear and the cross section is typically taken to be the Bondi radius (D'Angelo & Lubow 2008;Lissauer et al. 2009;Ormel et al. 2015). For the planet we consider, this corresponds to a mass flux of the order 10 3 M E /Myr, which is indeed in reasonable agreement with our measurement of the flow through the outer envelope (Fig.9). Previous studies (D'Angelo & Lubow 2008;Lissauer et al. 2009;Machida et al. 2010) have interpreted this mass flux as a measurement of the accretion rate, implying unrealistic fast envelope growth timescales. Instead, we have show here that most of this mass flux is caused by an unbound flow through the envelope.

Accretion of gas
A direct measure of the gas accretion rate, the net gas flux that remains bound to the core, is difficult to make. The flow of gas that transits the envelope dominates the much smaller fraction that is accreted. Directly measuring the mass growth of the envelope is difficult as this occurs on timescales that are numerically not feasible. In Appendix A, we show results of the longest integrations which are numerically feasible, of the order of 50 orbits. A accurate measure of any direct mass growth of the planetary envelope cannot be made. We can however indirectly place a constraint on the gas accretion rate. The generation of heat requires a net divergence of the flow (−P∇ · u), which we can link to the accretion flow (inset of Fig. 13). As the luminosity generated by compressional heating is ultimately the liberation of potential energy, we can estimate the gas accretion rate aṡ where r char is a characteristic radius, corresponding to depth of the potential. Taking L ∼ 5 × 10 26 erg/s, and r char ∼ r core we would find an accretion rate of the order of aboutṀ gas ∼ 1 M E /Myr. This slow accretion estimate is in line with 1Dmodels (Piso & Youdin 2014), and further illustrates the challenge of making direct measurements of the gas accretion rate.

opacity in protoplanetary discs
Inside the planet-forming regions of protoplanetary discs, dust particles smaller than mm in radius are the dominant cause of absorption and scattering of radiation. The opacity caused by dust can be expressed as, just outside the ice line at approximately 170 K. Here, we have taken a composition and dust-to-gas ratio Z dust similar to the interstellar medium (ISM) from which the star and protoplanetary disc formed. Inside the ice line, icy particles sublimate and opacities rapidly settle to that provided by silicate grains, A more complete description can be found in Bell & Lin (1994); D'Angelo & Bodenheimer (2013); Semenov et al. (2003). These values can best be interpreted as upper limits. As dust grows by collisions into larger particles, the ratio of small dust to gas, Z dust will decrease. Nevertheless, in special some regions in the disc, for example around ice lines, the complex interplay between growth, fragmentation and sublimation may lead again to higher values.
Inside the envelopes of planets, the evolution of the dust opacity is poorly constrained. As a first approximation, one could assume the envelope inherits the same opacity as the disc, valid when coagulation proceeds more slowly then advection and settling of dust inside the envelope. Alternatively, one can estimate how the dust opacity may be reduced by grain growth inside the planetary envelope (Movshovitz et al. 2010;Ormel 2014;Mordasini 2014). These models find atmospheres with ISM-like opacities in the outer parts of the envelope, that decrease radially inwards to values around 10 −3 cm 2 /g. Nevertheless, these models do not take into account bouncing and fragmentation of grains, as well as the large gas advection flows through the atmosphere, that will incorporate dust grains.
In the deep interiors of planetary envelopes, within the sublimation line for dust T dust ∼ 1000 K, opacities are set by the gas component. Such high temperatures are only reached very close to the core, within about 10 core radii (Piso & Youdin 2014;). Our simulations do not resolve such short distances. Therefore, we do not reach dust sublimation temperatures in the envelope, as can also be seen in the top panel of Fig. 2. Consequently, we are not required to include gas opacities in our calculations.
Given these uncertainties in dust opacity values, we have limited this study to opacities that are constant throughout the envelope. The previous results we showed were obtained for envelopes with κ = 0.01 cm 2 /g, which corresponds an envelope with moderate grain growth (as used in standard woks such as Hubickyj et al. 2005). We now will present envelopes where we assume no significant growth of dust has occurred and opacities are close to an ISM-like opacity for silicate grains (κ = 1 cm 2 /g). These unreduced opacities thus assume no significant dust depletion. We will again consider two cases, the case of a planet accreting solids (run RADL27k1), as well as the case where solid accretion has come to a halt (run RADL0k1).

Solid-accreting planet with unreduced opacity
In a higher opacity environment (κ = 1 cm 2 /g), gas flows through the entire envelope of a solid-accreting planet (L acc = 10 27 erg/s). This can be seen in Figure 16, which shows the azimuthallyaveraged velocity field and can be compared directly to Figure 10. Clearly, there is no trace left of the 3-layer envelope dis- Fig. 16. Azimuthally-averaged velocity field around an accreting planet in a high opacity environment (run RADL27k1). The gas flux from the poles enters even the deep interior, as is clarified in more detail in the inset. This stands in contrast to the low opacity case (run RADL27k0.01, Fig. 10), where rapid advection was limited to the outer envelope. cussed in Sec. 3.7. Instead gas circulates through the planet even at depths close to the core. There is thus no development of a more bound interior shielded by a thin radiative shell where t adv > t diff , as was the case in the reduced opacity models. Instead, radiative diffusion timescales remain of the same order as the advection timescales, as can be seen from the red curves in Fig. 17. Radiative diffusion is overall much more efficient, compared to the low opacity case, presented in Fig 11 (run RADL27k0.01). It might appear counter-intuitive that higher opacities promote more efficient radiative diffusion. However, the improved radiative transport is mostly just a consequence of the lower central densities. Indeed, the envelope has a much smaller mass in the interior, as can be seen in Fig. 8.
It is worthwhile to reflect on this unexpected critical transition obtained by simply increasing the opacity. A new equilibrium state emerged, where the envelope is dominated by the advection of gas that enter from the poles. This stands in contrast to spherical hydrostatic models which would have predicted that increasing the opacity would result in the envelope becoming simply more convection dominated. In such models the depth of the radiative zone is approximately inversely proportional to the opacity (P rcb ∝ (κL acc ) −1 , Eq. D.3). Indeed, we verified that this scaling approximately holds in the low κ-regime by recovering equivalent envelopes between run RADL27k0.01 and a run with κ = 0.1 cm 2 /g and L acc = 10 26 erg/s, which holds the product κL acc constant. However, as we have shown here, when the opacity is increased enough the radiative outer zone starts to shrink up to point where a sudden transition occurs. Then the envelope opens up and becomes advection-dominated.
Article number, page 13 of 20 A&A proofs: manuscript no. paper  Fig. 17. Advection and thermal diffusion timescales through the envelope, similar to Fig. 11, but for an unreduced opacity κ = 1 cm 2 /g. Red curves show results from run RADL27k1, which is the model of a solidaccreting core. Advection timescales are short, of the order of 10 Ω −1 p throughout the envelope. The case without energy deposition by solid accretion (run RADL0k1, green curves) shows advection and diffusion time scales are of the same order and peak towards the inner envelope where they become as long as 10 4 Ω −1 p .
Finally, we verified this result once more by analyzing a simulation with an initial condition that was different from the standard procedure obtained by gradually introducing the planet in the disc, as outlined in Sec. 2.4. Instead, the equilibrium state of RADL0k1 was used as the initial condition. After initialization, the envelope gradually got heated by accretion heating and we found an equivalent envelope structure this way as well. Therefore the fast advection result is not dependent on the method used to initialise the simulation.

After solid accretion, with unreduced opacity
We now consider the case where solid accretion has come to a halt in an unreduced opacity environment (κ = 1 cm 2 /g), corresponding to run RADL0k1. We recover an envelope similar to the case with lower opacity (run RADL0K0.01). Diffusion and advection time scales throughout the envelope are indeed similar, as can be seen from comparing the diffusion and advection time scales (green curves in Fig. 17) with those shown previously for run RADL0K0.01 in Fig. 11.
The luminosity released by compressional heating in the interior envelope is approximately L comp = 5×10 25 erg/s, as shown by the green curve in Fig.18. This is an order of magnitude lower compared to the case with κ = 0.01 cm 2 /g. This reduction can be understood approximatively, because we recover a similar depth of the radiative-convective boundary r rcb /r B , or equivalently P rcb /P mid when we force no accretion heat. From Eq. D.3, we find P rcb /P mid ∝ κ −1 L −1 P −1 mid T 4 mid (for more details see Appendix D). Therefore, one would expect L ∝ 1/κ and a reduction of the luminosity by a factor 100. However, changing the opacity also changes the disc structure, and therefore the values for the temperature and pressure. This softens the reduction in the Integrated energy delivery and removal by respectively compression and expansion, for the unreduced opacity environment (κ = 1 cm 2 /g). Red curves show the case where solids are accreted (run RADL27k1). Inside the Hill sphere, we dominantly see that inner envelope mostly acts as a sink of energy that helps gas parcels to expand as they move through the envelope. Green curves represent the energy deposition from compressional heating after solid accretion (run RADL0k1). Compressional heating now provides the luminosity that supports the envelope. Note the y-axis has been cut to represent the negative values on the log axis.
luminosity to a factor 10. It also explains why the deposition L acc = 10 27 erg/s onto the same 5 M E -planet with k = 1 cm 2 /g is a significant contribution that drastically changes the envelope structure, as discussed in Sec. 5.2. The interior of the envelope does not longer generate a supporting luminosity from compressed gas, instead gas heats up and expands leading to a net cooling effect, as can be seen from the red curve in Fig. 18. The lower luminosity we measure implies slower gas accretion rates. Repeating the order of magnitude analysis in Sec.4.4, we would find a gas accretion rate ofṀ gas ∼ 0.1 M E /Myr, for a luminosity of L ∼ 5 × 10 25 erg/s, and r char ∼ r core . This would lead to the formation of planets that have envelopes up to ∼10% of the total planetary mass, similar to the known Super-Earths and ice-giants.

Dependency on core mass
The mass of the core is, besides the accretion rate and the opacity, an important quantity that regulates the structure of the envelope. We have now seen that in environments where it is difficult for heat to diffuse out of the envelope, the atmosphere can get advection-dominated. We identified this effect around a lowmass core of 5 M E . Increasing the mass of the core should suppress this effect, as the potential of the planet becomes too deep to allow the flow of gas to go through unhindered.
recover a flow profile that shows great similarity to the one we observed around a 5 M E -core, in low opacity environment. Fig. 19 shows the azimuthally-averaged flow profile around a core in the process of accreting solids. Clearly, we do not longer see the large scale advecting flows seen in Fig. 10 which occurred around the 5 M E -core in a similar κ = 1 cm 2 /g environment. The envelope recovers an exterior radiative-advective region, with interior a shell where radiative transport dominates. Close to the center overturning motions become the dominant mode for energy transport. This can also be seen when inspecting the advection and diffusion timescales presented in Fig. 20.
Also, without the contribution of solid accretion, we note that gas accretion speeds up again (run RADL0k1M15). In this case, we measure a luminosity from compressional heating of approximately 5×10 26 erg/s. Following the order of magnitude analysis of Eq. 27, we find a gas accretion rate of around 1 M E /Myr, about an order of magnitude larger than around the 5 Earth mass core in a similar unreduced opacity region. This facilities the transition to runaway gas accretion, but the growth rate is still modest and disc dissipation might occur before growing a massive envelope.

Implications: delaying runaway gas accretion
Previous studies (Mizuno 1980;Pollack et al. 1996;Ikoma et al. 2000;Rafikov 2006;Piso & Youdin 2014) investigated the role of the core mass, opacity, solid accretion rate and envelope composition, but assumed hydrostatic envelopes. Here, we have seen that envelopes are not in hydrostatic balance, and large flows of gas enter and leave the planetary atmosphere. This result has several significant implications on the growth of the envelope.

Dust opacity
Previous studies have speculated that envelopes could acquire atmospheres up to 100 times enriched in opacity-providing grains, with respect to the star and the disc (Lee et al. 2014;Venturini et al. 2016). Such high opacities would delay the contraction of the envelope. However, this pathway can now be dismissed on three grounds. High dust-to-gas ratios do not remain homogeneously mixed (Lambrechts et al. 2016), coagulation aids the rain out of solid from dust-rich atmospheres (Ormel 2014;Mordasini 2014), and here we add that the rapid flux of gas through the atmosphere does not allow for such dust pile ups in the first place. Conversely, extremely low dust opacities are also ruled out, because fresh dust is continuously transported through the outer envelope. Therefore, we expect dust opacities in the radiative-advective outer region of planets to be, to a good approximation, comparable to the dust opacity in the disc at the location of the planet. A more precise determination of the opacity at the radiative-convective boundary would require to solve selfconsistently for the transport of dust particles and the opacity dependent depth of the advective outer shell.

Accretion of solids and envelope pollution
We expect the gas flows through the envelope to only have a moderate effect on the accretion rate of solids. Large planetesimals feel little gas drag, because of their large size. Smaller pebbles on the other hand could see their trajectories changed. However, accretion of pebbles dominantly occurs close to the disc midplane, from a direction diagonal to the direction of outflow along the spiral density waves (Lambrechts & Johansen 2012). Additionally, pebbles have friction and crossing times that are comparable, of the order of the orbital timescale, while the flux through the atmosphere occurs on 10 times as long timescales in the outer envelope. Nevertheless, we plan to study this in more detail in an upcoming work, which will also allow us to study at which depths in the envelope the accretion heat is released.
As solids sink in the envelope they can evaporate and pollute the envelope. Particles with a silicate composition only vaporize at high temperatures exceeding 1000 K that are found close to the core. Therefore, these dust species settle in the inner nearlyhydrostatic envelope and will be accreted. Ices on the other hand can deposit water vapour at much shallower depths in the envelope. If the sublimation point is located in the highly advective outer regions of the atmosphere, vapour will simply be transported away in tens of orbital timescales. However, in our simulations where the planet is located outside of the water ice line of the disc, we can see that this is not the case. In fact, the water ice sublimation point, at about 170 K, is located deep inside the convective interior, shielded from the outer layers where rapid gas advection takes place (Fig. 2). Therefore, not all water vapour will be immediately lost. Instead, the solid accretion rate, the redistribution of vapour by convection and the advection rate will set the fraction of solids that settle to the core, the fraction that is deposited in the envelope and the vapour fraction that is lost. It will be important to quantify this in more detail, because this balance influences core growth and the accretion of gas as well, since the deposition of high mean molecular weight material speeds up gas accretion (Hori & Ikoma 2011;Venturini et al. 2016).

Cooling the envelope
Finally, there is the impact on the secular cooling of the envelope. Gas in spherical hydrostatic models is bound. Therefore the outwards transport of heat is paired with cooling the envelope and increasing the envelope mass. In this way, a hydrostatic growth sequence can be constructed for planets (Pollack et al. 1996;Ikoma et al. 2000;Piso & Youdin 2014).
Rapid advection of gas through the envelope challenges this paradigm (Ormel 2014;Fung et al. 2015). The outward transport of heat is not extracted from a bound reservoir of gas, which slows down envelope cooling and associated mass growth compared to spherical hydrostatic models. However, the picture is complicated. After the solid-accretion phase, we have shown that planets have their strong advecting flows limited to the outer envelope (Sec. 4). In the inner envelope, we see mass advection diminish as we get closer to the core. Therefore, gas close to the core is in approximate hydrostatic balance and modified 1D prescriptions may be suitable, allowing an approach as in D 'Angelo & Bodenheimer (2013).
A more quantitative analysis determining the bound envelope will require an improved resolution in the inner shells. Resolving this region close to the core would also necessitate the use of a more accurate, non-constant, opacity prescription that takes the high temperatures in the inner envelope into account. An improved equation of state may also be necessary. Such improved measurements would aid in making secular cooling histories of young planets more accurate. Furthermore, it would be of great value to explore also a larger planetary mass range, since even planets as large as Jupiter see a large gas flux transit through their atmospheres (Gressel et al. 2013;Morbidelli et al. 2014;Szulágyi et al. 2014).

Summary
In summary, we have identified three different effects that stall envelope accretion. Dust opacity values will be similar to the disc, water pollution can be inefficient, and cooling times are delayed by the delivery of fresh gas.

Interpreting exoplanet Super-Earths and gas giants
Why do some planets build up massive envelopes and others do not? We have argued that the formation of gas giants with large atmospheres requires specific conditions: the rapid formation of cores larger than 15 M E in a disc environment with reduced opacities. These requirements can be met in the outer regions of the protoplanetary disc, where particles have grown to cm-sizes and radially drift inwards through the disc. Such rocky/icy pebbles can be efficiently accreted by planetary embryos (Ormel & Klahr 2010;Lambrechts & Johansen 2012), allowing for approximate 20 M E -cores to form outside of the iceline within disc lifetimes . Such fast core growth does generally require initial dust-to-gas ratios in the disc that are solar-like or higher Bitsch et al. 2015). This requirement seems to be in line with the observed lack of close-in gas giant planets around stars with sub-solar metallicities (Buchhave et al. 2012). We note that the required solar or higher initial dust-to-gas ratio does not imply that densities of opacity-providing grains must have been high as well. Particle growth outside the ice line reduces the opacity and local pebble surface densities fall rapidly by a factor 10 below initial values because of radial drift . Super-Earths and ice giants did not accrete such massive gas envelopes. Their smaller core masses helped prevent rapid gas accretion. Additionally, they may have formed dominantly in regions where disc opacities were not significantly reduced.
Ice giants have cores that formed slowly in very wide orbits and likely never experienced strongly reduced opacities, as dominant particle sizes remain small in the outer disc (Brauer et al. 2008).
Super-Earths on the other hand are commonly found in short orbits. They thus formed in situ or the cores migrated into this inner disc region. Possibly, they dominantly formed in a part of the disc where solids piled up, which would explain why Super-Earths multiples are common (Winn & Fabrycky 2015) and occurrence rates are only weakly dependent on stellar metallicity (Buchhave et al. 2012). Such a pile-up of solids would likely have been associated with a high opacity environment, because of dust delivery by radial drift and possibly collisional fragmentation, which would further limit gas accretion.
In summary, we propose that gas giant planets dominantly form outside of the ice line, where particles grew to pebble sizes, allowing the formation of large cores. Models of exoplanet compositions support that gas giants typically have substantial cores exceeding 10 M E (Thorngren et al. 2016). Reduced opacities may then have further facilitated rapid gas accretion. Super-Earths on the other hand have smaller cores, and may additionally have formed in high opacity environments associated with pile-ups of solids in the inner disc. This prevented rapid gas accretion. The specific conditions for the formation of massive gas envelopes may explain in part why occurrence rates for gas giants are below Super-Earths by about a factor five (Winn & Fabrycky 2015).

Summary
Atmospheres around growing planets in the Super-Earth to ice giant regime are not spherically symmetric envelopes in hydrostatic balance. Here, we presented global 3D radiative hydrodynamical simulations that show a steady state flow through planetary envelopes. Large amounts of gas enter through the poles of the envelopes and exit near the disc midplane. This is the result of density perturbations with gradients that are not aligned with the gravitational potential.
Our results imply that (1) opacities in the outer envelope are similar to disc opacities, (2) icy solids will sublimate in the convective interior and (3) the flow of gas delays the secular cooldown of the envelope. We numerically investigated the dependency of our results on three key parameters; the accretion rate, the dust opacity and the mass of the core.
Dust opacities can fluctuate widely in the protoplanetary disc, because of particle growth and drift. In regions of the disc with dust opacities of κ = 0.01 cm 2 /g, below the initial ISMlike opacity of about κ ≈ 1 cm 2 /g, planets develop a 3-layer atmosphere. The envelope consists of an advective outer shell, a radiative shell around a radius of 0.3 Hill radii, and an interior which shows convection-like overturning motions that transfer the heat of solid accretion. When solid accretion comes to a halt, gas close to the core slowly settles in the potential. The liberated heat from this convergent motion supports the envelope. Measurements of the accretion rate are challenging, but our estimates are of the order of 1 M E per Myr. Improved measurements will require long-term simulations exceeding the advection timescale which are performed at high resolution in order to allow smoothing lengths down to the core surface.
We also explored unreduced dust opacities that are comparable to values found in the ISM (κ = 1 cm 2 /g) in the disc and envelope, motivated by the the high gas advection rates in the outer atmosphere. For planets in the process of accreting solids, we identify a previously unreported transition into an advectiondominated envelope. Hydrostatic models would have predicted a nearly completely convective atmosphere. Instead, we found an envelope where advection of gas, from the poles to the deep interior, everywhere becomes the dominant mode of energy transport. After solid accretion comes to an end, gas accretion rates are low, of the order of 0.1 M E per Myr. Under these conditions, a 5-M E core would not go into runaway gas accretion and therefore remain a Super-Earth.
However, larger cores lead to more bound interior envelopes and an increase in the gas accretion rate. We found that around massive solid-accreting cores (M c 15 M E ) envelopes are no longer advection-dominated in a κ = 1 cm 2 /g environment. After solid accretion, the increase in the core mass boosts the envelope growth rate by a factor 10, compared to the 5 M E case. Therefore, a transition to runaway gas accretion becomes possible again.
These findings deserve further attention and future models of envelope evolution would benefit from a more quantitative analysis. However, based on our current results we can already argue that Super-Earths did not accrete massive gas envelopes because their cores are sufficiently small, provided they formed in disc environments where dust opacities were not strongly reduced below κ ∼ 0.01 cm 2 /g. We have verified that we obtain steady state envelope structures. Figure A.2 shows the evolution of the envelope mass with time. We find an envelope structure that remains stable on timescales of around 50 orbits. We also verified this by inspecting the stability of the total torque. The long term integrations presented in Fig. A.2, do not reveal any accretion of gas. Even longer simulations, exceeding 10 5 years, would likely be necessary in order to capture the secular cooling process of the envelope and the resulting growth in envelope mass. Inspection of Fig. A.3 shows that about 8 orbits at medium resolution are sufficient to capture the thermal structure and advection flow inside the planetary envelope to sufficient precision. However, we do note that high resolution simulations are required to capture the flow structure around r ≈ 0.2 r H . Therefore, about 4 orbits at high resolution are required to capture the effect of a reduced vertical flux. Demonstrating full convergence in t adv would demand running a high resolution simulation for a time t t adv , which is currently unfeasible. Nevertheless, since no qualitative changes are expected to occur physically at later times, a shorter run time is sufficient in the context of this work.

Appendix B: Methods: Accretion luminosity
We use a simple prescription to model the accretion heat released at the interior of the planetary envelope (Benítez-Llambay et al. 2015). All accretion heat L acc is assumed to be deposited uniformly in the volume prescribed by the inner grid cells around the planet core. Simulations are mirrored across the midplane, so we only prescribe the 4 grid cells of the upper midplane that are centered around the midpoint of the planetary potential. For each of the 4 cells, we set an energy deposit rate of l = L acc 8V cell . (B.1) Here, V cell is the grid cell volume, which is well described by a constant, since our non-uniform grids are to a very good accuracy equal-sized close to the position of the planet. Convergence of radiative and advective timescales. Radiative timescales in run RADL27k0.01LONG converge rapidly in about 4 orbits, outside the smoothing radius at 0.25 r H . Little change occurs between 8 to 50 orbits. After 8 medium resolution orbits, the restart at high resolution, labeled with hr in the legend, does not alter t diff outside 0.1 r H . Advective timescales also quickly settle. Note however that approximately 4 orbits at high resolution are required to capture the increased vertical advection timescale around 0.2 r H (run RADL27k0.01).

Appendix C: Methods: Isothermal Equation of state
In this section, we compare FARGOCA using an isothermal equation of state, with previous results obtained in the literature for the 3D flow in low viscosity discs. In this way we test the code

Fig. C.2.
Velocity field around a 5 M E -core in the co-rotating frame (run ISO). Displayed are the streamlines in the midplane (z = 0). One can identify the outer regions dominated by Keplerian shear, the horshoe orbits in front and behind the planet, and two distinct arms transporting mass away from the planet. The background shows the gas density in the midplane. and more specifically the implementation of the 3D non-uniform grid.
We first compare the morphology of the flow field in the disc midplane around the planet. We have used a parameter set very close to the one reported by Fung et al. (2015) (simulation ISO). Figure C.2 and C.1 show the motion of gas around the planet. We find the same morphology as reported in Fung et al. (2015), for their 3D results in a low viscosity disc, and Ormel et al. (2015). Gas enter vertically the Hill sphere ( Fig. C.1) and escapes along the midplane through two channels in the second and fourth quadrant (white streamlines in Fig. C.2). This represents the transient horseshoe flow reported in Fung et al. (2015).

Appendix D: Hydrostatic model of a spherical atmosphere
In this appendix, we briefly discuss the structure of a planetary envelope in hydrostatic balance (Mizuno 1980;Rafikov 2006;Piso & Youdin 2014). In the outer shell, energy transport occurs through radiation. The temperature gradient with respect to the pressure takes the form ∇ rad = d ln T d ln P = 3κ(r)P(r) 64πGσ SB T (r) 4 L(r) M int (r) . (D.1) We now will assume a constant opacity and luminosity through the envelope, L(r) = L. We consider low mass envelopes, therefore the potential is dominated by the core, M int (r) = M c . Then the temperature remains nearly constant in the outer envelope, T (r) = T out . Thus the pressure and density are exponentially increasing, P = P out exp r B 1 r − 1 r out .

(D.2)
Here P out is the pressure at an unperturbed point at r out , often taken to be the Bondi radius. The temperature gradient cannot grow unbounded towards the interior. Steep temperature gradients become unstable to convection. Energy transport starts occurring through overturning motions of gas, that are assumed to maintain a critical temperature gradient ∇ ad = γ−1 γ . Therefore the depth of the radiative zone is set by the point where ∇ ad = ∇ rad . This allows us the rewrite Eq. D.1 to find the pressure at the radiative convective boundary, Here the subscript rcb stands for values evaluated at the point where outward radiation transport starts, at the radiative convective boundary. The depth of the radiative zone can be expressed, using Eq. D.2, as This expression represents a powerful relationship between the depth of the radiative zone, the opacity and the luminosity of the planet. Note also that the depth of the radiative zone uniquely determines the energy stored in the envelope. Hot, nearly fully convective envelopes, have r rcb close to the outer boundary edge of the atmosphere. Cooling down the envelope goes paired with the radiative zone entering progressively deeper in the envelope.

Appendix E: Pebble accretion
In Fig. E.1, we present the accretion luminosity for a core growing by the accretion of pebbles, in a nominal pebble disc model . The efficient accretion of dragsensitive particles in the cm size range allows cores to grow to completion before gas disc dissipation. This growth rate can be considered to give an upper limit on the accretion luminosity, because pebbles be can accreted from the full Hill sphere and the whole mass reservoir of inwards drifting pebbles is available. Conversely, accretion rates cannot have been much lower in order to meet the constraint that cores form before disc dissipation. The red curve in Fig The black curve corresponds to the luminosity by pebble accretion, which can be considered the maximal possible accretion rate. The red dashed line gives a lower bound on the accretion luminosity from the requirement order to have sufficiently high accretion rates to finish core growth within the lifetime of the protoplanetary disc. The orange squares represent the accretion luminosities explored in this work.