The population of closein superEarths with gas mass fractions of up to 10 per cent represents a challenge for planet formation theory how did they avoid runaway gas accretion and collapsing to hot Jupiters despite their core masses being in the critical range of MSUBcSUB  10 MSUBSUB Previous threedimensional 3D hydrodynamical simulations indicate that atmospheres of lowmass planets cannot be considered isolated from the protoplanetary disc contrary to what is assumed in 1Devolutionary calculations. This finding is referred to as the recycling hypothesis. In this paper we investigate the recycling hypothesis for superEarth planets accounting for realistic 3D radiation hydrodynamics. Also we conduct a direct comparison in terms of the evolution of the entropy between 1D and 3D geometries. We clearly see that 3D atmospheres maintain higher entropy although gas in the atmosphere loses entropy through radiative cooling the advection of highentropy gas from the disc into the BondiHill sphere slows down KelvinHelmholtz contraction potentially arresting envelope growth at a subcritical gas mass fraction. Recycling therefore operates vigorously in line with results by previous studies. However we also identify an inner core  in size 25 per cent of the Bondi radius  where streamlines are more circular and entropies are much lower than in the outer atmosphere. Future studies at higher resolutions are needed to assess whether this region can become hydrodynamically isolated on long timescales.
['hydrodynamics', 'planets and satellites: atmospheres', 'planets and satellites: formation', 'protoplanetary discs', 'Astrophysics - Earth and Planetary Astrophysics']
Hydrodynamics of embedded planets first atmospheres  III. The role of radiation transport for superEarth planets
{"Hydrodynamics of embedded planets' first atmospheres III. The role of radiation transport for super-Earth planets": 'Nicolas P. Cimerman, 1 ? Rolf Kuiper, 1 Chris W. Ormel 2  \nInstitut fur Astronomie & Astrophysik, Universitat Tubingen, Auf der Morgenstelle 10, D-72076 Tubingen, Germany 2 Anton Pannekoek Institute, University of Amsterdam, PO Box 94249, NL-1090 GE Amsterdam, the Netherlands  \n- 1  \nAccepted 2017 July 22. Received 2017 July 18; in original form 2017 June 7', 'ABSTRACT': "The population of close-in super-Earths , with gas mass fractions of up to 10% represents a challenge for planet formation theory: how did they avoid runaway gas accretion and collapsing to hot Jupiters despite their core masses being in the critical range of M c ' 10 M GLYPH<8> ? Previous three-dimensional (3D) hydrodynamical simulations indicate that atmospheres of low-mass planets cannot be considered isolated from the protoplanetary disc, contrary to what is assumed in 1D-evolutionary calculations. This finding is referred to as the recycling hypothesis. In this Paper we investigate the recycling hypothesis for super-Earth planets, accounting for realistic 3D radiationhydrodynamics. Also, we conduct a direct comparison in terms of the evolution of the entropy between 1D and 3D geometries. We clearly see that 3D atmospheres maintain higher entropy: although gas in the atmosphere loses entropy through radiative cooling, the advection of high entropy gas from the disc into the Bondi/Hill sphere slows down Kelvin-Helmholtz contraction, potentially arresting envelope growth at a sub-critical gas mass fraction. Recycling, therefore, operates vigorously, in line with results by previous studies. However, we also identify an 'inner core' - in size GLYPH<25> 25 % of the Bondi radius - where streamlines are more circular and entropies are much lower than in the outer atmosphere. Future studies at higher resolutions are needed to assess whether this region can become hydrodynamically-isolated on long time-scales.  \nKey words: hydrodynamics - planets and satellites: atmospheres - planets and satellites: formation - protoplanetary discs.", '1 INTRODUCTION': "Giant planets form in gas rich discs either by gravitational instability (Kuiper 1951; Cameron 1978) or, as is assumed more commonly, by core accretion (Safronov 1969; Pollack et al. 1994). In the latter scenario the final rapid gas accretion proceeds only after a massive enough 'core' has been assembled. A typically quoted number for this critical core mass (CCM) is GLYPH<24> 10 Earth mass (Stevenson 1982; Pollack et al. 1994) - a number that follows from solving 1D stellarstructure like equations. Assuming radial symmetry, these models compute how the temperature, pressure, and density profiles evolve as function of time. As pointed out by Piso & Youdin (2014), these models come in different flavours: static models consider accretion of solids to provide a steady substantial luminosity to support the envelope (e.g. Stevenson 1982; Rafikov 2006), while in quasi-static models hydrostatic equilibrium is assumed at any given time, but the KelvinHelmholtz contraction of the atmosphere is considered as  \n- ? E-mail: npcphys@gmail.com (NPC)  \nsource of the luminosity (e.g. Ikoma et al. 2001; Lee & Chiang 2015).  \nThe GLYPH<24> 10 M GLYPH<8> for the CCM is typical for the formation of Jupiter assuming ISM-like grain opacities (in regions cool enough for silicate grains to be present), a steady influx of planetsimals, and a H/He-dominated envelope (all solids end up in the core) (Pollack et al. 1994; Hubickyj et al. 2005; Mordasini et al. 2009). However, in the last years it has been realized that some of these model parameters are chosen rather conservatively. In envelopes, grains will quickly coagulate and settle (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014). A zero grain opacity may be more realistic, especially in the case without solid accretion. It is also rather restrictive to assume continuous planetesimal accretion. Without an external source for luminosity and grain opacity Hori & Ikoma (2010) calculate that cores as low as 2 M GLYPH<8> will collapse into a gas giant after 1 Myr. A metalrich composition of the envelope will also reduce the CCM (Venturini et al. 2015).  \nThese considerations are especially important for superEarth planets. Super-Earths are a class of planet that lie by  \nmass in between terrestrial ('rocky') planets and giant planets. They have been found in great numbers by transit and radial-velocity surveys, in particular by the Kepler spacecraft (Fressin et al. 2013; Petigura et al. 2013), at short distances (typically GLYPH<24> 0.1 au) from their host star. With radialvelocity (RV) follow-up (Mayor et al. 2011; Marcy et al. 2014) or transit timing variations (TTV) (Hadden & Lithwick 2014; Jontof-Hutter et al. 2016; Mills & Mazeh 2017) a bulk density can be obtained, which, for the more massive super-Earths, is low enough to infer that these planets contain significant amounts ( GLYPH<24> 1% by mass) of hydrogen and helium (Lopez & Fortney 2014). These super-Earths are also referred to as mini-Neptunes. With these hydrogen/helium mass fractions, it is plausible to assume that these planets formed early, in a gas-rich disc, like giant planets do. Even the current rocky super-Earths ( i.e. , those planets of high bulk density) could have started with gas-rich envelopes, only to have their atmospheres stripped by photoevaporation (Lopez & Fortney 2013, 2014; Owen & Wu 2013) or lost in a 'boil-off' event, when the confining pressure of the disc is gone after dispersal (Owen & Wu 2016). A further analysis based on more precise radii from the CaliforniaKepler Survey by Fulton et al. (2017) provides some evidence for such a scenario which has been recently discussed by Owen & Wu (2017) and Jin & Mordasini (2017).  \nWhen super-Earths did assemble before the disc dispersal, the question is how they avoided breaching the CCM, which would have turned them into hot-Jupiters. A key point is that at GLYPH<24> 0.1 au planetesimals or any other solid particles are accreted on time-scales GLYPH<28> Myr (Chiang & Laughlin 2013). Modelling super-Earths envelopes including non-ideal effects (like hydrogen dissociation and silicates evaporation), Lee et al. (2014) calculated that the Kelvin-Helmholtz contraction took GLYPH<24> Myr. To avoid envelope collapse Lee & Chiang (2015) suggested super-Earths resided in gas-depleted transition discs. The little time that is left for the planet to accrete remnant gas during the disk dispersal phase is limiting in this scenario, rather than the reduced gas density. In such a disc, super-Earths could form from a bunch of planetary embryos, similar to well-established models for the inner solar system (Dawson et al. 2016). Other works have continued to search for a thermodynamic explanation to prolong the cooling time of super-Earth planets (Batygin et al. 2016; Coleman et al. 2017). When super-Earths are on eccentric orbits ( e GLYPH<24> 0 : 1 ) Ginzburg & Sari (2017) showed that stellar tidal heating would suppress the cooling.  \nIn contrast, we favour a hydrodynamical explanation. It was shown in Ormel et al. (2015b) (hereafter referred to as Paper II) that envelopes are not isolated from the disc; instead, they exchange gas with the disc. Unlike the 2D case (Ormel et al. 2015a; hereafter Paper I), where we saw closed streamlines around the planet, in 3D horseshoe and in-spiralling streamlines effectively recycle material between the disc and 'envelope' (in an open system the very definition of 'envelope' becomes arbitrary). Similar results were found by Wang et al. (2014) and Fung et al. (2015). Importantly, atmosphere recycling will arrest thermodynamical cooling when the atmosphere recycling timescale (the time it takes to replenish the planets' atmosphere; Paper II) is shorter than the Kelvin-Holmholtz contraction time-scale. In Paper II we argued that recycling should be particularly effective for super-Earths: they have rather modest gravitational po-  \ntentials (because of their low ratio of physical to Hill radii) and due to their large orbital speeds they encounter much disc gas per unit time.  \nHowever, the simulations of Paper II were conducted for a Bondi radius-to-gas scale-height of 1%, 1 which is representative of a Mars-sized embryo at 1 au, not for superEarth planets. In addition, the Paper II simulations considered an isothermal equation of state and, because of the computationally-intensive nature of the simulations, could only run the simulation for several dynamical times. Here, to test whether the recycling hypothesis is truly applicable for super-Earths, we conduct simulations that (i) involve super-Earth planets (up to 5 M GLYPH<8> cores at 0 : 1 au); (ii) utilize proper radiation transport to consider the cooling of their atmospheres; and (iii) run the simulations for much longer times.  \nWith these improved physics, our key goal is to test the recycling hypothesis for super-Earth planets. We do observe a very similar flow topology as in Paper II, with disc streamlines penetrating deep in the Bondi sphere. However, unlike Paper II, we also see a low-entropy'inner core'emerging where streamlines are more circular. Indeed, the recycling hypothesis can best be tested by following the entropy. Radiative cooling decreases entropy and in every 1D thermal evolutionary model entropy always decreases with time (since the planet is hydrodynamically isolated). However, recycling has the opposite effect: high entropy gas from the disc enters the envelope, whereas low-entropy gas flows back to the disc. To quantitatively test the influence of recycling we conduct a direct comparison of radiative hydrodynamical simulations between 1D and 3D geometries. Although the 3D simulations can only run up to GLYPH<24> 10 2 -10 3 orbits, this comparison unambiguously shows that recycling delays the thermal cooling of super-Earths' atmospheres. On the other hand, the evolution time is usually not sufficient to evolve simulations into a steady state ( t 4 recycling time of the interior atmosphere), where recycling would finally arrest radiative cooling.  \nThe paper is organized as follows: In Section 2 we discuss the physical setup of our model, highlighting differences to Paper I and Paper II. Initial conditions are explained in Section 3. In Section 4 we describe the numerical configuration. Results are presented in Section 5 and discussed in Section 6. To conclude, a summary is made and followed by a brief outlook in Section 7.", '2 THEORY AND NUMERICS': "In this section we first discuss the equations of RadiationHydrodynamics in the Flux-Limited-Diffusion (FLD) approximation for radiation transfer (Levermore & Pomraning 1981). Secondly, we explain the forces that operate in the  \n1 Throughout this work, we define the Bondi radius as R Bondi = GM GLYPH<157> c 2 s , where c s is the (isothermal) sound speed of the background disc gas. Also, we take the gas scale-height as H = c s GLYPH<157> GLYPH<10> where GLYPH<10> is the local orbital frequency, and the Hill radius is R Hill = r ' M p GLYPH<157> 3 M ? ' ' 1 GLYPH<157> 3 ' where r is the disc orbital radius, M p the planet mass and M ? the stellar mass. With these definitions we have the relation R 3 H = R B H 2 GLYPH<157> 3 . For the close-in super-Earth planets that we consider, usually we have R B ' R H .  \nadopted local frame co-orbiting with the planet. Then we describe the numerical methods used to solve the governing equations.", '2.1 Equations of Radiation-Hydrodynamics': "Following previous studies (Paper I; Paper II) we describe the disc gas as an inviscid 2 and compressible fluid. In this study, we drop the previously made isothermal assumption and solve for the radiative transport of energy. We employ an ideal equation of state:  \nP = ' GLYPH<13> GLYPH<0> 1 ' GLYPH<26> e gas GLYPH<17> exp ' s ' GLYPH<26> GLYPH<13> ; (1)  \nwith the dimensionless entropy s = S GLYPH<157> c V , where S is entropy, c V the heat capacity at constant volume and GLYPH<13> is the adiabatic exponent. We consider local thermal equilibrium (LTE) and use the one-temperature approximation, assuming the gas and radiation to have the same temperature T gas = T rad = T . The three basic equations that account for conservation of mass, momentum and energy for such a fluid are the continuity, Euler's and the energy equation:  \n@GLYPH<26> @ t + r GLYPH<1> ' GLYPH<26> v ' = 0 (2)  \nGLYPH<18> @ @ t + v GLYPH<1> r GLYPH<19> v = GLYPH<0> r P gas GLYPH<26> + GLYPH<213> i a i ; (3)  \n@ E @ t + r GLYPH<1> ' E v ' = GLYPH<0> r GLYPH<1> F GLYPH<0> P gas r GLYPH<1> v (4)  \nwhere GLYPH<26> is density, t time, P gas gas pressure, v the gas velocity E = E int + E rad the sum of internal and radiation energy and F the flux of radiation energy density. Accelerations due to external forces are represented by the a i . Decomposing the energy equation into a component for internal and radiation energy, we write the internal energy density as E int = c V GLYPH<26> T , and by virtue of the assumed LTE, the radiation energy density is given by E rad = aT 4 , where a is the radiation constant. Their time evolution is governed by the following equations (e.g. Kuiper et al. 2010):  \n@ E int @ t + r GLYPH<1> ' E int v ' = GLYPH<0> P gas r GLYPH<1> v + GLYPH<3> (5)  \n@ E rad @ t + r GLYPH<1> ' E rad v ' = GLYPH<0> r GLYPH<1> F GLYPH<0> GLYPH<3> ; (6)  \nwhere GLYPH<3> describes the coupling term of matter and radiation, which in LTE vanishes. The (FLD) flux of radiation energy density is then given by  \nF = GLYPH<21> c GLYPH<26> GLYPH<20> R r E rad ; (7)  \nwith the flux-limiter GLYPH<21> , the speed of light c and the Rosseland mean opacity GLYPH<20> R . For the choice of the flux-limiter, we follow Levermore & Pomraning (1981), ensuring that the flux of radiation energy is described correctly in the limiting cases of free streaming in the optically thin regime and pure diffusion in optically thick regions.  \nUsing operator splitting we solve for the transport term  \nr GLYPH<1> ' E v ' separately during the hydrodynamical step so that the equation to be solved in the radiation/FLD step reads  \n@ E @ t = GLYPH<0> r GLYPH<1> F GLYPH<0> P gas r GLYPH<1> v : (8)  \nUsing the relation  \n@ E int @ t = c V GLYPH<26> 4 aT 3 @ E rad @ t ; (9)  \nwe can rewrite Eqn. (8) as a diffusion equation  \n@ E rad @ t = GLYPH<0> r GLYPH<1> F GLYPH<0> P gas r GLYPH<1> v 1 + c V GLYPH<26> 4 aT 3 : (10)  \nFor more technical details, see e.g. Kuiper et al. (2010).  \nIn the temperature regime that we find in our simulations ( 3 glyph[lessorsimilar] log T » K … glyph[lessorsimilar] 4 ), the pressure due to radiative forces P rad = E rad GLYPH<157> 3 GLYPH<28> P gas = GLYPH<26> k B T GLYPH<157> GLYPH<22> m H is small compared to the gas pressure. Here, k B is Boltzmann's constant and m H is the mass of atomic hydrogen. In this regime, the total energy budget is dominated by the internal gas energy ( E rad GLYPH<28> E int ) and the transport of energy will be dominated by the flux of radiation energy F . We neglect forces due to radiation and do not account for ionization or dissociation of molecules, taking a constant mean molecular weight GLYPH<22> = 2 : 353 , corresponding to gas of solar metallicity. Since molecular hydrogen is the main constituent, the ratio of heat capacities or adiabatic index is taken as GLYPH<13> = 7 GLYPH<157> 5 . We note that most one-dimensional (e.g. Lee et al. 2014) and some threedimensional studies (e.g. D'Angelo & Bodenheimer 2013) include much more complex thermodynamics, for example by accounting for ionization and dissociation of molecules in their opacities and equation of state. However, since we are primarily interested in the effects of the flow structure, we chose a simple model.", '2.2 Forces in the local frame of reference': "For all our simulations we assume the planet to be on a circular orbit with a fixed separation a from the host star. We then adopt a local non-inertial frame of reference that is centred on the planet and follows its motion around the star. Given that the planet's mass for a low-mass gaseous envelope is determined by the core mass M p ' M c , and the stellar mass that we set as M ? = 1 M GLYPH<12> , the Keplerian angular velocity of the planetary core's orbit is  \nGLYPH<10> = s G ' M GLYPH<12> + M c ' a 3 ' r GM GLYPH<12> a 3 : (11)  \nWe set up our frame to be centred on the planet, with the x -direction pointing radially away from the star and the y -direction following the orbital direction of the planet's orbit. This 'shearing-box' approach is a standard tool to study effects in discs on a smaller, local scale (e.g. Goldreich & Lynden-Bell 1965; Goldreich & Tremaine 1978; Hawley et al. 1995). In Fig. 1 we show a schematic drawing of the setup. The axis of rotation is aligned with the direction of the z -axis, so that we can write the local Keplerian frequency as Ω = GLYPH<10> e z . Since our local frame is non-inertial, we have to consider accelerations  \nFigure 1. Schematic drawing of local shearing sheet approximation. a) The planet is fixed on a circular orbit with the seperation a from the star. b) We align our co-moving frame that is centred on the planet such that the x -direction is pointing radially away from the star and the y -direction follows the orbital direction of the planet. Arrows indicate the gas velocity of the shearing flow. Without a 'headwind' arising from a pressure gradient in the disc, this flow is symmetric with respect to the planet.  \n<!-- image -->  \ndue to the Coriolis force, a cor = GLYPH<0> 2 Ω GLYPH<2> v 3 and the centrifugal force a cent = GLYPH<0> Ω GLYPH<2> ' Ω GLYPH<2> r ' . In a linear approximation (for x ; y ; z GLYPH<28> a ), the latter can be put together with contributions of stellar gravity to give the tidal acceleration a tid = 3 x GLYPH<10> 2 e x GLYPH<0> z GLYPH<10> 2 e z . The vertical component of a tid gives rise to a density stratification of the disc in z and was not considered in the low-mass limit (Paper I; Paper II), because there R Bondi GLYPH<157> H GLYPH<28> 1 . Super-Earths, however, are more massive planets, so that this additional term cannot be neglected.  \nWe only consider a Keplerian disc, so we do not take into account the global pressure force, which results in a sub-Keplerian motion of the gas, seen as a 'headwind' in the planet's frame, as described in Paper I; Paper II. Lastly, we have to take care of the gravitational interaction of the disc gas with the planetary core. As in Paper II, we use an exponential smoothing term in space to modify the Newtonian force in a way that results in a force-free boundary at the core's surface:  \n˜ a 2b = GLYPH<0> GM c r 2 exp h GLYPH<0> A GLYPH<16> r 0 r GLYPH<17> p i e r ; (12)  \nwhere r is the distance to the planet center, e r is the radial unit vector and A = 10 , p = 8 are parameters that ensure a quick transition between the smoothed and pure Newtonian regime between r GLYPH<24> r 0 GLYPH<0> 2 r 0 . In this work, the smoothing length r 0 is given by the physical radius of the planetary core r c . It is calculated assuming a mean density of the core corresponding to a rocky composition, GLYPH<26> c = 5 g cm GLYPH<0> 3 . We have confirmed that slightly varying the inner radius does not substantially change the outcome of the simulations.  \nIn order to avoid artificial initial shocks, we introduce the planet's potential gently into the disc environment. We adopt the same time-dependent smoothing term as in Paper II. The injection time t inj will vary between simulations, depending on the core mass. After 2.5 t inj , the smoothing term  \nhas already reached 95% of it's limit at unity. We note that the chosen injection times are usually a few GLYPH<10> GLYPH<0> 1 and are thus much shorter than the actual time for the coagulation of the core. While for massive envelopes self-gravity of the gas is the important effect potentially leading to runaway gas accretion, we can neglect it while the atmosphere's mass is small compared to the core mass, M env GLYPH<28> M c .", '2.3 Numerical methods': 'We solve the above described hydrodynamic equations using a modified version of the Godunov-type hydrodynamical code PLUTO (Mignone et al. 2007), complemented with a module to solve for radiation transfer in the FLD approximation (Kuiper et al. 2010). We use a Roe solver (Roe 1981) to solve the Riemann problem. The FLD equation (10) is solved using an implementation of the Generalized Minimal Residual (GMRES) method (Saad & Schultz 1986) including restarts for faster convergence. For more technical details on the code, please refer to Kuiper et al. (2010).', '3 PHYSICAL INITIAL CONDITIONS': "In this section, we first describe the physical model of the protoplanetary disc that was used to determine the environmental conditions of the planet's surroundings such as the background mid-plane density, the local disc temperature and opacities. We then motivate the one-dimensional simulations that we conducted for direct comparison.", '3.1 Disc model': "Consistent with a Keplerian disc, we set the initial velocity field to the state of an unperturbed shearing flow that results from the differential rotation of the disc. In our local  \nTable 1. List of simulations. The first column gives the simulation name. The second column shows the core mass M c . The third column gives the injection time t inj for the planet's gravity and the fourth column shows the total simulation time t end . The spatial resolution for the domain » r c ; 5 H … GLYPH<2> » 0 ; GLYPH<25> GLYPH<157> 2 … GLYPH<2> » 0 ; GLYPH<25> … is given in the fifth column. Comments and references to figures are given in the sixth column. Standard resolution simulations were conducted for three different opacities.Table 2. List of parameters for different core masses M c in units of Earth masses. We also give the core radius, Bondi and Hill radius in Earth radii. The fifth column gives the dimensionless thermal mass m = R B GLYPH<157> H of the planet, which depends on orbital location and the local disc properties. In the last two columns, we give the radius at which we find the transition between transient and close to circular streamlines in the mid-plane and the envelope mass fraction GLYPH<31> env GLYPH<17> M env GLYPH<157> M c , both after t = 1000 GLYPH<10> GLYPH<0> 1 . Values of GLYPH<31> env range from results for standard opacity to VLO opacity.  \napproximation it is given by (see also Fig. 1b)  \nv 1 ' x ' = GLYPH<0> 3 2 x GLYPH<10> e y : (13)  \nIn the frame where the planet is at rest, the co-rotation orbits at x = 0 are so as well. Gas that orbits closer to the star is overtaking the planet due to faster angular motion, while the gas further out lags behind.  \nWe follow Lee et al. (2014) in adopting the minimummass extrasolar nebula (MMEN, Chiang & Laughlin 2013) to provide parameters for the unperturbed state of the disc. Taken at a = 0 : 1 au , this yields values for the midplane gas density GLYPH<26> 0 = 6 GLYPH<1> 10 GLYPH<0> 6 g cm GLYPH<0> 3 and the temperature T 1 = 1000 K (Lee et al. 2014, Eqns. (12) & (13)). The subscript 1 stands for unperturbed mid-plane quantities. The obtained density is a factor of a few larger compared to the minimum-mass solar nebula (MMSN, Weidenschilling 1977; Hayashi 1981). The exact value of the background density should not affect the qualitative outcome of our simulations, since we cover a range of optical depths by using different opacities. However, it might introduce an additional scaling factor for the density profiles that we find.  \nAssuming an initially isothermal hydrostatic structure, the vertical density stratification in absence of the planetary core is due to the vertical component of stellar gravity. For this configuration the vertical density profile is given by (e.g. Armitage 2010)  \nGLYPH<26> 1 ' z ' = GLYPH<26> 0 exp GLYPH<20> GLYPH<0> 1 2 GLYPH<16> z H GLYPH<17> 2 GLYPH<21> ; (14)  \nwhere H is the vertical pressure scale-height. Due to the isothermal nature of this configuration, the density decreases, while the entropy of the gas is increasing towards higher altitudes.  \nWe confirmed that without including the planet's gravitational potential, the code accurately maintains the initial equilibrium state for the full simulation time.", '3.2 Opacities': "As mentioned above, a dust-free environment is a viable scenario for super-Earth planets that we are focusing on. Due to the close-in location, the atmospheres easily stay hot enough ( T glyph[greaterorsimilar] 1500 K ) to effectively evaporate dust grains, rendering their contribution to the opacity negligible (see e.g. Lee & Chiang 2015, their Fig. 6). Another reason to consider a dust-free scenario is that the planetary core has already formed and incorporated most of the dust in its surroundings. Where the background disc is below the sublimation temperature of dust, coagulation of dust grains and settling is another effective mechanism for grain removal (Ormel 2014; Mordasini 2014). These mechanisms result in opacities that are reduced by a few orders of magnitude compared to typical ISM opacities. More generally, the opacity of hot, pre-planetary atmospheres will be determined by the gas and exhibit a complex dependence on temperature, pressure and metallicity (Ferguson et al. 2005; Freedman et al. 2014). Since in this work we are mainly interested in qualitative results - not, for example, in absolute contraction timescales we consider only constant Rosseland mean opacities. Three values are explored: GLYPH<20> R = 10 GLYPH<0> 2 ; 10 GLYPH<0> 3 ; 10 GLYPH<0> 4 cm 2 g GLYPH<0> 1 . Simulations with the lower opacities are tagged as low opacity ( LO ) or very low opacity ( VLO ), respectively. For the configuration RT-M1 we additionally evolved a model with artificially increased opacity of GLYPH<20> R = 10 cm 2 g GLYPH<0> 1 , which we will refer to as VHO .  \nIn the top panel of Fig. 2 we show the azimuthally averaged optical depth GLYPH<28> ' GLYPH<18> ' of our domain for radial rays as a  \nFigure 2. Top: Azimuthally averaged total optical depth of the simulation domain (out to r = 5 H ) along radial rays of different polar angles GLYPH<18> . Different colors correspond to the three constant values for the opacity ranging from highest (blue) to lowest (green). Solid lines show the state of the unperturbed disc while dashed lines correspond to the final state of simulations for a 1 M GLYPH<8> core. Bottom: Optical depth as a function of radius along the polar axis ( GLYPH<18> = 0 ) for the very low opacity model. The meaning of line color and style is the same as in the top panel.  \n<!-- image -->  \nfunction of the polar angle. Solid lines mark the state of the unperturbed disc and dashed lines show the state at the end of our simulations for an Earth-mass core. In the bottom panel we show the optical depth along the polar direction as a function of radius. Initially, the vertical direction ( GLYPH<18> = 0 ) is less optically thick but with the contraction of the atmosphere, this feature vanishes. Even for the lowest of our opacities of GLYPH<20> R = 10 GLYPH<0> 4 cm 2 g GLYPH<0> 1 , our domain and the Bondi sphere are optically thick, justifying the FLD approximation.", '3.3 One-dimensional runs': 'For a better understanding of the implications of the 3D nature of the flow for the thermal evolution of the atmosphere, we also perform 1D spherically symmetric simulations using the same code. This allows us to look for the key differences that arise only from the assumption of spherical symmetry, instead of comparing our findings to other 1D studies that  \nusually cover a different time domain and different physical effects. In these simulations we cover the same range in radius and initialize the simulation with constant density, which is given by the mid-plane density, GLYPH<26> = GLYPH<26> 0 . We have also removed all non-radial forces from the setup and initialize the simulation with zero velocity. For this 1D configuration, we expect an atmosphere that is steadily cooling and contracting, because a simultaneous in- and outflow of gas is not possible.', '4.1 Spherical Grid, Domain and Boundary conditions': "In Paper I we showed that already in two dimensions, a polar grid outperforms a Cartesian grid in terms of accuracy and computational efficiency. For the three-dimensional problem at hand, we again adopt a spherical grid with coordinates ' r ; GLYPH<18>; GLYPH<30> ' , centred on the planet, following Paper II. Not only is the geometry more natural to the expected structure of the flow around the planetary core, by using a linear spacing of log r in the radial dimension, it allows us to sample the region of the atmosphere achieving high resolution while still covering a large domain in the radial dimension. Starting at the radius of the core our local frame extends out to five pressure scale-heights r 2 » r c ; 5 H … . At the planet surface, we chose a reflective boundary-condition (BC), so there is no transport of matter or radiation. For the outer radial boundary, we set ghost-cells to the unperturbed shearing-flow state of the disc, fixing GLYPH<26> = GLYPH<26> 1 ' z ' ; T = T 1 and v = v 1 ' x ' . Parameters such as background density and temperature depend on the disc model that is discussed above and the velocity field is given by the shearing flow in a Keplerian disc (Eqn. 13).  \nOne disadvantage of this grid is the singularity at the pole ' GLYPH<18> ! 0 ' : moving towards the pole, the d GLYPH<30> length of grid cells quickly decreases, forcing a very small time-step according to the Courant condition. In Section 2.3 of Paper II, we used a linear spacing of cos GLYPH<18> . This partially relaxes the constraint on the time-step by effectively increasing the size of the cells that lie at high latitudes. While this method is feasible for the study of low-mass planets, we realized that it is not suited for more massive cores. As the vertical disc stratification and vertical component of the tidal force become relevant, the coarse polar sampling of the first grid cells quickly resulted in visible grid artefacts during test runs. Even in test runs with doubled resolution in GLYPH<18> , these effects were clearly visible. Thus, we decided to trade a longer runtime of the simulation for higher accuracy and chose a linear sampling of the polar angle GLYPH<18> .  \nSince we only consider a Keplerian disc, the equations that describe this system have a symmetry under a rotation of GLYPH<25> in azimuth ( GLYPH<30> ! GLYPH<30> + GLYPH<25> ), which in Cartesian coordinates corresponds to ' x ; y ; z ' ! 'GLYPH<0> x ; GLYPH<0> y ; z ' . In spherical coordinates it can be expressed as A ' r ; GLYPH<18>; GLYPH<30> ' = A ' r ; GLYPH<18>; GLYPH<30> + GLYPH<25> ' for any scalar or vector quantity A . This allows us to cover the azimuthal dimension only in the forward direction of the planet GLYPH<30> 2 » 0 ; GLYPH<25> … , adopting periodic BCs. The system is also symmetric with respect to the mid-plane ' z = 0 ' , so that we only cover the upper hemisphere, GLYPH<18> 2 » 0 ; GLYPH<25> GLYPH<157> 2 … , that has z > 0 and take the boundary condition in the mid-plane ' GLYPH<18> = GLYPH<25> GLYPH<157> 2 '  \nto be reflective, ensuring that no hydrodynamic or radiative fluxes cross this boundary. At the polar singularity, we use the so-called GLYPH<25> -boundary-condition to address ghost-cells that effectively have negative GLYPH<18> : A ' r ; GLYPH<0> GLYPH<18>; GLYPH<30> ' = A ' r ; GLYPH<18>; GLYPH<30> + GLYPH<25> ' for any quantity A .", '4.2 Resolution and simulation time': "Compared to previous isothermal simulations, we follow the evolution of the system roughly a factor one hundred times longer than the isothermal simulations of Paper II in order to accurately quantify the cooling and recycling behaviour of the planet atmosphere. For this reason, we conducted the 3D simulations in our standard resolution, given by N r GLYPH<2> N GLYPH<18> GLYPH<2> N GLYPH<30> = 128 GLYPH<2> 32 GLYPH<2> 64 cells in the domain » r c ; 5 H … GLYPH<2> » 0 ; GLYPH<25> GLYPH<157> 2 … GLYPH<2> » 0 ; GLYPH<25> … for a total simulation time of t end = 1000 GLYPH<10> GLYPH<0> 1 ' 5 yr . One exception to the rule is simulation RT-M1-VLO , which we evolved up to t = 2000 GLYPH<10> GLYPH<0> 1 in order to judge the convergence to a steady state. As an example for an Earth-mass core, our standard resolution allows for 67 cells inside the Bondi radius in the radial direction. In a resolution study (see Appendix A), we conducted simulations with half/twice as many cells ( loRs / hiRs ) in each dimension in order to find and understand resolution effects. However, hiRs runs were computationally too costly to cover the same time domain. A list of simulations can be found in Table 1.", '5 RESULTS': 'In this section we present and analyse the results of the radiative hydrodynamical models. First we discuss the general structure of the atmospheric region and dynamical interaction with flows from the disc in Section 5.1. We then describe the thermal properties, focusing on the transport of entropy by advection and radiation in Section 5.2. A direct comparison between 1D and 3D simulations is made in Section 5.3.', '5.1 Three-dimensional flow structure': "For all our simulations we find that the 'atmosphere' - material within the Bondi and Hill sphere - is connected to the disc by flows both entering and leaving: recycling always operates. In Fig. 3 we show the time evolution of the flow and density structure in a mid-plane section for a two Earth mass core as a typical example. The blue and green dashed circles mark the Bondi and Hill radius, respectively. At j x j GLYPH<24> R B we recognize the shearing flow, streamlines further in are bent and pass closer to the core. We also recognize typical horseshoe streamlines in the co-rotation region around x = 0 which are slightly skewed and not symmetric with respect to the y -axis.  \nThe resulting mass fluxes GLYPH<219> M in GLYPH<157> out ' r ' through a spherical surface at radius r and cumulative mass M < r inside this sphere are plotted in Fig. 4 for two simulations. 4 The top panel shows results for RT-M1-VLO , the bottom panel for RTM2-VLO . Mass fluxes are higher in the outer parts of the  \nBondi sphere and drop off towards the center. The rates of in- and outflow match very closely, resulting in a net accretion that is orders of magnitudes smaller. By taking the ratio of mass fluxes and the cumulative mass, we can define a recycling time-scale  \nt recycling GLYPH<17> M < r GLYPH<219> M in ; (15)  \nsimilar to that in Paper II. However, here we do not sort the mass fluxes by the minimum radius they reach by integrating their streamlines, but look at the local mass flux. In other words, this does not tell us directly, how deep the streams penetrate the atmosphere but still gives us an understanding of how large the mass exchange is at each radius. This quantity is shown in Fig. 5 for RT-M1-VLO at different times as a function of radius. For r glyph[greaterorsimilar] 0 : 4 R B , we have t replenish GLYPH<24> 1 GLYPH<0> 10 GLYPH<10> GLYPH<0> 1 i.e., days, while in the inner regions we usually find t recycling GLYPH<24> 100 GLYPH<10> GLYPH<0> 1 i.e., years. With time, the recycling time in the outer part slightly increases, while we see a decrease close in. The latter trend has to be taken with care however, since due to our local definition of the recycling time, circulation of material that has some radial component will affect this quantity, without effectively transporting in gas from the disc. Thus, recycling times might be underestimated for layers close to the core. The fraction of a Bondi radius at which this transition in recycling time happens is smaller for a two Earth mass core compared to a one Earth mass. This means that the layers of the envelope are not isolated but steadily exchange gas on a larger scale than the net accretion that is allowed for by radiative cooling.  \nOver time, due to the atmosphere's contraction, mass accumulates close to the core. As angular momentum is transported to its inner parts, rotation increases. The retrograde nature of this rotation is different to the results of Paper II, where we found outward spiralling streams with a prograde sense in the mid-plane.  \nTo assess the importance of rotation we show the shell averaged density profile and azimuthally averaged azimuthal velocity in Fig. 6. The top panel compares the simulation results (blue) of 3D (solid) and 1D (dashed) models with analytical hydrostatic solutions corresponding to the planet's smoothed gravity, adopting an isothermal (black dashed) and adiabatic (red dash-dotted) equation of state. The gravitational smoothing which causes the vanishing density gradient at the inner boundary, hardly affects the atmosphere structure. Overall, the simulation results agree with an adiabatic structure with deviations close in, where the atmosphere has contracted the most. We confirm that right after the injection of the planet, before the atmosphere has a chance to cool, results of RT runs agree closely with adiabatic profiles. For higher opacity runs, this remains true for longer times than VLO runs. Our very high opacity test VHO maintains virtually constant adiabatic structure up to t = 1000 GLYPH<10> GLYPH<0> 1 .  \nThe bottom panel shows the azimuthal velocity (blue) with positive/negative values drawn as a solid/dashed line. Comparing the simulation results to velocities of Keplerian rotation, we find that the any of the rotation seen in Fig. 7 is sub-Keplerian and dynamically unimportant in supporting the atmosphere. Clearly, there is no evidence for a circumplanetary disc that is often found in simulations of very massive planets (e.g. Ayliffe & Bate 2012; Tanigawa et al.  \nFigure 3. Time evolution of the density and the flow structure in the mid-plane for simulation RT-M2-VLO . Dashed circles mark the Bondi (blue) and Hill (green) radius. The arrows show gas velocity and the grey lines are corresponding streamlines. By symmetry the vertical velocity component vanishes here.  \n<!-- image -->  \nFigure 4. Time series of cumulative gas mass M < r (solid blue) inside radius r together with accretion rates corresponding to in- (dashed) and outflow (dotted) scaled with the Keplerian frequency. Top: RT-M1-VLO , bottom: RT-M2-VLO . While mass fluxes are high in the outer regions of the Bondi sphere, they decreased closer in. Incoming and outgoing streams compensate each other, resulting in much smaller net accretion. Over time the inner regions allow for more exchange of gas. The pileup of gas close to the planet due to cooling and contraction is visible as an increase in M < r close in.  \n<!-- image -->  \nFigure 5. Time evolution of the recycling time t recycling (Eqn. 15) for RT-M1-VLO as a function of radius. While gas in the outer regions is quickly replenished, the time-scale on which recycling occurs closer in is usually larger, as the atmosphere becomes denser and mass fluxes reduce (compare Fig. 4).  \n<!-- image -->  \n2012; Szul'agyi et al. 2016), at this point. The sign change in azimuthal velocity from positive (prograde) close to the planet surface to negative (retrograde) at the point, where gravitational smoothing kicks in is most likely caused by this simplified treatment of gravity, since varying the smoothing parameter r smooth slightly changed its location (see Appendix A).  \nIn Table 2, we give the radius r circ , where we find the transition between transient and circular, (quasi)bound streamlines in the mid-plane as a fraction of the Bondi radius. It is defined as the largest radius inside the atmosphere, where the azimuthally averaged ratio of radial to azimuthal velocities falls below a threshold of ten percent, so where hj v r GLYPH<157> v GLYPH<30> ji GLYPH<30> GLYPH<20> 0 : 1 . Although the choice of this threshold value is somewhat arbitrary, it helps in identifying a trend for the dependence of this radius on the core mass. We emphasize that this crude analysis is limited to the mid-plane and does not directly tell us whether there is still massexchange happening at higher latitudes, which we will address in Section 5.2. We also give the envelope mass fraction GLYPH<31> env GLYPH<17> M env GLYPH<157> M c at t = t end (i.e. the mass inside R B ; H ) for different core masses for standard and VLO opacity. While the circulation radius decreases with increasing core mass, the envelope mass fraction increases and is typically of several percent for super-Earth mass cores. The envelope masses show little dependence on the adopted opacity at this point. Since we consider a gas-rich environment, even fully adiabatic atmospheres have significant mass fractions of some percent.  \nWe note however, that generally our simulations have not yet reached steady state. Thus, by comparing envelopes for different core masses at a fixed time we are viewing them at different stages of their thermal evolution.", '5.2 Entropy transport': "The implications of the previously-discussed recycling on the cooling of the planet's atmosphere can best be described in  \nFigure 6. Simulation results for RT-M1-VLO at t = 2000 GLYPH<10> GLYPH<0> 1 . Top: shell-averaged density profile (blue) for 3D (solid) and 1D (dotted) runs compared to analytical hydrostatic solutions employing an isothermal (black dashed) and adiabatic (red dash-dotted) EoS. The simulation results closely follows the adiabatic profile, deviating only close in. Bottom: azimuthally averaged azimuthal velocity in the mid-plane (blue) split into positive (solid) and negative (dashed) component. We also show rotational profiles of Keplerian rotation in a Newtonian (grey dashed) and in our smoothed potential (grey dotted). We recognize that the retrograde rotation from Fig. 7 is well sub-Keplerian.  \n<!-- image -->  \nterms of entropy. In the following we express the difference in specific entropy in units of k B per baryon as  \nGLYPH<1> S k B GLYPH<157> baryon = 1 ' GLYPH<13> GLYPH<0> 1 ' GLYPH<22> ln P GLYPH<157> P 1 ' GLYPH<26> GLYPH<157> GLYPH<26> 1 ' GLYPH<13> ; (16)  \nwhere k B is Boltzmann's constant and again the subscribed quantities are those of the unperturbed mid-plane. In Fig. 7 we show the mid-plane flow for RT-M1-VHO and RT-M1-VLO , color coded by GLYPH<1> S . The two bottom rows show the same data with different ranges of the colormap. A clear distinction between two regions can be made: (i) an outer nearisentropic region where S is the same as that of the disc; and (ii) an inner core where S is clearly lower. In the outer regions, because of continuous exchange of gas from the disc, most of the gas is able to maintain entropy close to the initial value. This picture remains almost unchanged between the snapshots at t = 500 GLYPH<10> GLYPH<0> 1 and t = 1000 GLYPH<10> GLYPH<0> 1 , suggesting rapid convergence to steady state. However, the evolution proceeds differently within the hot, close-in part of the atmosphere, where gas cools efficiently by radiation. Initially, the extent of the dark area of low entropy gas grows. With  \nFigure 7. Mid-plane flow structure for RT-M1-VHO (top) and RT-M1-VLO (bottom), similar to Fig. 3. Dashed circles mark the Bondi (blue) and Hill (red) radius. Coded with colors is diffrence in specific entropy to the unperturbed mid-plane value. The two bottom rows show the same data, but with different scaling of the colormap. The bottom row emphasizes the mixing of gas and resulting lower entropy outflow in the mid-plane.  \n<!-- image -->  \ntime however, the inner region is also affected by shearing and horseshoe orbits, which inject high entropy gas while removing low entropy gas from the very inner regions. This process is emphasized in the bottom panel of Fig 7. Comparing the results for different opacities, we see that more efficient cooling (and thus lower entropy) promotes the transition towards more circular - and arguably more bound motion around the core. We find a similar structure for simulation RT-M2-VLO .  \nWe emphasize however that 3D effects are key. Fig. 8 displays the mass flux F m (in units of Earth masses per inverse orbital frequency per surface of the Bondi sphere) and the specific entropy GLYPH<1> S on a spherical surface at r =  \n1 R B ; 0 : 5 R B and 0 : 35 R B . All quantities are averaged over a time of one GLYPH<10> GLYPH<0> 1 . Due to the symmetry of our model, we describe the longitude GLYPH<30> of only one feature, the other will be found at GLYPH<30> 0 = GLYPH<30> + GLYPH<25> . The longitude of GLYPH<30> = 3 GLYPH<25> GLYPH<157> 2 corresponds to looking at the planet from the GLYPH<0> y direction, following its orbital direction. On the surface of the Bondi sphere the mass flux through the atmosphere is dominant in the mid-plane, where mainly the shearing and horseshoe flow is passing perpendicularly through it. The intensity of the mass flux is highest there and reduces towards higher latitudes, where the outward flow dominates. Comparing the flow pattern to the entropy plot below, we recognize the lowS outflow from Fig. 8 in the mid-plane. For this feature at GLYPH<30> ' 3 GLYPH<25> GLYPH<157> 4 , the  \nFigure 8. Mass flux F m and specific entropy GLYPH<1> S on spherical surfaces at r = R B (top), r = 0 : 5 R B (middle) and r = 0 : 35 R B (bottom) for RT-M1-VLO at t = 1000 GLYPH<10> GLYPH<0> 1 . All quantities are averaged over 1 GLYPH<10> GLYPH<0> 1 . The mass flux has similar units as accretion rates above, but is additionally scaled with the surface of the Bondi sphere. It is defined such that positive values (blue) correspond to inward mass flux. We only show the upper hemisphere because of the mid-plane symmetry. The longitude of GLYPH<30> = 3 GLYPH<25> GLYPH<157> 2 corresponds to the view that is following the planets orbital direction. On the surface of the Bondi sphere, we still recognize the vertical stratification in entropy giving GLYPH<1> S > 0 at high altitudes. There, we also find a clear trend for leaving gas (red) having lower entropy compared to gas entering the atmosphere (blue).  \n<!-- image -->", '12 N. P. Cimerman et al.': "lowest entropy flows are close to the mid-plane. However, it extends vertically around this longitude, so that even at high latitudes the out-flowing streams show lower entropy compared to their unperturbed state. In general we find that on the Bondi sphere, regions of accretion overlap with regions of higher entropy compared to areas of outflow. Making the same surface cut at r = 0 : 5 R B , we find a similar structure, with the main accretion streams having GLYPH<1> S GLYPH<24> 0 , the disc's mid-plane entropy. At this distance from the core mass fluxes are still the highest close to the mid-plane, while at higher latitudes they virtually vanish. Closer in at r = 0 : 35 R B , the mass flux in the mid-plane almost vanishes, corresponding to the circular streamlines that we find there, allowing only for very little outflow. However, the inner region is still exchanging gas at higher latitudes. Gas is leaking the strongest at GLYPH<30> = GLYPH<25> GLYPH<157> 4 , which corresponds to the line of x = y , regions of highest accretion are shifted by GLYPH<25> GLYPH<157> 2 . In contrast to the outer regions, at r = 0 : 35 R B , we don't find any gas that has GLYPH<1> S ' 0 and the correlation between the direction of the flow and its entropy has blurred. Together with the reduced mass flux (Fig. 4), this suggests that advective entropy transport is less efficient by this radius.", '5.3 Comparison 3D vs. 1D': "We next illustrate the effects of recycling on the entropy with a direct comparison to 1D models. In Fig. 9 we show the spherically averaged radial entropy profile for an Earthmass core, directly comparing 1D (solid) and 3D (dashed) in the top row and showing their difference GLYPH<1> 3D ; 1D ' GLYPH<1> S ' in the bottom row. Results are given for different opacities (columns) at different times (colours). Starting close to the unperturbed mid-plane entropy, the innermost regions that cool the most due to radiation lose entropy, which happens most efficiently for the lowest opacity. Initially, 1D and 3D simulations closely follow each other. 5 However, after t = 500 GLYPH<10> GLYPH<0> 1 , we already see a difference emerging for the VLO runs. While the 1D simulation continues to lose entropy by radiative cooling, deviating from the initial value all the way out to r > R Bondi , the 3D atmosphere effectively maintains its initial entropy in regions r glyph[greaterorsimilar] 0 : 4 R Bondi due to the combined effect of advection and radiative cooling, showing that most of the gas in the Bondi sphere is not losing entropy as efficiently as in the 1D case. In regions where the entropy is higher in 3D compared to 1D ( r ' 0 : 2 GLYPH<0> 1 R B ), the density is lower as can be seen in Fig. 6. This nicely illustrates how high entropy counteracts contraction. We find similar results the for higher mass case RT-M2-VLO .  \nWith time this discrepancy between 1D and 3D atmospheres grows mainly due to the fact that the evolution of the 3D atmosphere slows down, suggesting convergence to a quasi-steady state that is close to adiabatic in most of the Bondi sphere. Clearly, the effects of the low-opacity run are most pronounced, because the radiative diffusion timescale scales inversely with opacity. But for the intermediate-  \nopacity ( LO ) run we also see the first signs of a mismatch emerging after t = 10 3 GLYPH<10> GLYPH<0> 1 . In fact, the results show selfsimilarity: the GLYPH<20> R = 10 GLYPH<0> 3 and 10 GLYPH<0> 4 cm 2 g GLYPH<0> 1 runs in panel (a) match those of GLYPH<20> R = 10 GLYPH<0> 2 and 10 GLYPH<0> 3 cm 2 g GLYPH<0> 1 in panel (c). Also note that at 0.1 au, 10 3 GLYPH<10> GLYPH<0> 1 only corresponds to a physical time of 5 yr ! We therefore expect that the behaviour seen in the VLO run will (i) continue to diverge from the 1D evolution; and (ii) will be representative for the higher opacities.  \nTo illustrate the convergence towards a steady state for simulation RT-M1-VLO we have evolved this simulation for twice as long, up to t = 2000 GLYPH<10> GLYPH<0> 1 . The results are shown as the cyan lines in the right panel of Fig. 9. Clearly, dashed green and cyan lines follow each other closely with deviations only in the deepest layers of the atmospheres, which are not as effectively replenished.", '6.1 Recycling for radiative atmospheres': "Similar to Paper II, we found that the Bondi and Hill sphere are dynamically interacting with the parent protoplanetary disc by constantly exchanging gas. Accounting for radiative transport processes, we found a flow structure similar to the isothermal 3D atmospheres. Most importantly, there is a continuous supply of high entropy gas from the disc that is brought deep into the atmosphere along horseshoe and shearing orbits. Thus, we confirm that such atmospheres cannot be described as isolated. For the higher mass planets considered in this work, a larger fraction of the Bondi sphere is dynamically dominated by the shearing flow (compare Fig. 3 to Paper II, Fig. 2). Different from Paper II, however, we do not find out-spiralling streamlines in the mid-plane, but circular streamlines close to the planet (see Sect. 6.2).  \nRecently, Wang et al. (2014) and Fung et al. (2015) also reported that isothermal, embedded atmospheres represent open systems. Studying planets with lower thermal masses ( R Bondi GLYPH<157> H GLYPH<24> 0 : 1 ), Wang et al. (2014) found circumplanetary discs forming. Their larger orbital distance at 5 : 2 au and the assumed isothermality are most likely beneficial for the formation of such a disc, since we do not find this feature for planets with even higher thermal masses. In a global simulation using refined grids towards the planet and also solving for radiation transport, D'Angelo & Bodenheimer (2013) identified bound envelopes of sizes GLYPH<24> 0 : 4 GLYPH<0> 0 : 9 R B , far larger than our case. Their parameter space is still considerably different from ours as they used a global domain, explicit gas viscosity, non-ideal effects (hydrogen dissociation and ionization) and planetesimal accretion. Therefore, it is still unclear which mechanism is responsible for isolating the preplanetary atmosphere. But it is not radiation transport.", '6.2 Emergence of inner, low-entropy core': "Although we cannot identify a strictly bound region in any of our simulations, we do find that a low entropy core emerges around the planet and that this region is associated with more circular streamlines (Fig. 7) and longer replenishment time-scales (Fig. 4). This somewhat discrete picture (low S , more bound, region interior to r = 0 : 2 R Bondi ; highS , thoroughly recycled region outside) nonetheless differs from  \nFigure 9. Time evolution of the radial profile of the spherically averaged specific entropy for a 1 M GLYPH<8> core. Colours mark different times while the three columns display results for decreasing opacity from left to right. The top row shows both 1D (solid) and 3D (dashed) results while the bottom row shows their difference. As expected, the radiative cooling is fastest for the lowest optical depth. Quickly, the entropy in the 1D atmosphere deviates from the 3D result, which effectively remains at the initial entropy for r glyph[greaterorsimilar] 0 : 4 R B . This discrepancy grows as the cooling stalls in the 3D case, while the 1D atmosphere continues to lose entropy. To illustrate the convergence to a steady state for the 3D model RT-M1-VLO , the cyan lines in the right panel show the state at twice our usual simulation time t = 2 t end . Note that the gradient inversion close to r = 0 is an artefact of the smoothed gravity.  \n<!-- image -->  \nPaper II. Apart from the more realistic thermodynamics, it should also be realized that in this work we focused on the long-term behaviour, while sacrificing resolution (which mattered according to Paper II to be important). In addition, the adopted gravitational smoothing in our model renders the study of the innermost regions somewhat problematic; a pure Newtonian potential is far preferable (but runs into problems with PLUTO , see Paper II). Runs at higher resolution (Fig. A1) and different values for r smooth (Fig. A3) do show a slightly different flow pattern for the inner core region, mainly affecting the rotational feature, but also that those parameters do not affect the main results. This is discussed further in the Appendix A.  \nA relevant question is how this inner lowS region will evolve on evolutionary time-scales (recall again that the GLYPH<24> 10 3 orbits only corresponds to 5 yr). From Fig. 7 one observes that the lowS region expands with time. Possibly, this may imply that the core becomes more isolated, although Fig. 4 illustrates that radial recycling rates remain robust. 6 Likewise, the low-opacity run in Fig. 9 also shows that the point where the 3D entropy profile detaches from the 1D curve moves inwards with time: r = 0 : 19 at t = 500 GLYPH<10> GLYPH<0> 1 to r = 0 : 14 at t = 10 3 GLYPH<10> GLYPH<0> 1 . These are clear manifestations of continued recycling. Finally, physical effects, like dissociation of H 2 and evaporation of silicate grains or pebbles, may decide the evolution of the inner core (possibly already during the assembly of the super-Earth planet). From the above results  \nit is clear that the efficiency of recycling strongly depends on the fraction of atmospheric gas that is replenished with gas from the disc. Given the above described topology of the flow, the distribution of mass as a function of distance from the core plays an important role. From Fig. 4, we see that in our simulations, where we fixed GLYPH<13> = 7 GLYPH<157> 5 for simplicity, the mass of the atmosphere is dominated by the outer parts that show significant recycling. However, Piso & Youdin (2014) and Lee et al. (2014) argue that close to the planet, where the temperature is sufficiently high ( T GLYPH<21> 2500 K ), accounting for dissocation of H 2 will result in more realistic values of GLYPH<13> GLYPH<20> 4 GLYPH<157> 3 . Depending on the ratio of the core radius to the outer radius of the atmosphere r c GLYPH<157> R B ; H , this value of GLYPH<13> can represent a transition to an atmospheric structure, where most of the mass accumlates close to the core. If this is the case, recycling time-scales deep in the atmosphere could dramatically increase. However, by integration of the exact hydrostatic density profile (not the power-law approximation) we find that for values of r c GLYPH<157> R B ; H (see Table 2) corresponding to our planets, this transition occurs at an adiabatic index of only GLYPH<13> ' 1 : 2 . Still, non-ideal effects as H-dissociation (resulting in non-uniform GLYPH<13> ) could yet render the atmosphere centrally condensed, but this we leave as an exercise to be investigated by future models. For a more realistic treatment, future studies should also consider tabulated, non-constant opacities. Altogether, given the short physical times covered, it is hard to predict the long-term thermal and dynamical evolution of this critical region based on the simulations of this work.", '6.3 Implication: Survival of Super-Earths planets': "But these considerations only apply to the inner core; in most of the super-Earth's atmosphere recycling will be robust, consistent with the results of Paper II. In Paper II, we showed that the time-scale for atmospheric recycling is a steep function of disk orbital radius a : t recycling GLYPH<24> a 2 : 75 . 7 Physically, the strong dependence on orbital radius follows from the fact that the planet encounters more gas at a higher rate in the inner disk. Recycling will hence be more efficient for close-in orbits, where t recycling GLYPH<28> t cool , and provides in particular for super-Earth planets an avenue to prevent the premature collapse of their envelopes. With the radiationtransport simulations of this work, the recycling hypothesis, first suggested by Paper II, has been put on a firmer ground.  \nIt is encouraging that our simulations readily yield envelopes with mass fractions of a few percent, consistent with current predictions from observations (e.g. Wolfgang & Lopez 2015). On time-scales of a disc lifetime, we expect these values to be subject to change. Further growth or atmospheric loss due to above mentioned mechanisms will affect the envelope. Here, we just argue that recycling does not restrict atmosphere masses to values that are too low, given that super-Earths formed in a gas-rich environment.  \nThe maximum mass of the planets that we tested has been fixed at 5 M GLYPH<8> . At this mass the ratio of the Bondi radius-to-gas scale-height, R Bondi GLYPH<157> H = 1 : 9 , which, being larger than unity, suggests that non-linear effects become important. Because the thermal mass reduces with higher temperature discs, we expect that 10 M GLYPH<8> planets in hotter disks will be consistent with the local result of this paper. On the other hand, for a thermal mass much beyond unity we expect that the gravitational backreaction of the disc on the planet atmosphere matters; effects such as shocks at the spiral arms and gap opening are likely to affect the evolution of gaseous envelopes (e.g. Kley & Nelson 2012). While the high resolution around the planet is key to study atmospheric structure, global models will be needed to properly capture the response of the protoplanetary disc in order to study long-term evolution. Approaches such as nested grids (e.g. D'Angelo & Bodenheimer 2013) or a coupling of local and global simulations might present a way to achieve this goal in the future for planets in this regime.", '7 SUMMARY': "In this study, we have used hydrodynamical simulations including radiative transport and directly account for the thermal effects of atmospheric recycling. Adopting the same spherical geometry in a local shear frame as in Paper I & Paper II, we studied the evolution and structure of the emerging gaseous envelope for embedded planets with higher thermal masses, approaching m GLYPH<24> 1 . In particular, we focused on low-mass planets on close-in orbits at 0.1 au, close to where many super-Earths are found today and added the vertical component of stellar gravity and the resulting disc stratification to our setup.  \nOur key findings are:  \n- (i) Allowing for radiation transport, we find that preplanetary atmospheres represent open systems that show significant rates of recycling, qualitatively and quantitatively (Fig. 4) similar to the isothermal runs for very low mass planets of Paper II,with the important exception of the more isolated inner core that we find here.\n- (ii) Recycling equilibrates the entropy between disk and envelope: lows gas is removed from the envelope, whereas highs enters the atmosphere, in particular from high latitudes.\n- (iii) Through a comparison of the entropy profile between 1D and 3D geometries, we found that recycling suppresses cooling of the envelope. Already after a few years the 3D flow simulations show significant departure from 1D (hydrodynamically-isolated) atmospheres.\n- (iv) Contrasting the results of isothermal 2D simulations (Paper I) but in agreement with isothermal 3D simulations (Paper II), rotation does not play a significant role in the emerging atmospheres, which are fully pressure supported.  \nAccounting for more realistic thermodynamics via radiation transport, we have shown that recycling constitutes an important mechanism to change the thermal evolution of low-mass planets - super-Earths of thermal mass m glyph[lessorsimilar] 1 ( M p glyph[lessorsimilar] G GLYPH<157> c 2 s H ' - compared to isolated 1D atmospheres. Recycling is therefore an attractive way for super-Earths to survive in gas-rich discs. Staying in the sub-critical regime they can avoid runaway gas accretion and becoming hot Jupiters even when the atmosphere opacities are low.", 'ACKNOWLEDGEMENTS': 'We thank many colleagues for helpful discussions that have benefited this work. Also, we would like to thank the anonymous referee for thought- and helpful comments that have substantially improved the quality of this manuscript. NPC and RK acknowledge financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1. CWO is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum fur Datenverarbeitung of the University of Tubingen, the state of Baden-Wurttemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG.', 'REFERENCES': "Armitage P., 2010, Astrophysics of Planet Formation. Cambridge University Press  \n- Ayliffe B. A., Bate M. R., 2012, MNRAS, 427, 2597 Batygin K., Bodenheimer P. H., Laughlin G. P., 2016, ApJ, 829, 114\n- Cameron A. G. W., 1978, Moon and Planets, 18, 5\n- Chiang E., Laughlin G., 2013, MNRAS, 431, 3444\n- Coleman G. A. L., Papaloizou J. C. B., Nelson R. P., 2017, preprint, ( arXiv:1705.08147 )  \nD'Angelo G., Bodenheimer P., 2013, ApJ, 778, 77  \nDawson R. I., Lee E. J., Chiang E., 2016, ApJ, 822, 54  \nFerguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik  \n- J. G., Hauschildt P. H., Heffner-Wong A., Tamanai A., 2005, ApJ, 623, 585  \nFreedman R. S., Lustig-Yaeger J., Fortney J. J., Lupu R. E.,  \nMarley M. S., Lodders K., 2014, ApJS, 214, 25  \nFressin F., et al., 2013, ApJ, 766, 81  \nFulton B. J., et al., 2017, ArXiv e-prints:1703.10375,  \nFung J., Artymowicz P., Wu Y., 2015, ApJ, 811, 101  \nGinzburg S., Sari R., 2017, MNRAS, 464, 3937  \nGoldreich P., Lynden-Bell D., 1965, MNRAS, 130, 97  \nGoldreich P., Tremaine S., 1978, ApJ, 222, 850  \nHadden S., Lithwick Y., 2014, ApJ, 787, 80  \nHawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742  \nHayashi C., 1981, Progress of Theoretical Physics Supplement,  \n70, 35  \nHori Y., Ikoma M., 2010, ApJ, 714, 1343  \nHubickyj O., Bodenheimer P., Lissauer J. J., 2005, Icarus, 179,  \n415  \nIkoma M., Emori H., Nakazawa K., 2001, ApJ, 553, 999  \nJin S., Mordasini C., 2017, preprint,  \n(  \narXiv:1706.00251  \n)  \nJontof-Hutter D., et al., 2016, ApJ, 820, 39  \nKley W., Nelson R. P., 2012, ARA&A, 50, 211  \nKuiper G. P., 1951, Proceedings of the National Academy of Sci-  \nence, 37, 1  \nKuiper R., Klahr H., Dullemond C., Kley W., Henning T., 2010,  \nA&A, 511, A81  \nLee E. J., Chiang E., 2015, ApJ, 811, 41  \nLee E. J., Chiang E., Ormel C. W., 2014, ApJ, 797, 95  \nLevermore C. D., Pomraning G. C., 1981, ApJ, 248, 321  \nLopez E. D., Fortney J. J., 2013, ApJ, 776, 2  \nLopez E. D., Fortney J. J., 2014, ApJ, 792, 1  \nMarcy G. W., Isaacson H., Howard A. W., Rowe J. F., Jenkins  \nJ. M., et al. 2014, ApJS, 210, 20  \nMayor M., et al., 2011, preprint,  \n(  \narXiv:1109.2497  \n)  \nMignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O.,  \nZanni C., Ferrari A., 2007, ApJS, 170, 228  \nMills S. M., Mazeh T., 2017, ApJ, 839, L8  \nMordasini C., 2014, A&A, 572, A118  \nMordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139  \n- Movshovitz N., Bodenheimer P., Podolak M., Lissauer J. J., 2010, Icarus, 209, 616  \nOrmel C. W., 2014, ApJ, 789, L18  \n- Ormel C. W., Kuiper R., Shi J.-M., 2015a, MNRAS, 446, 1026  \nOrmel C. W., Shi J.-M., Kuiper R., 2015b, MNRAS, 447, 3512  \nOwen J. E., Wu Y., 2013, ApJ, 775, 105  \nOwen J. E., Wu Y., 2016, ApJ, 817, 107  \n- Owen J. E., Wu Y., 2017, preprint, ( arXiv:1705.10810 )  \nPetigura E. A., Marcy G. W., Howard A. W., 2013, ApJ, 770, 69  \nPiso A.-M. A., Youdin A. N., 2014, ApJ, 786, 21  \nPollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush  \nT., Fong W., 1994, ApJ, 421, 615  \nRafikov R. R., 2006, ApJ, 648, 666  \nRoe P. L., 1981, Journal of Computational Physics, 43, 357  \nSaad Y., Schultz M. H., 1986, SIAM J. Sci. Stat. Comput., 7, 856  \nSafronov V. S., 1969, Evolution of the Protoplanetary Cloud and  \nFormation of Earth and the Planets. Moscow: Nauka. Transl.  \n1972 NASA Tech. F-677  \nStevenson D., 1982, Planetary and Space Science, 30, 755  \n- Szul'agyi J., Masset F., Lega E., Crida A., Morbidelli A., Guillot T., 2016, MNRAS, 460, 2853  \nTanigawa T., Ohtsuki K., Machida M. N., 2012, ApJ, 747, 47  \nVenturini J., Alibert Y., Benz W., Ikoma M., 2015, A&A, 576,  \nA114  \nWang H.-H., Bu D., Shang H., Gu P.-G., 2014, ApJ, 790, 32 Weidenschilling S. J., 1977, Ap&SS, 51, 153 Wolfgang A., Lopez E., 2015, ApJ, 806, 183  \nFigure A1. Spherically averaged radial entropy profile for RT-M1VLO run in different resolutions. Only the first two panels include the hiRs results, since this is the last available output time for that resolution.  \n<!-- image -->", 'APPENDIX A: EFFECTS OF RESOLUTION AND SMOOTHING ON INNER CORE': 'To understand the effects of resolution on the main result of our study, Fig. A1 shows the spherically averaged entropy profile that was obtained for three different resolutions. While initially all resolutions agree closely, a small deviation in the entropy profile is visible at t = 1000 GLYPH<10> between the low and standard resolution. Mainly above described the inner core is affected, as the loRs maintains slighlty higher entropy. However, the difference is small and the overall structure remains valid. Results of hiRs runs agree more closely to those of our standard resolution, suggesting convergence. Most importantly: Independent of resolution, all runs yield similar values for r circ , which means that in all runs of the resolutin study, the upper layers are affected by recycling in the same way. Similar to the top panel of Fig. 7, we show the mid-plane structure for the lower resolution run RT-M1-VLO-loRs at the end of the simulation in Fig. A2  \n<!-- image -->  \n]  \n[  \nFigure A2. Mid-plane entropy and flow structure for RT-M1-VLOloRs . Similar to Fig. 7.Figure A3. Mid-plane density and flow structure for RT-M2-VLO comparing larger (top) and smaller (bottom) values of r smooth compared to our standard setup.  \n<!-- image -->  \nwith the same scaling of the color bar. While the overall structure of the flow is similar, the rotational feature has moved slightly outward and rotation has decreased. These findings suggests that higher resolution is important for the evolution of the inner core and such simulations should be conducted also for longer times, once this is computationally feasible.  \nFor the simulation RT-M2-VLO we conducted test runs with half of the standard smoothing length r smooth = 0 : 5 r c ( Ssm ) and twice r smooth = 2 r c ( Lsm ). A mid-plane view of the flow and densities is shown in Fig. (A3). As expected, the weaker smoothing results in increased densities close to the core, while stronger smoothing allows for lower densities compared to the standard case. Similar to our standard smoothing, rotation always stays sub-Keplerian. However, the sense of the closed streamlines that are embedding the core changes with different smoothing parameters. For both lower and higher values we find prograde rotation at a similar distance from the core, contrasting the retrograde feature reported above. Clearly, further efforts should be made to implement a true Newtownian potential in future studies of the inner core. We note that the overall atmosphere structure is not substantially affected by the choice of r smooth .  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}