Issue 
A&A
Volume 646, February 2021



Article Number  L11  
Number of page(s)  5  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202040220  
Published online  12 February 2021 
Letter to the Editor
Steady state by recycling prevents premature collapse of protoplanetary atmospheres
^{1}
Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: tobias.moldenhauer@unituebingen.de
^{2}
Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, Haidian District, Beijing, PR China
Received:
23
December
2020
Accepted:
24
January
2021
Context. In recent years, space missions such as Kepler and TESS have discovered many closein planets with significant atmospheres consisting of hydrogen and helium: miniNeptunes. This indicates that these planets formed early in gasrich disks while avoiding the runaway gas accretion that would otherwise have turned them into hotJupiters. A solution is to invoke a long KelvinHelmholtz contraction (or cooling) timescale, but it has also been suggested that thermodynamical cooling can be prevented by hydrodynamical planet atmospheredisk recycling.
Aims. We investigate the efficacy of the recycling hypothesis in preventing the collapse of the atmosphere, check for the existence of a steady state configuration, and determine the final atmospheric mass to core mass ratio.
Methods. We use threedimensional radiationhydrodynamic simulations to model the formation of planetary protoatmospheres. Equations are solved in a local frame centered on the planet.
Results. Ignoring small oscillations that average to zero over time, the simulations converge to a steady state where the velocity field of the gas becomes constant in time. In a steady state, the energy loss by radiative cooling is fully compensated by the recycling of the low entropy gas in the planetary atmosphere with high entropy gas from the circumstellar disk.
Conclusions. For closein planets, recycling naturally halts the cooling of planetary protoatmospheres, preventing them from contracting toward the runaway regime and collapsing into gas giants.
Key words: hydrodynamics / protoplanetary disks / planets and satellites: atmospheres / planets and satellites: formation
© ESO 2021
1. Introduction
Closein planets are a new exoplanet population that were discovered in abundance by the Kepler mission and now also by the Transiting Exoplanet Survey Satellite (TESS). They tend to be small, up to several Earth radii in size, and are located close to their host stars in systems with a fair degree of uniformity (Weiss et al. 2018). Above all, they are very common: one out of three solartype stars may host these types of planets (Winn & Fabrycky 2015; Zhu et al. 2018). Important to this work, their atmospheres are inferred to harbor significant amounts of hydrogen and helium (Wu & Lithwick 2013), strongly suggesting that they formed early within the primordial disk, hence the name “miniNeptune”. The observed valley in the planet size histogram (Fulton et al. 2017) further suggests that even those planets without an atmosphere were born in the disks, with the smallest and most irradiated of them having their atmospheres stripped over time (Jin & Mordasini 2018).
From a theoretical perspective, there is no consensus on how exactly the atmosphere was accreted. Because these planets orbit so close to their host star, they likely assembled in a very short time (Lee et al. 2014). Even if they did not form in situ but migrated to their present orbital radii, it means that their protoatmospheres had copious time to contract on the KelvinHelmholtz timescale, increasing their atmosphere mass to eventually become selfgravitating and form hotJupiters (Batygin et al. 2016). The rapid thermal evolution follows from our expectation that these atmospheres are characterized by a low opacity due to the fact that grains settle quickly (Ormel 2014; Mordasini 2014) and are not replenished by planetesimal infall (Movshovitz et al. 2010). But for miniNeptunes (and superEarths), which by far outnumber hotJupiters, this runaway gas accretion must have been avoided. Proposed solutions include a high opacity atmosphere (Lee et al. 2014), late formation (Lee & Chiang 2015), atmospheric recycling (Ormel et al. 2015; AliDib et al. 2020), and gap opening (Fung & Lee 2018; Ginzburg & Chiang 2019).
In this Letter, we report that planet atmospheredisk interaction yields a steady state and results in the complete recycling of the atmosphere. The recycling hypothesis states that atmospheres are not isolated from their disks, which may be inferred (erroneously) from onedimensional or twodimensional modeling. Isothermal threedimensional simulations by Ormel et al. (2015), Fung et al. (2015), and Béthune & Rafikov (2019) instead have shown that disk gas flows deep into the atmosphere and vice versa. However, the flow structure was seen to strongly depend on the thermodynamical properties of the gas, that is, its cooling rate (Kurokawa & Tanigawa 2018; Popovas et al. 2018). Therefore, this problem is best addressed with radiation hydrodynamical (RHD) simulations (Ayliffe & Bate 2009; Cimerman et al. 2017). Here, we present the results of new threedimensional high resolution RHD simulations of an Earthmass core embedded in a disk and demonstrate that the atmosphere reaches an equilibrium. A true steady state is retrieved, in which the atmosphere mass reaches an asymptotic value, because gas elements are transported back to the disk before they have the chance to cool.
2. Model setup
To study the physics of protoplanetary atmospheres, we followed Cimerman et al. (2017) and ran RHD simulations on a threedimensional grid in spherical coordinates in a local shearing frame centered on the planet. We inserted a planetary core with a mass M_{c} = 1 M_{Earth} at a distance a = 0.1 au from the star in the protoplanetary disk and followed the formation of an atmosphere around the core. Following Lee et al. (2014), assuming a minimummass extrasolar nebula (MMEN, Chiang & Laughlin 2013), we adopted an ambient midplane density of ρ_{0} = 6 × 10^{−6} g cm^{−3} and an ambient temperature of T_{0} = 1000 K. The initial temperature is constant, and the density is vertically Gaussian. For the rocky core, we took a mean density of ρ_{c} = 5 g cm^{−3}. In units of the thermal mass (Rafikov 2006), the planetary core mass corresponds to a dimensionless mass of m = M_{c}/M_{th} = 0.18. For the time unit, we used the inverse Keplerian orbital frequency of the planet with respect to the star, .
In our numerical setup, we followed Cimerman et al. (2017) and used the open source PLUTO code (Mignone et al. 2007) and the fluxlimited diffusion solver from Kuiper et al. (2010, 2020). Similar to Ormel et al. (2015) and Cimerman et al. (2017), we adopted a spherical grid centered on the planet, but contrary to Cimerman et al. (2017), we simulated the full upper hemisphere around the core, assuming symmetry with respect to the equatorial plane. Previously, we had used the additional symmetry of the local shear and restricted the computational domain to half of the upper hemisphere. However, in this case, gravitational smoothing was necessary to avoid numerical instabilities caused by the boundaries together with the radiative transfer. By simulating the full upper hemisphere in our new setup, the uncorrected Newtonian potential for the gravity of the planetary core can be used throughout, which stabilizes the flow structure of the inner atmosphere. Our new simulations therefore no longer feature these unphysical flow artifacts.
For a steady state to emerge, the simulations have to exceed the KelvinHelmholtz timescale. Hence, in order to keep the cooling timescale short and thus save on computation time, we used a low constant opacity of κ = 10^{−4} cm^{2} g^{−1}. At that opacity, the protoplanetary atmosphere is still optically thick. For the spherical grid, we used 128 logarithmically spaced grid cells in the radial direction, 32 equidistant cells in the polar direction, and 128 cells in the azimuthal direction. We tested for convergence by comparing the density profile against a simulation with half resolution. In the lower resolution simulation, the density profile differs only slightly in its evolution and settles into the same steady state. As a comparison, and to demonstrate the impact of recycling, we also ran a onedimensional simulation with otherwise identical parameters and physics.
For the inner boundary at the planetary core, we used a closed boundary condition where all fluxes, including the radiative, vanish (i.e., no energy is exchanged with the planetary core). At the outer boundary, the hydrodynamical quantities were kept constant at their initial values to account for a circumstellar disk that evolves on much longer timescales. The coordinate system is oriented in such a way that the vertical, θ = 0, axis points in the direction of the disk’s angular momentum vector. In a local Cartesian system, this corresponds to the zaxis, where the xaxis points to the radial direction (away from the star) and the yaxis in the direction of the disk’s motion.
In order to prevent initial shocks at the planetary core surface and allow for a smooth simulation start, we multiplied the planet’s gravity by a switchon factor, , where is the injection time. For later usage, we defined the Hill radius as R_{Hill} = a(M_{c}/(3M_{⋆}))^{1/3}.
3. Results
We first present the comparison of the threedimensional simulation to a onedimensional simulation with identical physical parameters to highlight that a threedimensional simulation behaves fundamentally differently. In Sect. 3.2, we then demonstrate the steady state found in the threedimensional simulation by analyzing the entropy evolution and show that the atmosphere eventually stops cooling. In Sect. 3.3, we then show, using tracer particles, that the steady state is indeed maintained by the recycling of the atmospheric gas with gas from the circumstellar disk, in agreement with the recycling hypothesis.
3.1. Onedimensional versus threedimensional comparison
Figure 1 displays the mass inside the Hill sphere as a function of time since the start of the simulation for the threedimensional simulation and a onedimensional simulation with identical parameters. The onedimensional simulation starts with a slightly higher mass because, unlike the threedimensional simulation, it does not include the vertical stratification of the circumstellar disk. Initially, during the switchon time t_{inj}, both simulations behave similarly. However, afterward the onedimensional simulation continues to accrete mass, while the threedimensional simulation stops accreting. In the threedimensional simulation, the final ratio of the mass inside the Hill sphere and the core mass, χ = M_{final}/M_{c}, converges to χ = 2.76%, a value that is typical for a miniNeptune planet. In the onedimensional simulation, gas accretes on a timescale shorter than 1000 years. A mass ratio of χ = 100% is reached at , and χ = 200% is reached at . This behavior suggests that a pure threedimensional mechanism is responsible for counteracting the radiative cooling and prevents the atmosphere from contracting.
Fig. 1.
Mass evolution comparison inside the Hill sphere for the onedimensional (blue) simulation and the threedimensional (orange) simulation. The onedimensional simulation starts with a slightly higher mass than the threedimensional simulation because it does not account for the vertical stratification of the circumstellar disk. While the onedimensional simulation continues to accrete, the threedimensional simulation reaches a steady state, where the final mass, χ = M_{final}/M_{c} = 2.76%, is typical for a miniNeptune planet. 
3.2. Steady state
In Fig. 2, we display the thermal evolution of the threedimensional simulation. The left panel shows the evolution of the radial entropy profile from the start of the simulation up until the maximum simulation time of . We show the entropy profile because adiabatic compression and expansion do not change the entropy of the gas, and therefore the entropy is a direct indicator for the cooling of the gas. The entropy profile lines are spaced logarithmically in time. Initially, the entropy drops rapidly because of the radiative cooling. As explained in Sect. 2, the adopted opacity is very low in order to facilitate a cooling timescale shorter than the simulation time. After approximately , the entropy profile stops decreasing and becomes constant in time, ignoring small oscillations. These oscillations are displayed in the right panel of Fig. 2. There we see three distinctive radii, r ∈ [0.1, 0.3, 1] R_{Hill}, and compare the temperature to the mean temperature averaged over the steady state (). With a maximum amplitude of 0.15% in temperature and a maximum amplitude of 0.3% in density, these oscillations are very weak and therefore not important for the evolution of the atmosphere. Toward the inner boundary, the oscillations are strongly periodic, in contrast to the outer regions where the oscillations are more irregular. The low frequency of the oscillations, ∼0.038 Ω_{K}, is not consistent with the high local Brunt–Väisäla frequency (∼100 Ω_{K}), while the period is comparable to the sound crossing time through the computational domain. Hence, we conclude that they are caused by small global residual perturbations in our inviscid simulation. The statistical trend of these oscillations is consistent with zero. Hence, we can conclude that the threedimensional simulation is in a true steady state in which the radiative cooling is fully compensated.
Fig. 2.
Thermal evolution of the threedimensional simulation. Left: time evolution of the spherically averaged entropy. The lines are spaced logarithmically in time and the last four lines overlap, indicating the emergence of a radial entropy profile that is constant in time. Right: relative temperature deviation to the mean temperature in the steady state at three distances to the protoplanetary core. The trend line slopes of (− 0.84 ± 4.21) × 10^{−6}% Ω_{K} at 0.1 R_{Hill}, (− 4.21 ± 1.12) × 10^{−6}% Ω_{K} at 0.3 R_{Hill}, and (− 2.28 ± 0.31) × 10^{−5}% Ω_{K} at 1 R_{Hill} are consistent with a true steady state. 
3.3. Recycling hypothesis
Energy may only be transported by radiation and advection. Radiation transports energy outward from the hotter inner regions of the atmosphere toward the cooler circumstellar disk. Therefore, advection must counter the radiation energy flux and transport energy back inward to the planetary core. In a steady state, where the net mass flux is zero, adiabatic compression and expansion cancel each other out. Even though the net mass flux over a shell vanishes, this does not imply that the advection energy flux vanishes as well. The energy of the gas that enters the atmosphere is higher than the energy of the gas that leaves the atmosphere, resulting in a net advection energy flux inward. In one dimension, there is only one possible direction for the gas to flow, which makes recycling physically impossible.
To further investigate the recycling of the atmospheric gas, we used tracer particles to measure the time the gas has spent inside the Hill sphere, which we will refer to as the “recycling time”. The recycling time was calculated by inserting tracer particles throughout the Hill sphere. In order to remove the small oscillations shown in Fig. 2, we timeaveraged the velocity field in the steady state before introducing the tracer particles. We then integrated them backward in time along streamlines and measured the time it takes for them to leave the Hill sphere. This recycling time measures how much time the tracer particle has spent inside the Hill sphere; it quantifies how rapidly atmospheric gas and disk gas are exchanged. Because the trajectory, and thus the measured recycling time, is highly sensitive to the starting location of the particles, we used many more particles than grid cells and took the median of adjacent particles. This approach is justified because the particles represent a continuous fluid. Figure 3 displays twodimensional cuts of this recycling time. As expected, the recycling time is significantly higher in the inner part of the Hill sphere, where the gas rotates around the protoplanetary core, than in the outer regions where it flows past the planet. Most importantly, however, almost the whole atmosphere recycles in a finite time frame. Even the region directly at the surface of the planetary core is recycled. Only a small region at the boundary of the vertical influx and the outflow in the midplane, the yellow ring in Fig. 3, did not recycle even after . However, this small remaining part is negligible for the thermodynamic structure of the atmosphere. Therefore, the gas that forms the planetary atmosphere is eventually replenished with fresh, high entropy gas from the circumstellar disk.
Fig. 3.
Recycling time and density of the gas inside the Hill sphere. The length unit is R_{Hill}. The velocity of the gas is measured in units of the Keplerian orbital velocity of the planet around the star. The recycling time is defined as the time the gas has already spent inside the Hill sphere. The arrows indicate the velocity of the gas. The density contour lines (cyan) are drawn at [2, 4, 8, 16] ρ_{0}. These lines are jagged because of the resolution of the grid. Left: xy plane or midplane. Middle: xz plane. Right: yz plane. The atmosphere is recycled up to the surface of the protoplanetary core. Gas around the zaxis has a lower recycling time than the gas in the midplane. This indicates that gas enters from the poles and leaves the atmosphere through the midplane. 
In addition to the recycling time, the figure also displays the density distribution using isodensity lines. The density at the inner boundary reaches a maximum value of ρ_{max} = 117 ρ_{0}. We can see that both the long recycling time region and the high density region are squished vertically and along the orbital direction. However, the long recycling time region is also elongated along the xdirection because of the stellar tidal forces. Close to the vertical axis, the recycling time is also significantly shorter than in the midplane. This shows that fresh gas enters the atmosphere through the poles and leaves the atmosphere by spiraling outward in the midplane before it gets picked up by the shearing flow.
Figure 4 shows the massweighted recycling time (i.e., the recycled mass) and compares it to the total mass inside the Hill sphere. Eventually, over 99.85% of the mass inside the Hill sphere is recycled with fresh gas from the circumstellar disk. The remaining part is heated radiatively by the surrounding gas and therefore does not continue to cool either. This illustrates that the radiative cooling is fully compensated by the recycling of the atmospheric gas with fresh gas from the circumstellar disk.
Fig. 4.
Recycled mass in the Hill sphere as a function of time after a steady state has emerged. Most mass in the Hill sphere is part of the shearing flow and is therefore recycled very quickly (). For the inner parts of the atmosphere to be recycled, it takes approximately . A small part of the atmosphere is not recycled even after our maximum integration time of . 
4. Discussion and summary
It has been customary to model planets’ thermal evolution like stars, using the onedimensional stellar structure equations. But young planets, embedded in their natal disks, cannot be regarded as isolated from their environment. Facing the disk gas flow, planets need to “process” material at a rate of – equivalent to ∼1 Earth mass per year (!) – for the parameters used in this study. This number, although crude and not accounting for the spatial variations in the recycling efficiency, in itself underscores the fact that embedded planets cannot be treated as onedimensional isolated objects. In this Letter, with improved threedimensional RHD simulations, we have shown that the entire envelope is affected by the vigorous recycling, preventing the adiabatically heated gas from cooling radiatively. We found this despite our choice of very low opacity, for which the onedimensional setup would enter runaway gas accretion in a mere 10^{3} yr. In a future work, we will investigate the dependence of our result on parameters such as the planet mass, location, and disk headwind. In particular, for typical disk parameters with increasing orbital radii, the impact rate of disk gas () decreases, while the gravitational potential at the planet surfaces steepens (∼R_{Hill}/R). Therefore, we can expect recycling to become less vigorous in the outer disk. Nonetheless, the ubiquity of superEarths and miniNeptunes close to their host star is most naturally explained by the recycling mechanism. Gas hydrodynamics trumps radiative cooling.
Acknowledgments
TWM and RK acknowledge funding through the German Research Foundation (DFG) under grant no. KU 2849/61 as well as support through travel grants under the SPP 1992: Exoplanet Diversity program. RK acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the DFG under grant no. KU 2849/31 and KU 2849/32. We acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of BadenWürttemberg through bwHPC and the DFG through grant no. INST 37/9351 FUGG.
References
 AliDib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440 [Google Scholar]
 Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49 [Google Scholar]
 Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114 [Google Scholar]
 Béthune, W., & Rafikov, R. R. 2019, MNRAS, 488, 2365 [Google Scholar]
 Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 [Google Scholar]
 Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
 Fung, J., & Lee, E. J. 2018, ApJ, 859, 126 [Google Scholar]
 Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101 [Google Scholar]
 Ginzburg, S., & Chiang, E. 2019, MNRAS, 487, 681 [Google Scholar]
 Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., Yorke, H. W., & Mignone, A. 2020, ApJS, 250, 13 [Google Scholar]
 Kurokawa, H., & Tanigawa, T. 2018, MNRAS, 479, 635 [Google Scholar]
 Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41 [Google Scholar]
 Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mordasini, C. 2014, A&A, 572, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [Google Scholar]
 Ormel, C. W. 2014, ApJ, 789, L18 [Google Scholar]
 Ormel, C. W., Shi, J.M., & Kuiper, R. 2015, MNRAS, 447, 3512 [Google Scholar]
 Popovas, A., Nordlund, Å., Ramsey, J. P., & Ormel, C. W. 2018, MNRAS, 479, 5136 [Google Scholar]
 Rafikov, R. R. 2006, ApJ, 648, 666 [Google Scholar]
 Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [Google Scholar]
 Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [Google Scholar]
 Wu, Y., & Lithwick, Y. 2013, ApJ, 772, 74 [Google Scholar]
 Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101 [Google Scholar]
