Free Access
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/0004-6361/202040220
Published online 12 February 2021

© ESO 2021

1. Introduction

Close-in 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 solar-type 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 “mini-Neptune”. 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 proto-atmospheres had copious time to contract on the Kelvin-Helmholtz timescale, increasing their atmosphere mass to eventually become self-gravitating and form hot-Jupiters (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 mini-Neptunes (and super-Earths), which by far outnumber hot-Jupiters, 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; Ali-Dib et al. 2020), and gap opening (Fung & Lee 2018; Ginzburg & Chiang 2019).

In this Letter, we report that planet atmosphere-disk 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 one-dimensional or two-dimensional modeling. Isothermal three-dimensional 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 three-dimensional high resolution RHD simulations of an Earth-mass 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 three-dimensional grid in spherical coordinates in a local shearing frame centered on the planet. We inserted a planetary core with a mass Mc = 1 MEarth 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 minimum-mass 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 T0 = 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 M th = c s 3 / ( G Ω K ) $ M_{\rm th}=c_{\rm s}^{3}/(G\Omega_{\rm K}) $ (Rafikov 2006), the planetary core mass corresponds to a dimensionless mass of m = Mc/Mth = 0.18. For the time unit, we used the inverse Keplerian orbital frequency of the planet with respect to the star, Ω K 1 = a 3 / ( G M ) $ \Omega_{\mathrm{K}}^{-1} = \sqrt{a^3 / (G M_\star)} $.

In our numerical setup, we followed Cimerman et al. (2017) and used the open source PLUTO code (Mignone et al. 2007) and the flux-limited 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 Kelvin-Helmholtz 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 cm2 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 one-dimensional 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 z-axis, where the x-axis points to the radial direction (away from the star) and the y-axis 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 switch-on factor, α inj ( t ) = 1 exp ( 0.5 t 2 / t inj 2 ) $ \alpha_{\rm inj}(t) = 1-\exp(-0.5 t^2/t_{\rm inj}^2) $, where t inj = 4 Ω K 1 $ t_{\mathrm{inj}} = 4\,\Omega_{\mathrm{K}}^{-1} $ is the injection time. For later usage, we defined the Hill radius as RHill = a(Mc/(3M))1/3.

3. Results

We first present the comparison of the three-dimensional simulation to a one-dimensional simulation with identical physical parameters to highlight that a three-dimensional simulation behaves fundamentally differently. In Sect. 3.2, we then demonstrate the steady state found in the three-dimensional 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. One-dimensional versus three-dimensional comparison

Figure 1 displays the mass inside the Hill sphere as a function of time since the start of the simulation for the three-dimensional simulation and a one-dimensional simulation with identical parameters. The one-dimensional simulation starts with a slightly higher mass because, unlike the three-dimensional simulation, it does not include the vertical stratification of the circumstellar disk. Initially, during the switch-on time tinj, both simulations behave similarly. However, afterward the one-dimensional simulation continues to accrete mass, while the three-dimensional simulation stops accreting. In the three-dimensional simulation, the final ratio of the mass inside the Hill sphere and the core mass, χ = Mfinal/Mc, converges to χ = 2.76%, a value that is typical for a mini-Neptune planet. In the one-dimensional simulation, gas accretes on a timescale shorter than 1000 years. A mass ratio of χ = 100% is reached at t = 35 000 Ω K 1 $ t=35\,000\,\Omega_{\mathrm{K}}^{-1} $, and χ = 200% is reached at t = 65 000 Ω K 1 $ t=65\,000\,\Omega_{\mathrm{K}}^{-1} $. This behavior suggests that a pure three-dimensional mechanism is responsible for counteracting the radiative cooling and prevents the atmosphere from contracting.

thumbnail Fig. 1.

Mass evolution comparison inside the Hill sphere for the one-dimensional (blue) simulation and the three-dimensional (orange) simulation. The one-dimensional simulation starts with a slightly higher mass than the three-dimensional simulation because it does not account for the vertical stratification of the circumstellar disk. While the one-dimensional simulation continues to accrete, the three-dimensional simulation reaches a steady state, where the final mass, χ = Mfinal/Mc = 2.76%, is typical for a mini-Neptune planet.

3.2. Steady state

In Fig. 2, we display the thermal evolution of the three-dimensional 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 2000 Ω K 1 $ 2000\,\Omega_{\mathrm{K}}^{-1} $. 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 200 Ω K 1 $ 200\,\Omega_{\mathrm{K}}^{-1} $, 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] RHill, and compare the temperature to the mean temperature averaged over the steady state ( t > 500 Ω K 1 $ t > 500\,\Omega_{\mathrm{K}}^{-1} $). 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 three-dimensional simulation is in a true steady state in which the radiative cooling is fully compensated.

thumbnail Fig. 2.

