Formation of Super-Earths by Tidally-Forced Turbulence

The Kepler observations indicate that many exoplanets are super-Earths, which brings about a puzzle for the core-accretion scenario. Since observed super-Earths are in the range of critical mass, they would accrete gas efficiently and become gas giants. Theoretically, super-Earths are predicted to be rare in the core-accretion framework. To resolve this contradiction, we propose that the tidally-forced turbulent diffusion may affect the heat transport inside the planet. Thermal feedback induced by turbulent diffusion is investigated. We find that the tidally-forced turbulence would generate pseudo-adiabatic regions within radiative zones, which pushes the radiative-convective boundaries (RCBs) inwards. This would decrease the cooling luminosity and enhance the Kelvin-Helmholtz (KH) timescale. For a given lifetime of protoplanetary disks (PPDs), there exists a critical threshold for the turbulent diffusivity, $\nu_{\rm critical}$. If $\nu_{\rm turb}>\nu_{\rm critical} $, the KH timescale is longer than the disk lifetime and the planet would become a super-Earth rather than a gas giant. We find that even a small value of turbulent diffusion has influential effects on evolutions of super-Earths. $\nu_{\rm critical}$ increases with the core mass. We further ascertain that, within the minimum mass extrasolar nebula (MMEN), $\nu_{\rm critical}$ increases with the semi-major axis. This may explain the feature that super-Earths are common in inner PPD regions, while gas giants are common in the outer PPD regions. The predicted envelope mass fraction (EMF) is not fully consistent with observations. We discuss physical processes, such as late core assembly and mass loss mechanisms, that may be operating during super-Earth formation.


Introduction
One of the widely accepted mechanisms of planet formation is the core-nucleated instability theory (Perri & Cameron 1974;Haris 1978;Mizuno et al. 1978;Stevenson 1982). According to this scenario, the massive gaseous atmosphere would be accumulated in a runaway manner when the core mass reaches a critical value. In static models, when heating balances cooling, runaway accretion occurs when the planet is beyond the critical mass, because the envelope fails to hold hydrostatic equilibrium. Rafikov (2006) found a broad range of critical mass (0.1M ⊕ ≤ M critical ≤ 100M ⊕ ) due to various disk properties and planetesimal accretion rate. However, in dynamic or quasi-static models, the thermal disequilibrium rather than the hydrostatic disequilibrium plays the dominant role. The runaway accretion occurs because the envelope becomes thermally unstable as the cooling timescale becomes catastrophically shorter. In this case, the runaway accretion is driven by runaway cooling (Bodenheimer & Pollack 1986;Lee et al. 2014;Piso & Youdin 2014). Three stages are involved in the formation process. In the first stage, rocky cores grow by rapid planetesimal accretion. In the second stage, the core's feeding zone is depleted of solids and the atmosphere grows gradually, regulated by the KH contraction. Finally, when the atmosphere reaches the crossover mass, gas runaway takes place and the planet gets inflated into a gas giant. The timescale of second stage is the longest among the three and dominates whole formation process (Pollack et al. 1996).
About 20% of Sun-like stars host super-Earths with radii of 1-4 R ⊕ at distance 0.05-0.3 AU (Howard et al. 2010;Batalha et al. 2013;Petigura et al. 2013). Radial velocity measurements (Weiss & Marcy 2014) and transit timing variations (Wu & Lithwick 2013) manifest that masses of these super-Earths are in the range of 2−20 M ⊕ . The abundance of super-Earths presents a puzzle for the core instability theory. This theory indicates that when a protoplanet reaches the super-Earth size, two physical processes make the survival of super-Earths difficult, leading to a planetary "desert" in this size range (Ida & Lin 2004). Super-Earths would excite density waves in PPDs and give rise to rapid type I migration. This type of migration would cause the planet to be engulfed by its host star if the disk inner edge touches the stellar surface. Recent studies have sought various remedies for type I migration (Yu et al. 2010;Fung & Chiang 2017). PPDs are expected to have an inner edge at the stellar magnetosphere (e.g., Long et al. 2005). For planets undergoing disk-driven migration, they are expected to pile up near this edge. They would stay either at the edge because the gas runs out, or inside the edge down to 2:1 resonance because that's where the tidal torque will taper off, or outside the edge as the standing waves generated by wave reflection off the inner edge stall planet migration (Tsang 2011). In this paper we would focus on another threat for super-Earths. Super Earths have low mean density, which suggests that they must be surrounded by gas envelopes (Rogers & Seager 2010). Since these observed super-Earths are in the range of critical mass, they would trigger efficient gas runaway and accumulate massive gas envelope. They would become gas giants. As a result, super-Earths are supposed to be rare. However, the Kepler's discovery wreck these predictions. Lee et al. (2014) has proposed metallicity gradient inside the PPDs or late assembly of cores to resolve the puzzle of super-Earth formation. Lee & Chiang (2016) stressed that the late core assembly in transitional PPDs is more consistent with observations. In gas-poor environments, gas dynamical friction has weakened to allow proto-cores to stir one another and merge. In addition, this formation scenario ensures that super-Earth cores accrete mass with a few percent envelope mass fraction (EMF). Guillot & Showman (2002) argued that the dissipation of kinetic energy of atmospheric wind, driven by intense irradiation, could bury heat inside the planet. Many studies extend this idea to explain the radius anomaly of hot Jupiters (Youdin & Mitchell 2010;Ginzburg & Sari 2015;Komacek & Youdin 2017). These investigations focus on the late evolution after the disk dispersal. Unfortunately, this is invalid for the early evolution of super-Earths because they are still embedded within disks. The irradiation may not penetrate the disk and is not able to bury heat in the exoplanets.
However, we note that tidal interactions between the host star and planet can periodically perturb the planet and generate mechanical forcing of the fluid motions (Zahn 1977;Goldreich & Nicholson 1989). Heating by tidal dissipation in primordial super-Earth envelope can inhibit the gas cooling (Ginzburg & Sari 2017). This mechanism requires the orbital eccentricity of super-Earths be continuously pumped. But super-Earths may not be massive enough to clear a clean gap to excite orbital eccentricity (Goldreich & Sari 2003). Another important aspect about tidal interaction is that tidally-forced turbulent mixing would induce heat transport inside the planets. Recent laboratory experiment shows that turbulence could penetrate deep inside the planet interior (Cabanes et al. 2017). By combining laboratory measurements and high resolution simulations, Grannan et al. (2017) confirmed the generation of bulk filling turbulence inside planet driven by tidal forcing. Turbulent mixing plays an essential role in heat transport in strongly stratified environments (Garaud & Kulenthirarajah 2016). This motivates us to study the effects of turbulent diffusion on the planet's thermal evolution. Prior studies have noticed that the turbulent mixing induced by mechanical forcing leads to heat transport inside hot Jupiters (Youdin & Mitchell 2010). These tides would produce appreciable thermal feedback and may lead to interior radiative zones, enhancing g-mode dissipations with a wide spectrum of resonances (Jermyn et al. 2017). We find that the thermal feedback associated with the externally-forced turbulent stirring may greatly alter the accretion history of super-Earths.
It is well known that the timescale of gas accretion is dictated by the KH timescale. In other words, accretion is determined by the planet's ability to cool (Lee & Chiang 2015). In this paper, we note that the tidally-forced turbulent diffusion influences the heat transport inside the planet's envelope. Thermal feedback would be induced by turbulent diffusion. The heat transport associated with tidally-forced turbulent diffusion would reduce the cooling luminosity and enhance the KH timescale. We find that turbulent diffusion may have significant effects on the planet accretion history 1 . Based on our calculations, we propose that tidally-forced turbulent diffusion would effectively help super-Earths evade growing into gas giants. This paper is structured as follows. In section 2, we provide a brief description of the accreting planet envelope with tidally-forced turbulent diffusion. In section 3, we compare the planet interior thermal profile with and without turbulent diffusion, discussing the thermal feedback induced by turbulent diffusion, especially the shift of RCBs. In section 4, we depict the cooling luminosity variations and onset of gas runaway. The quasi-static Kelvin-Helmholtz evolution and critical turbulent diffusivity are discussed in Section 5. In Section 6, we discuss the mass loss mechanisms for super-Earths and the limitation of super-Earth formation by the turbulent diffusion. Summary and conclusions are given in Section 7.