All Figures
Fig. 1.
Mass evolution comparison inside the Hill sphere for the onedimensional (blue) simulation and the threedimensional (orange) simulation. The onedimensional simulation starts with a slightly higher mass than the threedimensional simulation because it does not account for the vertical stratification of the circumstellar disk. While the onedimensional simulation continues to accrete, the threedimensional simulation reaches a steady state, where the final mass, χ = M_{final}/M_{c} = 2.76%, is typical for a miniNeptune planet. 

In the text 
Fig. 2.
Thermal evolution of the threedimensional simulation. Left: time evolution of the spherically averaged entropy. The lines are spaced logarithmically in time and the last four lines overlap, indicating the emergence of a radial entropy profile that is constant in time. Right: relative temperature deviation to the mean temperature in the steady state at three distances to the protoplanetary core. The trend line slopes of (− 0.84 ± 4.21) × 10^{−6}% Ω_{K} at 0.1 R_{Hill}, (− 4.21 ± 1.12) × 10^{−6}% Ω_{K} at 0.3 R_{Hill}, and (− 2.28 ± 0.31) × 10^{−5}% Ω_{K} at 1 R_{Hill} are consistent with a true steady state. 

In the text 
Fig. 3.
Recycling time and density of the gas inside the Hill sphere. The length unit is R_{Hill}. The velocity of the gas is measured in units of the Keplerian orbital velocity of the planet around the star. The recycling time is defined as the time the gas has already spent inside the Hill sphere. The arrows indicate the velocity of the gas. The density contour lines (cyan) are drawn at [2, 4, 8, 16] ρ_{0}. These lines are jagged because of the resolution of the grid. Left: xy plane or midplane. Middle: xz plane. Right: yz plane. The atmosphere is recycled up to the surface of the protoplanetary core. Gas around the zaxis has a lower recycling time than the gas in the midplane. This indicates that gas enters from the poles and leaves the atmosphere through the midplane. 

In the text 
Fig. 4.
Recycled mass in the Hill sphere as a function of time after a steady state has emerged. Most mass in the Hill sphere is part of the shearing flow and is therefore recycled very quickly (). For the inner parts of the atmosphere to be recycled, it takes approximately . A small part of the atmosphere is not recycled even after our maximum integration time of . 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.