Thermal evolution of the three-dimensional 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 RHill, (− 4.21 ± 1.12) × 10−6% ΩK at 0.3 RHill, and (− 2.28 ± 0.31) × 10−5% ΩK at 1 RHill 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 time-averaged 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 two-dimensional 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 10 4 Ω K 1 $ 10^{4}\,\Omega_{\mathrm{K}}^{-1} $. 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.

thumbnail Fig. 3.

Recycling time and density of the gas inside the Hill sphere. The length unit is RHill. 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 z-axis 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 iso-density 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 x-direction 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 mass-weighted 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.

thumbnail 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 ( 10 Ω K 1 $ {\approx}10\,\Omega_{\mathrm{K}}^{-1} $). For the inner parts of the atmosphere to be recycled, it takes approximately 10 3 Ω K 1 $ 10^{3}\,\Omega_{\mathrm{K}}^{-1} $. A small part of the atmosphere is not recycled even after our maximum integration time of 10 4 Ω K 1 $ 10^{4}\,\Omega_{\mathrm{K}}^{-1} $.

4. Discussion and summary

It has been customary to model planets’ thermal evolution like stars, using the one-dimensional 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 R Hill 3 Ω K ρ 0 $ R_{\rm Hill}^3\,\Omega_{\rm K} \rho_0 $ – 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 one-dimensional isolated objects. In this Letter, with improved three-dimensional 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 one-dimensional setup would enter runaway gas accretion in a mere 103 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 ( R Hill 3 Ω K ρ 0 $ {\sim}R_{\rm Hill}^3\,\Omega_{\rm K} \rho_0 $) decreases, while the gravitational potential at the planet surfaces steepens (∼RHill/R). Therefore, we can expect recycling to become less vigorous in the outer disk. Nonetheless, the ubiquity of super-Earths and mini-Neptunes 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/6-1 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/3-1 and KU 2849/3-2. 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 Baden-Württemberg through bwHPC and the DFG through grant no. INST 37/935-1 FUGG.

References

  1. Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440 [Google Scholar]
  2. Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49 [Google Scholar]
  3. Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114 [Google Scholar]
  4. Béthune, W., & Rafikov, R. R. 2019, MNRAS, 488, 2365 [Google Scholar]
  5. Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 [Google Scholar]
  6. Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662 [Google Scholar]
  7. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
  8. Fung, J., & Lee, E. J. 2018, ApJ, 859, 126 [Google Scholar]
  9. Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101 [Google Scholar]
  10. Ginzburg, S., & Chiang, E. 2019, MNRAS, 487, 681 [Google Scholar]
  11. Jin, S., & Mordasini, C. 2018, ApJ, 853, 163 [Google Scholar]
  12. Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Kuiper, R., Yorke, H. W., & Mignone, A. 2020, ApJS, 250, 13 [Google Scholar]
  14. Kurokawa, H., & Tanigawa, T. 2018, MNRAS, 479, 635 [Google Scholar]
  15. Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41 [Google Scholar]
  16. Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95 [Google Scholar]
  17. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  18. Mordasini, C. 2014, A&A, 572, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [Google Scholar]
  20. Ormel, C. W. 2014, ApJ, 789, L18 [Google Scholar]
  21. Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 [Google Scholar]
  22. Popovas, A., Nordlund, Å., Ramsey, J. P., & Ormel, C. W. 2018, MNRAS, 479, 5136 [Google Scholar]
  23. Rafikov, R. R. 2006, ApJ, 648, 666 [Google Scholar]
  24. Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [Google Scholar]
  25. Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [Google Scholar]
  26. Wu, Y., & Lithwick, Y. 2013, ApJ, 772, 74 [Google Scholar]
  27. Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101 [Google Scholar]

All Figures

thumbnail Fig. 1.

Mass evolution comparison inside the Hill sphere for the one-dimensional (blue) simulation and the three-dimensional (orange) simulation. The one-dimensional simulation starts with a slightly higher mass than the three-dimensional simulation because it does not account for the vertical stratification of the circumstellar disk. While the one-dimensional simulation continues to accrete, the three-dimensional simulation reaches a steady state, where the final mass, χ = Mfinal/Mc = 2.76%, is typical for a mini-Neptune planet.

In the text
thumbnail Fig. 2.

Thermal evolution of the three-dimensional 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 RHill, (− 4.21 ± 1.12) × 10−6% ΩK at 0.3 RHill, and (− 2.28 ± 0.31) × 10−5% ΩK at 1 RHill are consistent with a true steady state.

In the text
thumbnail Fig. 3.

Recycling time and density of the gas inside the Hill sphere. The length unit is RHill. 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 z-axis 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
thumbnail 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 ( 10 Ω K 1 $ {\approx}10\,\Omega_{\mathrm{K}}^{-1} $). For the inner parts of the atmosphere to be recycled, it takes approximately 10 3 Ω K 1 $ 10^{3}\,\Omega_{\mathrm{K}}^{-1} $. A small part of the atmosphere is not recycled even after our maximum integration time of 10 4 Ω K 1 $ 10^{4}\,\Omega_{\mathrm{K}}^{-1} $.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.