Accreting Envelope with Tidally-Forced Turbulence
Super-Earths are susceptible to runaway accretion (Pollack et al. 1996). The ability to accrete is determined by the planet's power to cool (Lee & Chiang 2015). How super-Earths avoid rapid gas runaway depends critically on the cooling history of planet, which is closely related to the thermal structure of the envelope. In the convectively stable region, the turbulent diffusion would induce heat transport within the planet. In this section we will concentrate on the thermal feedback caused by tidally-induced turbulent diffusion.

Thermal structure of Gaseous Envelope
Since the planet's ability to cool depends on planets' thermal structure of the envelope, we first study the gaseous envelope structure of planets, i.e., the distribution of pressure, temperature, and mass around a protoplanetary core with mass M c embedded within the protoplanetary nebular. The planet envelope (or, interchangeably "atmosphere") structure is governed by the following equations of mass conservation, hydrostatic equilibrium, thermal gradient, and energy conservation (Kippenhahn et al. 2012) : where G is the gravitational constant, P is the pressure, ρ is the density, T is the temperature, L is the luminosity, and M r is the mass, including the core mass and the atmosphere mass, enclosed inside the radius r, M r = M atm + M c . The symbol "∇" denotes the temperature gradient inside the envelope. The energy generation ǫ is set to zero since there is no nuclear reaction inside the planet. The above equations implicitly indicate that the envelope quickly adjusts and dynamical timescale is shorter than the accretion timescale (Rafikov 2006). Note that the right hand side term, −T ∂s ∂t , in the energy equation dictates the cooling process. Replacing the local energy equation by a global energy equation would greatly reduce the numerical tasks and we need only deal with ODEs rather than PDEs (Piso & Youdin 2014;Lee et al. 2014). Details will be discussed in Section 3.
The energy transport in the convective region is very efficient and the temperature gradient is 2 The convective and radiative layers of the envelope are specified by the Schwarzschild criterion: the atmosphere is stable against convection when ∇ < ∇ ad and convectively unstable when ∇ ≥ ∇ ad . Since the convective energy transport is efficient, ∇ = ∇ ad in the convective region. The actual temperature gradient can be expressed as In this paper, we adopt a polytropic index γ = 7/5 for an ideal diatomic gas and the adiabatic gradient ∇ ad = (γ − 1)/γ. Note that the realistic equation of state (EOS) would change the value of ∇ ad and the effects of realistic EOS will be left for future studies.
The radiative temperature gradient where κ is the opacity. In the upper part of the atmosphere, the exact value of κ is highly uncertain because the amount of dust and the dust size distribution are not well constrained in PPDs. Lee et al. (2014) studied both dusty and dust-free atmosphere and found that the radiative-convective boundaries (RCBs) are determined by H 2 dissociation at an almost fixed temperature ∼2500 K for dusty atmosphere. They also found the for dust-free atmosphere, the radiative region keeps an almost isothermal temperature fixed by the envelope outer surface. Technically, the opacity laws can be written as a power law as a function of pressure and temperature whether or not the total opacity is dominated by dust grains. For these reasons, we adopt a power law opacity (Rafikov 2006;Piso & Youdin 2014;Ginzburg et al. 2016), by assuming that Here we choose κ 0 = 0.001cm 2 g −1 , which allows our fiducial model without turbulent diffusion to possess properties of more sophisticated super-Earth models (Lee et al. 2014). What is important is the opacity near the RCB. In that sense, it is important to keep in mind that the power-law indices α and β can change significantly within the envelope (and with distance from the star). We have tried different choices of α and β. We find that, as long as the parameter α and β satisfy ∇ 0 ≡ 1+α 4−β > ∇ ad , our results are robust and insensitive to the choices we made 3 .
Conventionally, it is believed that solid cores accrete planetesimal and gas simultaneously (Pollack et al. 1996;Bodenheimer et al. 2000). However, estimation shows that the termination epoch of accretion of solids is well before the accretion of gas. The dust coagulation timescale can be as short as t coagulate ∼ 10 4 yr especially when the planet is close to the central host (Lee et al. 2014). This timescale is much shorter than typical disk dispersal timescale (∼ 0.5−10 Myr). In addition, calculations by Lee & Chiang (2015) showed that planetesimal accretion does not generically prevent runaway. As a result, it is physically valid to set the planetesimal accretion rate to zero (L acc = 0) when we study accreting super-Earths within the disk. In this case, the core is free to cool and contract, and it is extremely susceptible to the gas runaway.
Note that the above differential equations are essentially identical to the usual planet interior structure equations. The distinction is the thermal feedback generated by tidallyforced turbulent mixing inside the stably stratified region. More specifically, ∇ rad is affected by the turbulent diffusion, which will be further discussed in the next section.

Thermal Feedback by Tidally-Forced Turbulent Mixing
How do super-Earths evade becoming gas giants? In this paper, we propose a robust mechanism to avoid runaway accretion. Due to the tidal forcing, the planet's gas envelope would be stirred and the turbulent motion may be initiated. Detailed analyses of these processes are rather complex and beyond the scope of this paper (e.g., Garaud & Kulenthirarajah 2016;Grannan et al. 2017). In this paper, we try to constrain the turbulent diffusion that is necessary to influentially affect the planet accretion timescale. We find that even weak turbulence would affect the planet accretion history significantly.
Since the sound-crossing time is much shorter than the time for heat to diffuse across the fluid blob, the blob conserves entropy (i.e. adiabatically) and keeps pressure equilibrium with the ambient environments when it displaces over a radial distance ℓ. The temperature difference between the blob and its surroundings is The heat excess associated with these fluid blobs can be written as δq = ρc p δT and the corresponding turbulent heat flux is F turb = vδq, where v is the characteristic speed of turbulent eddies. The entropy gradient can be put down as where g is the gravitational acceleration. This equation indicates that in the stably stratified region (∇ < ∇ ad ), the entropy gradient is positive (ds/dr > 0). The heat flux by turbulent mixing is then The flux is negative for stable stratification. For a thermal engine without external forcing, heat always flows from hot to cold regions. However, with external mechanical forcing by tides, heat flows from cold to hot regions becomes feasible (Youdin & Mitchell 2010). Note that the turbulent diffusion coefficient µ turb ≡ vℓ 4 and the corresponding luminosity is The total luminosity is carried by two components, the radiative and the turbulent We note that the temperature gradient in the radiative region can be arranged in a compact form as (see Appendix A for details), In the above equation, and where the superscript "(0)" indicates the radiative temperature gradient without turbulence 5 and M c is the mass of the solid core. It can be readily shown that the following inequality holds in radiative region ∇ (0) rad < ∇ < ∇ ad (see Figure 3 for the pseudo-adiabatic region). Here we stress that it is the turbulent diffusion driven by external tidal forcing that makes ∇ steeper than ∇ (0) rad . This inequality has significant implications for the thermal feedback induced by tidally-forced turbulent diffusion. An interesting issue is that radiative zones would be enlarged and the cooling luminosity would be greatly reduced.
Here we define two dimensionless parameters The two parameters represent the strength of turbulence. In the definition of ζ, H p ≡ −dr/d ln P and c s are pressure scale height and sound speed, respectively. It is obvious that, if the turbulence in the radiative region is negligible, i.e., η = 0, the temperature gradient recovers its usual definition, ∇ rad → ∇ rad . In section 5.1, we will give a physical estimation of the parameter ζ based on our calculations. We will see that small value of ζ ∼ 10 −6 −10 −5 has already appreciable effects on the formation of super-Earths. This mechanism is robust in the sense that even weak turbulence is adequate for it to operate. We should keep in mind that one limitation is that the turbulence strength is parameterized, not physically specified. This is an important issue which still remains to be addressed, i.e., forcing turbulence induced by tides should be investigated in further detail (Barker 2016; Grannan et al. 2017).

Boundary Conditions
The density and temperature at the outer boundary of the atmosphere are given by the nebular density and temperature. We adopt the minimum mass extrasolar nebula (MMEN) model of Chiang & Laughlin (2013). According to MMEN, the disk structure reads, The inner boundary lies at the surface of the inner core. The core density is assumed to be ρ core = 7g cm −3 , the core mass is 5 M ⊕ and the core radius is R core = 1.6 R ⊕ . The outer boundary condition is chosen at the smaller of the Bondi radius and Hill radius, which are respectively.

Thermal Properties of Gas Envelopes
Since the thermal cooling timescale is intimately related to the planet interior structure, we first describe the interior structure of the gaseous envelope. To avoid the complication induced by sandwiched convection-radiation structure inside the planet interior (Ginzburg & Sari 2015;Jermyn et al. 2017), we simply consider a two-layer model, i.e., a convective interior and a radiative exterior (Piso & Youdin 2014).
We adopt the assumption that the luminosity, L, is spatially constant, which is valid in radiative region if the thermal relaxation timescale is shorter than thermal times in the rest of the atmosphere. The validation of such assumption is corroborated by Piso & Youdin (2014) and Lee et al. (2014). To get thermal profiles within the envelope, a luminosity L is required to obtain ∇ rad before we numerically integrate the structure equations. The spatially constant L is treated as an eigenvalue of the ODEs. To get the eigenvalue numerically, we first give a guess value of L and re-iterate the integration until the mass at the core, m(R c ), matches the actual mass M c . Note that, once the luminosity is found, the location of radiative-convective boundary (RCB) can be specified accordingly.

Envelopes without Heat Transport by Turbulent Mixing
For the convenience of comparison, we first consider a fiducial model, i.e., an envelope without turbulence (ν turb = 0). In Figure 1, we show the radial profiles of pressure, temperature, and density of the envelope for a 5M ⊕ core with increasing envelope mass during atmospheric growth. The green, cyan and yellow curves denote the envelope mass fraction (EMF) = 0.1, 0.4, 0.8, respectively. The thicker and thinner parts stand for the convective and radiative region, respectively. The boundaries of the thicker and thinner part are the radiative-convective boundaries (RCBs). The convective region is adiabatic. The radiative region connects the lower entropy interior to the higher entropy exterior. In Figure 1, we note that the pressure in the convection zone increases with envelope mass, but the temperature only varies slightly. Since the entropy is ∝ ln(T 1/∇ ad /P ), it is clear that, with increasing envelope mass, the steady-state envelopes evolve in order of decreasing entropy (Marleau & Cumming 2014). This is consistent with the cooling process that the envelope experiences, which allows the atmosphere to accrete more gas. Lee et al. (2014) found that, for dusty atmosphere, the location of RCBs lies at an roughly fixed temperature where H 2 dissociates (∼ 2500K). In Figure 1, the RCB lies at the bottom of the outermost radiative region and the temperatures at the RCBs are no longer 2500K. This is because we adopt a grain-free atmosphere due to efficient grain coagulation (Ormel 2014). According to the middle panel of Figure 1, we find that grain-free atmosphere behaves differently from grain-rich atmosphere. The outer radiative region is nearly isothermal, which implies that T RCB ∼ T out . Such features have also been identified in Lee & Chiang (2015 and Inamdar & Schlichting (2015), which can be readily understood in terms of the following relation (Rafikov 2006;Piso & Youdin 2014) The term on the right hand side of this equation is around the order of unity. This explains why the temperature at RCB, T RCB ∼ T out . We would stress that, the above relation is only valid for atmosphere without turbulence. When heat transport by turbulent mixing is taken into account, the RCB is pushed inwards, and the temperature at the RCB (T RCB ) becomes higher.
At the early stage of accretion, the envelope mass is small and the envelope can be well treated as non-self-gravitating. In this case, simple analytic results can be derived (Rafikov 2006;Piso & Youdin 2014). Though the envelope we consider in this paper is self-gravitating, these analytical results are still very instructive to understand atmospheric evolution and interpret our numerical results. How the position of the RCBs varies with envelope mass can be understood with the following relations (Piso & Youdin 2014), where ξ is a variable on the order of unity and P M is the characteristic pressure that is related to the core mass (Piso & Youdin 2014). In the early stage of planet accretion, with the increase of envelope mass, the pressure at RCB would increase as well. Accordingly, the cooling luminosity would be reduced. When the self-gravity becomes important, the above relations no longer hold. The stronger luminosity is necessary to support the more massive envelope. With the increase of luminosity, the RCB would be shifted outward as shown in Figure 1 (Ginzburg & Sari 2015).

Envelopes with Heat Transport by Turbulent Mixing
In this section, we explore how turbulence (ν turb = 0) changes the structure of the planet envelope. The most interesting feature is that the turbulence would push the RCBs inwards and diminish the cooling luminosity. In Figure 2, we show the planet thermal profiles for envelope mass fraction, M atm /M c = 0.2, 0.4 and 0.8. The core mass M c = 5M ⊕ . In Figure  2, we find that the difference is that a pseudo-adiabatic region appears. Explicitly, we point out the location of the pseudo-adiabatic region in the middle panel of Figure 2. In such regions, the temperature gradient is very close to adiabatic gradient, but still smaller than adiabatic gradient (see Figure 3).
From middle panel of Figure 1, we see that, when the heat transport by turbulent diffusion is not included, the RCB lies around the isothermal radiative region, T RCB ∼ T out . When turbulent diffusion is included, the temperature gradient would deviate from the isothermal approximation, which is most obvious by comparing middle panels of Figure 1 and Figure 2. We can identify from Figure 3, that the temperature gradient near RCBs is approaching ∇ ad and clearly deviates from isothermal temperature gradient. Due to  Fig. 1.-Thermal profiles around a planet core with mass M c = 5M ⊕ at 0.1 AU. Turbulence is not included, ν turb = 0. The pressure, temperature, and density are shown in the upper, middle, and lower panels, respectively. In each panel, the green, cyan, yellow lines stand for M atm /M c = 0.1, 0.4, 0.8, respectively. With the increase of envelope mass, the pressure at RCBs always increases. However, the position of RCBs inside the planet first decreases and then increases. This non-monotonic behavior is due to the effects of self-gravity (Piso & Youdin 2014). Note in particular that no pseudo-adiabatic region appears in the envelope (cf. Figure 2). this temperature gradient deviation, a pseudo-adiabatic region appears. As a result, the temperature at RCB becomes higher and RCBs would penetrate deeper inside the envelope.
To better understand the effects of heat transport by turbulent mixing, we compare the profiles of planet envelope with and without turbulence. The results are shown in Figure 3 as red solid and blue dashed lines, respectively. The upper panel shows the global variation of temperature with pressure within the envelope. In this panel, the difference between the two cases with and without turbulence is hardly discernible. The middle panel shows again the variation of temperature with pressure but focuses on the localized region around the radiative-convective transition region. It shows that the turbulent mixing smoothes the transition toward the adiabat. There would appear a pseudo-adiabatic region above the actual adiabatic region. This pseudo-adiabatic region pushes the RCB inward to higher pressure. Turbulent mixing leads to a more gradual approach to adiabat.
The turbulent diffusion in stably stratified region provides heating, instead of cooling so it is natural to expect that with turbulent diffusion taken into account, the total cooling rate of envelope will decrease and KH contraction timescale would be prolonged (see Figure  4 for details).

Onset of Gas Runaway and Cooling Luminosity Variations
Since we are interested in the planet accretion history, it is necessary to investigate the luminosity with increasing envelope mass. In the deep atmosphere, heat is advected by convective eddies. Near the surface, this could be achieved by diffusion. The surface temperature gradients would become shallower and a radiative region shows up. The variations of luminosity with envelope mass is shown in Figure 4. With the accumulation of envelope mass, the luminosity reaches a minimum. Beyond this minimum, the luminosity L increases. As a result, the planet begins to cool at a very short timescale and the envelope mass would grow super-linearly after this epoch. Physically, it is natural to adopt the epoch when the minimum L is reached as the onset of gas runaway, t run .
On the right hand side of luminosity minimum, the luminosity-mass relation is relatively easy to understand. At this late stage of mass growth, the self-gravity of envelope appears to be prominent, and greater luminosity is necessary to support stronger gravity. However, on the left hand side of the luminosity minimum, the mass of envelope is small and the planet is at the its early stage of mass growth. At this early stage (envelope's self-gravity can be ignored), the luminosity diminishes with a thicker radiative outer layer and more massive envelope. This reduction in cooling luminosity is intimately related to the shift of RCBs. The temperature at RCBs becomes higher when heat transport by tidally-forced turbulent mixing is taken into account. To identify their differences, we show the two profiles near the radiative-convective transition region. The red curve shows a more gradual transition from the radiative to adiabatic region. The region between the blue dot and red dot is the pseudo-adiabatic region. Bottom panel: The ratio of temperature gradient to adiabatic gradient. The region with ∇/∇ ad = 1 is the convection zone. The region with ∇/∇ ad < 1 is the radiative zone. In the pseudo-adiabatic region, ∇ is very close to ∇ ad , but still smaller than ∇ ad . The RCBs shift inwards when heat transport by turbulent mixing is taken into account. The RCBs penetrate deeper with stronger turbulent mixing.
When the envelope self-gravity can be ignored, the luminosity at RCB can be written as (Piso & Youdin 2014) where M RCB and L disk reads The above equations can be written in terms of known properties if the envelope mass is centrally concentrated (see, e.g., Lee & Chiang 2015). This central concentration is physically expected since in deeper layers where temperatures rise above ∼ 2500K, hydrogen molecules dissociate. As energy is spent on dissociating H 2 molecules rather than heating up the gas, the adiabatic index drops below 4/3, to approach 1. The upshot is that both the densities at the RCB and the radiative luminosity can be written in terms of core properties and the temperature at the RCB.
As RCB deepens, the RCB becomes even more optically thick so it is harder to radiate energy away; as a result, the envelope cools more slowly.
In Figure 4, we stress that two important aspects of thermal evolution during the planet accretion would be affected by turbulent mixing. The first is that it influences the luminosity. In Figure 4, we know that when the turbulent diffusivity (ν turb ) is enhanced, the cooling luminosity would be reduced globally. That is, for any particular value of envelope mass, the cooling luminosity for an envelope with turbulence is always below that without turbulence. When the turbulence is stronger, the luminosity becomes even smaller. The second is that it changes the EMF at which the gas runaway occurs. In Figure 4, our calculations show that, when the turbulence becomes stronger, the onset of gas runaway takes place at higher envelope mass fraction (EMF).

Quasi-Static KH Evolution and Critical Turbulent Diffusivity
Since we ignore the accretion luminosity from the planetesimals, the gravitational KH contraction is the only source for the cooling. The gas accretion is regulated by the KH timescale. Our time evolution model can follow the envelope mass growth up to the very early epoch of runaway growth around the crossover mass. Fortunately, Pollack et al. (1996) found that the timescale spent in the runaway accretion stage is orders of magnitude smaller than the KH timescale. The mass growth timescale is actually dominated by the KH stage. When the envelope mass is small, the increase of envelope mass causes the luminosity to decrease. When the envelope mass is sufficiently large, the self-gravity of gas envelope become important, and bigger luminosity L is necessary to support stronger gravity. We choose the luminosity minimum as the epoch when the runaway accretion sets in. We note that two important aspects of thermal evolution during the planet accretion would be affected. With the enhanced turbulence, the cooling luminosity is reduced globally. When the turbulence becomes stronger, the onset of gas runaway occurs at a higher envelope mass fraction.
For this reason, our model can get rather accurate estimation of mass growth timescale of an accreting planet. In this section, we will explore how the turbulent mixing affect the KH contraction timescale. For strong turbulent diffusion, the heat transport may even inflate the planet (Youdin & Mitchell 2010). We are not interested in planet inflation induced by strong turbulence in this paper. We find that even weak turbulence can already play an essential role to delay the KH contraction.

Time evolution: Temporally Connecting Snapshots
In the previous section, we have obtained snapshots of envelope structure for different envelope masses. To estimate the accretion timescale, we need to connect them temporally in order of increasing mass. The gas accretion history can be followed by the cooling process (Piso & Youdin 2014). Detailed estimation shows the luminosity generated in the radiative region can be safely ignored (Lee et al. 2014). It is physically valid to assume the luminosity of the envelope is generated in the convective zone and the luminosity can be treated as constant in the outer radiative zone (Piso & Youdin 2014). This would greatly simplify our evolutionary calculations. Under such circumstances, we only need to solve a set of ordinary differential equations and connect the solutions in time. Lee & Chiang (2015) shows that it is physically valid to omit planetesimal heating during the gas accretion of super-Earths. When there is no planetesimal accretion to power the gas envelope, the time interval between two adjacent hydrostatic snapshots is the time it spends to cool between them. In addition to internal energy variations, gas accretion and envelope contraction also bring about changes to the global energy budget. Specifically, the time interval between two steady state solutions can be written as (Piso & Youdin 2014) Note that the symbol ∆ designates the difference between the two adjacent states and the bracket denotes the average of them. The total energy E consists of the internal energy and the gravitational potential energy, which reads where u is the specific internal energy, u = c v T . The second term in equation (26) stands for contribution from gas accretion. The specific energy of the accreting gas is e = −GM r /r + u. The third term in equation (26) accounts for P dV work done by the envelope contraction. All terms are calculated at the RCB. Note in particular that the volume difference between two adjacent snapshots are performed at fixed mass. We choose the fixed mass as the average of the masses at the RCB (Piso & Youdin 2014).
In the upper panel of Figure 5, we shown the planet mass growth history for different turbulent diffusivity. In our fiducial model without turbulence, t run ∼ 4.04 Myrs. Beyond this epoch, the gas runaway occurs. The gas runaway is due to the fast increase of L beyond t run , which leads to a rapid cooling process on a shorter timescale. The most intriguing feature is that the runaway time is delayed and accretion timescale is prolonged when heat transport by tidally-forced turbulent mixing is taken into account. For instance, when ν turb = 0.0016, 0.005, 0.016, the runaway time, t run = 10, 18.4, 48.3 Myr, respectively. The stronger the turbulence, the longer the gas runaway timescale.
In our calculations, we find that a small value of ν turb , on the order of 10 −3 , can already appreciably affect the cooling timescale of super-Earths. Since ν turb is dimensionless, it is better to recover its physical value according to Equation (17). Typically, luminosities for super-Earths are L ∼ 10 26 erg/s, M ⊕ = 5.97 × 10 27 g, and ρ 0 = 6 × 10 −6 g cm −3 . Then the term, L/(GM ⊕ ρ 0 ), defined in Equation (17) is approximately ∼ 4.2 × 10 10 cm 2 s −1 . For the dimensionless diffusivity ν turb = 0.0016, the physical diffusivity is approximately µ turb ∼ 4.2 × 10 7 cm 2 s −1 . For even larger ν turb , the K-H contraction timescale can be enhanced by orders of magnitude. According to Figure 5, it is evident that the turbulent diffusivity on the order ∼ 10 7 −10 8 cm 2 s −1 can already enhance the runaway timescale by an oder of magnitude. The pressure scale height inside the planet is H p ∼ 10 9 cm and the sound speed is c s ∼ 10 5 cm s −1 . We can get a physical sense how large the turbulent diffusivity is by estimating the dimensionless parameter ζ in Equation (17). In our calculation, the parameter ζ is pretty small, on the order of 10 −7 ∼ 10 −6 . This means that the turbulent diffusion necessary to prolong the cooling timescale needs not to be very strong.

Critical Turbulence Diffusivity ν critical and Super-Earth Formation
A gas giant would be formed if the protoplanetary disk is still full of gas when the planet enters the runaway accretion stage. However, if the runway time t run is longer than the disk lifetime t disk , the disk gas is depleted and the planet is unable to accrete sufficient gas to become a gas giant, then a super-Earth may be formed. Two timescales, t run and t disk , determine the ultimate destiny of the planet, i.e., whether the planet becomes a super-Earth or a gas giant. If t run < t disk , gas runaway occurs within the lifetime of disk. The planet would get inflated by the runaway gas accretion and become a gas giant. On the contrary, if t run > t disk , the disk disperses before the gas runaway takes place. Because there is not enough gas material for the planet to accrete, the planet is unable to become a gas giant. The accretion history for ν turb = 0, 0.0016, 0.005, and 0.016 is shown as cyan dot-dashed, blue solid, green dotted, red dashed lines, respectively. The initial time for the accretion is estimated as t 0 = |E|/L. The slightly different starting time is due to the luminosity decrease by the inclusion of turbulence (see Figure 4). The initial EMF is around 6%, where the planet is nearly fully convective. Different color dots in the upper panel denote the epoch, t run , when the gas runaway takes place. The runaway time is t run = 4.04, 10, 18.4, and 48.3 Myrs, respectively. The solid blue curve shows the critical solution, where t run = t disk . The critical diffusivity for M core = 5M ⊕ is ν critical ∼ 1.6 × 10 −3 if t disk = 10 Myrs. Lower panel : The critical ν critical for various core mass. For higher core mass, the critical ν critical is higher. We note that a weak turbulence with small diffusivity, µ turb ∼ 10 7 − 10 8 cm 2 s −1 , can already enhance the runaway timescale and delay the gas runaway.
Usually the disk life is about 5 − 10 Myr. To be specific, we take the disk lifetime as t disk = 10 Myrs throughout this paper.
In the upper panel of Figure 5, the core mass fixed at M core = 5M ⊕ . We find that there exists a critical diffusivity ν critical = 1.6 × 10 −3 . When ν turb > ν critical , the K-H contraction timescale becomes longer than the disk lifetime and the core would not be able to experience the gas runaway. In this case, the formation of gas giants can be avoided and the formation of super-Earths becomes viable. In the lower panel of Fig. 5, we show the variations of ν critical with M core . The critical diffusivity becomes larger when the core mass increases. Specifically, for a 10 Earth mass core, the critical dimensionless diffusivity is approximately ν critical = 3.2 × 10 −2 . The actual diffusivity is about ∼ 10 9 cm 2 s −1 .

Variations of ν critical with Planet Location in PPDs
Observationally, the Kepler statistics show that ∼20% of Sun-like stars harbors super-Earths at distance of 0.05-0.3 AU. By contrast, the occurrence rate for hot Jupiters inside ∼ 0.1 AU is only 1%. To explain these observational features, we consider how the turbulence affects the thermal evolution for planets at different locations in PPDs. The turbulent mixing considered in this paper is driven by the tides raised by the host star. We believe that the tidally-induced turbulent mixing inside the planet would become weaker when the planet is farther away from the host star. Lee et al. (2014) found that, for dusty disk, the runaway timescale is independent of the orbital location. However, since dust can not persist in the envelope due to coagulation and sedimentation (Ormel 2014;Mordasini 2014), the runaway timescale is no longer independent of the orbital location. In the upper panel of Figure 6, we show the accretion history for planets at three different locations. The core mass is M c = 6M ⊕ . The blue solid, green dot-dashed, and red dashed curves designate the temporal variations of envelope mass for a = 0.1AU, 1AU, and 5AU, respectively. The turbulent diffusivity is ν turb = 0.013. The gas runaway occurs at t run = 33.1, 3.2, and 1.7 Myrs. It is clear that gas accretion onto cores is hastened for planets that are farther away from the central star. This behaviour can be understood from the decrease in opacity at the RCB which makes the envelope more transparent, enhancing the rate of cooling (Lee & Chiang 2015;Inamdar & Schlichting 2015). The planet at a = 1AU and 5 AU would become a gas giant due to runaway accretion (t run < t disk ). However, the planet in the inner region a = 0.1 AU would become a super-Earth (t run > t disk ). The fact that atmospheres cool more rapidly at large distances as dust-free worlds has been used to explain the presence of extremely puffy, low mass planets (Inamdar & Schlichting 2015;Lee & Chiang 2016).
We explore the critical diffusivity, ν critical , for planets at different locations inside the minimum mass extrasolar nebula (MMEN). The results are shown in the lower panel of Figure 6. It shows that the critical diffusivity increases with the semi-major axis. When the planet is farther from the central star, ν critical becomes larger. This means that the more distant planet requires stronger turbulence to lengthen the KH timescale and avoid gas runaway. For tidally-induced forcing, we believe that the turbulent diffusion ν turb is determined by the tides inside the planet raised by host star. The tides become weaker if the planet is farther away from the host star. Our proposed mechanism can naturally explain the formation of close-in super-Earth, while still ensuring the gas giant formation at larger orbital distance. When the planet is near the host star, tidally-forced turbulent mixing is stronger and ν turb would be larger. According to Figure 6, the required threshold ν critical is smaller, As a result, the inequality ν turb > ν critical can be more readily to satisfy and formation of super-Earth becomes possible. On the contrary, when the planet is far from the host star, ν turb becomes smaller as the stirring by tides becomes weaker. The required threshold ν critical becomes larger. The threshold to avoid gas runaway is more difficult to satisfy. This indicates that, in the in-situ planet formation scenario, it is more readily to form close-in super-Earths and gas giants are more prone to appear in the outer region of PPDs. The above implication is consistent with occurrence rate inferred from observations.

Mass Loss Mechanisms
Observation shows that super-Earths possess hydrogen and helium envelopes containing only several percent of the planet's mass. However, we can see in Figure 5 that the planets accrete very massive gas envelopes. The planet core with ν turb = 0 reaches an envelope mass fraction (EMF) of ∼ 0.8 at the epoch of gas runaway. The envelope mass is considerably higher than the mass inferred from observations. These primordial super-Earths may experience significant mass loss during the post-formation evolution.
How super-Earths lose their mass still remains an open question. Here we briefly mention some possible ways to lose the envelope mass. The first possibility is that close-in planets are exposed to intense XUV (extreme UV and X-ray) irradiation from their host stars. Photoevaporation can significantly modify the structure of their atmosphere. Over the timescale of ∼ 100 Myrs, X-rays from host stars can photoevaporate the super-Earth envelopes from initial EMF ∼ 1 to EMF of ∼ 0.01 − 0.1, which may naturally explain the differences between the theoretical predictions and observational facts (e.g., Murray-Clay et al. 2009;Owen & Wu 2013;Owen & Wu 2017;Gaudi et al. 2017). Giant impact is the second possible mechanism to explain the mass loss, which is ex- The turbulent diffusivity is ν turb = 0.01. The blue solid, green dot-dashed, red dashed lines denote mass growth history for planets at 0.1 AU, 1 AU, and 5 AU, respectively. The critical mass ratio at the epoch of runaway decreases for more distant planet. The runaway time for the three different cases are 33.1, 3.2, and 1.7 Myr, respectively. It is expected that for more distant planets, larger turbulent diffusivity is required to prevent runaway gas accretion within t disk ∼ 10 Myrs. Lower panel : The critical diffusivity, ν critical , for different orbital locations, required to prevent gas runaway for disk lifetime t disk ∼ 10 Myrs. Beyond ν critical , the KH timescale is longer than the disk lifetime. The formation of super-Earths becomes possible.
pected to be common because they are needed to provide long-term orbital stability of planetary systems (Cossou et al. 2014). Hydrodynamical simulations show that a single collision between similarly sized exoplanets can easily reduce the envelope-to-core-mass ratio by a factor of two. Super-Earths' asymptotic mass can be achieved by one or two giant impacts. Under certain circumstances, almost 90% of the gas envelope can be lost during impact process (Liu et al. 2015;Inamdar & Schlichting 2016).
Mass transfer between the close-in planet and host star via Roche lobe represent the third way to reduce the planet mass (Valsecchi et al. 2015;Jia & Spruit 2017;Jackson et al. 2017). Tidal dissipation can drive orbits of these primordial super-Earths to decay toward the Roche limit. The mass transfer is quite rapid, potentially leading to complete removal of the gaseous envelope in a few Gyr, and leaving behind a super-Earth. Many gaseous exoplanets in short-period orbits are on the verge or are in the process of Roche-lobe overflow (RLO). The coupled processes of orbital evolution and RLO likely shape the observed distribution of close-in exoplanets and may even be responsible for producing some of the short-period rocky planets. But recent calculations by Dosopoulou et al. (2017) challenged this idea by claiming that, for high eccentric planets or retrograde planets, self-accretion by the planet would slow down the mass loss rate via Roche lobe overflow.
Super-Earth envelope mass fractions range just 1-10% and more typically just ∼1% (see Rogers & Seager 2010, Lopez & Fortney 2014, Wolfgang & Lopez 2015. The mechanism discussed in this paper overpredicts the envelope mass fraction of super-Earths, often beyond 80%. Photoevaporation, even around Sun-like stars, are only effective out to 10 days and many super-Earths lie beyond this (see, e.g., Figure 8 of Owen & Wu 2013). Removal of >90% of the envelope by giant impact requires impact velocity that exceeds the escape velocity (see, e.g., Figure 3 of Inamdar & Schlichting 2016). Finally, Roche lobe overflow only works within 2 stellar radii where the Roche radius is. Lee & Chiang (2016) proposed that the late-time formation of cores ensures that super-Earth cores accrete a few percent envelope mass fraction, in agreement with the observations. There is a clear difference in the expected final envelope mass fraction between their work and ours.
Very recent works have revealed that planetary envelopes embedded within PPDs may not be in hydrostatic balance, which slows down envelope growth. It is possible for a steady state gas flow enters through the poles and exits in the disc mid-plane (Lambrechts & Lega 2017). In the presence of a magnetic field and weakly ionizing winds, ohmic energy is dissipated more readily for lower-mass planets. Ohmic dissipation would make super-Earths more vulnerable to atmospheric evaporation (Pu & Valencia 2017). These findings may offer new explanations for the typical low-mass envelopes around the cores of Super-Earths. In addition, we also note that the turbulent diffusion mechanism may be still operating in the late core assembly scenario. In the late core assembly scenario without turbulent diffusion, the asymptotic EMF is about 3-5% (Lee & Chiang 2016). When turbulent diffusion is taken into account, the EMF can be further reduced to 1%.

Summary and Conclusion
In this paper, we propose a new mechanism to avoid gas runaway for planet cores within the lifetime of disks. The mechanism proposed in this paper is not subject the κ or µ catastrophe (Lee & Chiang 2015). Tidal heating (Ginzburg & Sari 2017) requires orbital eccentricity be continuously pumped up during super-Earth formation. Our mechanism does not depend on the orbital eccentricity of super-Earth. Incorporating this model into a population synthesis model may better constrain our understanding of the exoplanet formation (Ida &Lin 2004;Jin & Mordasini 2017).
We have explored the effects of heat transport induced by tidal stirring on the thermal structure of stably stratified, radiative layers of super-Earths, focusing on their influences on the KH timescale. When we take turbulent stirring into account, pseudo-adiabatic regions would show up within the radiative zone. This may push the RCBs inwards. The temperature, pressure at RCBs becomes higher and the cooling luminosity would be reduced. As a result, the KH timescale would be enhanced. We find that there exist a critical turbulent diffusivity ν critical . When ν turb > ν critical , the runaway time is greater than the disk lifetime (t run > t disk ). Under such circumstances, the onset of the planet gas runaway lags behind the disk gas depletion. Since the planet has not enough gas to accumulate, it can no longer grow into a gas giant and become a super-Earth instead. In addition, we also investigate the variations of ν critical with planet's semi-major axis in MMEN. Our calculations show that the condition for turbulence-induced formation of super-Earths is more readily satisfied in the inner disk region, but is harder to satisfy in the outer disk region. The occurrence rate of super-Earths and gas giant is consistent our calculations.
issues will be addressed in a further study.
A limitation of this work is that the turbulence strength is not specified from first principle. As a compromise, we parameterize the turbulence diffusion as a free parameter. We try to constrain the turbulence strength in terms of the planet thermal evolution. Interestingly, we find the turbulence in the radiative region have substantial effects on the planet accretion history. How turbulence is initiated during the planet formation and how strong the turbulent diffusion is involve very complicated physical processes, which are worth further investigations.
Realistic opacities and EOS have influential effects on the planetary thermal structure and the core accretion process (e.g. Stevenson 1982;Ikoma et al. 2000;Rafikov 2006), especially for timescale of the KH timescale (Lee et al. 2014;Piso & Youdin 2014). Our simple prescription of opacity needs to be improved. Guillot et al. (1994) showed that an convective layer lies between two adjacent radiative region due to the opacity window near ∼ 2000K. A relevant caveat is the existence of radiative zones sandwiched inside convective interior. Such radiative windows are ignored in our two-layer models. It would be interesting to consider how a downward turbulent heat flux would interact with such a sandwiched region. In summary, how super-Earth envelope cooling history responds to more realistic opacities and EOS needs to be further investigated. Calculation with realistic EOS and opacity are underway and will be reported elsewhere.
We have found that the epoch of runaway accretion can be effectively delayed by the turbulent diffusion within the stably stratified region. But we should be cautious that the envelope mass fraction predicted by this mechanism is not fully consistent with observations. The envelope mass fraction for planet embedded within the gas-rich MMEN is greater than 80%, much higher than the typical super-Earth envelope. It is difficult for the turbulent diffusion alone to make the envelope mass fraction be consistent with observations. Additional physical process, such as giant impact, photo-evaporation, Roche-lobe overflow may be operating to reduce the envelope mass fraction during the formation of super-Earth. But these mass loss processes either operate on distances shorter than most super-Earths or are applicable under certain circumstances. A promising mechanism for super-Earth formation is the late core assembly within the transitional PPDs. In this scenario, with the reduction of the PPD mass density, the envelope mass fraction can be as low as 3-5% (Lee & Chiang 2016). We note that the turbulent diffusion may be still working in the late core assembly scenario. How turbulent diffusion affect the envelope mass fraction within transitional PPDs is an interesting issue worth further investigation.
We thank the anonymous referee for the thoughtful comments that greatly improve this paper. Discussions about heat transport inside planet interior with Yanqin Wu and Re'em Sari are highly appreciated. This work has been supported by National Natural Science Foundation of China (Grants 11373064, 11521303, 11733010), Yunnan Natural Science Foundation (Grant 2014HB048) and Yunnan Province (2017HC018).

A. Temperature Gradient with Turbulence
When the turbulent flux is taken into account, the luminosity can be written as Note that in the above equation g = GM r /r 2 and ∇ = d ln T d ln P . This equation can be expressed as dT dr = L + 4πr 2 µ turb ρ GMr If we take the pressure P as the independent variables, the temperature gradient becomes where η = 4πµ turb GM r ρ/L.