Suppression of atmospheric recycling of planets embedded in a protoplanetary disc by buoyancy barrier
MMNRAS , 1–15 (2017) Preprint 6 June 2018 Compiled using MNRAS L A TEX style file v3.0
Suppression of atmospheric recycling of planets embeddedin a protoplanetary disc by buoyancy barrier
Hiroyuki Kurokawa, (cid:63) Takayuki Tanigawa Earth-Life Science Institute, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro, Tokyo 152-8550, Japan National Institute of Technology, Ichinoseki College, Takanashi, Hagisho, Ichinoseki-shi, Iwate 021-8511, Japan
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The ubiquity of super-Earths poses a problem for planet formation theory to explainhow they avoided becoming gas giants. Rapid recycling of the envelope gas of plan-ets embedded in a protoplanetary disc has been proposed to delay the cooling andfollowing accretion of disc gas. We compare isothermal and non-isothermal 3D hy-drodynamical simulations of the gas flow past a planet to investigate the influenceon the feasibility of the recycling mechanism. Radiative cooling is implemented byusing the β cooling model. We find that, in either case, gas enters the Bondi sphereat high latitudes and leaves through the midplane regions, or vice versa when disc gasrotates sub-Keplerian. However, in contrast to the isothermal case where the recyclingflow reaches the deeper part of the envelope, the inflow is inhibited from reaching thedeep envelope in the non-isothermal case. Once the atmosphere starts cooling, buoyantforce prevents the high-entropy disc gas from intruding the low-entropy atmosphere.We suggest that the buoyancy barrier isolates the lower envelope from the recyclingand allows further cooling, which may lead runaway gas accretion onto the core. Key words: hydrodynamics – planets and satellites: atmospheres – planets andsatellites: formation – protoplanetary discs.
The
Kepler mission has revealed that about half of Sun-likestars harbour super-Earths, here defined as planets whoseradii are 1–4 R ⊕ , with periods <
85 days (e.g., Fressin etal. 2013). Those super-Earths have the masses of 2–20 M ⊕ (e.g., Weiss & Marcy 2014) and their mass-radius relationsindicate that the gas-to-core mass ratios are ∼ ∼
1% (e.g., Fressin etal. 2013).The ubiquity of super-Earths poses a problem for planetformation theory: how did they avoid becoming gas giants?A solid core embedded in a protoplanetary disc captures thesurrounding disc gas within the Bondi or Hill spheres. Ra-diative cooling and contraction of the envelope gas allowsfurther accretion of disc gas (e.g., Lee et al. 2014). If themass of the core exceeds 5–20 M ⊕ , it undergoes runaway gasaccretion to become a gas giant (e.g., Mizuno 1980; Pollacket al. 1996; Ikoma et al. 2000). The critical core mass be-comes even smaller if the envelope was polluted with heavyelements (Hori & Ikoma 2011; Venturini et al. 2015).Late-stage core formation has been proposed to solve (cid:63) E-mail: [email protected] (HK) the problem of ubiquitous super-Earths. Given the shorttime-scale of the runaway gas accretion ( ∼ M ⊕ core), it has been proposed that super-Earth cores co-agulated just as disc gas was about to disappear completelyso that they did not have sufficient time to undergo run-away accretion (Ikoma & Hori 2012; Lee et al. 2014). Thefinal assembly during disc dispersal is the expected result inconventional planet formation theory due to mutual gravi-tational interactions of proto-cores (e.g., Kominami & Ida2002; Inamdar & Schlichting 2015).The feasibility of the late-stage core formation sce-nario depends on the planet formation regimes. The clas-sical planet formation theory assumed that the proto-coresare formed by planetesimal accretion. The isolation mass ofthe planetesimal accretion is small at ∼ . (Kokubo &Ida 2000), which is enough to avoid the runaway gas accre-tion before the final assembly. On the other hand, in thepebble accretion theory (Ormel & Klahr 2010; Lambrechts& Johansen 2012), where proto-cores accrete pebbles drift-ing from the outer regions of the disc, proto-cores can growup to a mass heavy enough to carve a gap in the pebble disc(Lambrechts et al. 2014). The inward drift of pebbles hasbeen suggested observationally from the disc sizes, whichare found to be more compact in the continuum (pebbles)than in the lines (gas) (Ormel 2017, and references therein), © a r X i v : . [ a s t r o - ph . E P ] J un H. Kurokawa and T. Tanigawa and from dust rings, which might be formed by the drift andsintering of icy dust (Okuzumi et al. 2016).Other mechanisms have been suggested to inhibit super-Earths from becoming gas giants by delaying the cooling andcontraction of their proto-envelopes. When super-Earths areon eccentric orbits, stellar tidal heating would suppress thecooling (Ginzburg & Sari 2017). Yu (2017) has proposedthat the tidally-forced turbulent diffusion may affect theheat transport inside the planet and may explain the factthat super-Earths are commonly found in the inner proto-planetary disc.Another mechanism to explain the ubiquity of super-Earths utilises a hydrodynamical explanation – a rapid recy-cling of the atmospheric gas embedded in a protoplanetarydisc may prevent the envelope from cooling (Ormel et al.2015b). They conducted isothermal hydrodynamical simu-lations of flow past a planet embedded in a protoplanetarydisc. The atmosphere inside the Bondi radius was shown tobe an open system where disc gas enters from high latitude(inflow) and leaves through the midplane region (outflow) incases where shear flow was adopted. They also found thatthe topography changes when the rotation of disc gas is sub-Keplerian, but the recycling was observed in any cases. Theyargued that the recycling is faster than the cooling of theenvelope gas for close-in super-Earths, and so that furtheraccretion of disc gas is prevented.The recycling mechanism works only if the high-entropydisc gas can replace the low-entropy atmospheric gas. Thetreatment of the thermal structure is critical to assesswhether the recycling can halt the cooling and further accre-tion of disc gas. Rapid recycling has been observed in isother-mal simulations (Ormel et al. 2015b; Fung et al. 2015). Adi-abatic simulations have also found no dynamical boundarybetween the atmosphere and disc gas (Popovas et al. 2018).In contrast, non-isothermal, radiative-hydrodynamical sim-ulations have found a bound atmospheric region around aplanet below the recycling region (D’Angelo & Bodenheimer2013; Cimerman et al. 2017; Lambrechts & Lega 2017). Whythe efficiency of recycling differs between isothermal andnon-isothermal models is poorly understood, and is the focusof this study.We conduct a comparison between isothermal andnon-isothermal 3D hydrodynamical simulations of the flowaround a planet embedded in a protoplanetary disc to in-vestigate the effect of thermal profiles on the feasibility ofthe recycling mechanism. Section 2 presents the model. Sec-tion 3 shows the results. The implications for the formationof super-Earths are discussed in Section 4. We conclude inSection 5.
Flow around a planet embedded in a protoplanetary discwas calculated in this study. We followed the setup of the hy-drodynamical simulations of Ormel et al. (2015b) except forthree major differences. i) Whereas Ormel et al. (2015b) con-ducted isothermal simulations, we performed both isother-mal and non-isothermal simulations where the radiativecooling was implemented by using the β cooling model (sub-section 2.1). ii) Ormel et al. (2015b) implemented a modi-fied gravitational potential of a planet to avoid numerical difficulties found in their study. We conducted simulationsusing a standard gravitational potential (subsection 2.1). iii)Whereas Ormel et al. (2015b) omitted the vertical stratifi-cation of the circumstellar disc for simplification, we consid-ered the cases both with and without the stratification, forlower and higher mass planets, respectively (subsection 2.2). A compressible, inviscid fluid of an ideal gas was assumed inthis study. We used dimensionless units where velocities areexpressed in terms of the sound speed c s , times in terms ofthe reciprocal of the Keplerian frequency Ω − , and lengthsin terms of the disc scale height H ≡ c s / Ω of the backgroundgas. The dimensionless mass of a planet m follows the re-lation m = R Bondi / H , where R Bondi is the Bondi radius. TheHill radius is given by R Hill = ( m / ) / . As the gravity ofthe gas itself was neglected in this study, we can define an-other dimensionless unit to measure the densities. Here, thebackground gas density at the midplane ρ was set as unity.The scaling gives the planetary mass M p = mc /( G Ω ) ,where G is the gravitational constant. Assuming a solar-masshost star and a disc temperature profile T = ( a / ) − / (the minimum-mass solar nebula model, Weidenschilling1977; Hayashi et al. 1985), M p is given by, M p (cid:39) m (cid:18) a (cid:19) M ⊕ , (1)where a is the orbital radius.The continuity, Euler’s, and energy conservation equa-tions are given by, ∂ ρ∂ t + ∇ · ( ρ v ) = , (2) ∂∂ t ( ρ v ) + ∇ · ( ρ vv ) = ∇ p + ρ ( F cor + F tid + F hw + F p ) , (3) ∂ E ∂ t + ∇ · [( E + p ) v ] = ρ v · ( F cor + F tid + F hw + F p )− U ( ρ, T ) − U ( ρ, T ) β , (4)where v is the velocity, ρ is the density, p is the pressure.The internal and total energy densities U and E are givenby, U = p γ − , (5) E = U + ρ v , (6)where γ is the ratio of specific heat. We assumed γ = / .The last term in Equation 4 is the thermal relaxationapproximated by the β cooling model where the tempera-ture T relaxes towards the background temperature T withthe dimensionless time-scale β (e.g., Gammie 2001). The β cooling was implemented in order to replicate the influenceof radiative cooling on the recycling flow (Section 1). As wewill show in the following sections, our model is useful to un-derstand the mechanism which determines the efficiency ofatmospheric recycling. The advantages to using the β cool-ing model are, i) it is computationally less expensive than MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Name Cooling time Mass Domain size Headwind speed Injection time Termination time β m r out M hw t inj t end shear-is-m001 isothermal 0.01 0.5 0 0.5 15 shear-B4-m001 − shear-B3-m001 − shear-B2-m001 − shear-ad-m001 adiabatic 0.01 0.5 0 0.5 15 shear-is-m01 isothermal 0.1 5 0 0.5 15 shear-B2-m01 − shear-is-m1 isothermal 1 5 0 1 30 shear-B2-m1 − shear-B0-m1 headw-is-m001 isothermal 0.01 0.5 0.1 0.5 10 headw-B2-m001 − heado-is-m001 isothermal 0.01 0.5 0.1 0.5 10 heado-B2-m001 − Table 1.
List of the simulations. The first column gives the name where shear , headw , and heado depict shear flow, shear flow withheadwind, and headwind only, respectively. The second to seventh columns show values of the cooling time, planetary mass, computationaldomain size, headwind speed, injection time of gravity of the planet, and termination time of computation, respectively. radiative-transfer calculations, and ii) we can obtain insightsinto the effect of radiative cooling by changing β values ar-tificially. The disadvantages are, i) this approximation be-comes unrealistic for optically thick and non-linear regimes,and ii) we assume one β value even though the relaxationtime may vary from place to place: the low-density disc gasand dense atmosphere have different relaxation time scales.We will change β values as a free parameter in Section 2in order to understand the influence. In subsection 4.2, wewill discuss the implications for gas accretion onto planetsby assuming that the β value represents the relaxation timeof disc gas descending into the Bondi radius.The forces in the non-inertial frame read as follows: theCoriolis force F cor = − e z × v , the tidal force F tid = x e x − z e z (where e i is the unit vector in the i -direction), the globalpressure force due to the sub-Keplerian motion of the gas F hw = M hw e x (where M hw is the Mach number of the head-wind), and the gravitational force due to the planet F p . Herewe assumed that the centre of the planet is located at theorigin of the coordinates. For the lowest mass simulations( m = . ), the tidal force in the z -direction was neglectedfollowing Ormel et al. (2015b) (subsection 2.2).The gravity term implementing a smoothing in spacewas gradually inserted over time: F p = −∇ (cid:18) m (cid:113) r + r (cid:19) · (cid:26) − exp (cid:20) − (cid:18) tt inj (cid:19) (cid:21)(cid:27) , (7)where m is the dimensionless mass of the planet, r is thedistance from the centre of the planet, r sm is the smooth-ing length, and t inj is the injection time. We assumed r sm = − × m . We note that Equation 7 differs from the formused in Ormel et al. (2015b) because they adopted a mod-ified gravitational potential of a planet to avoid numericaldifficulties found in their study. We solved the equations described in subsection 2.1 by usingthe hydrodynamical simulation code Athena++ (White et al. 2016, Stone et al. in prep.), an extended re-write of theAthena astrophysical magnetohydrodynamics code (Stoneet al. 2008). We adopted a spherical grid with coordinates( r , θ, φ ), centred on the planet. In order to resolve flow insidethe Bondi radius, a logarithmic grid for the radial dimensionwas employed. We adopted polar mesh spacing proportionalto ( ψ + ) where ψ is the angle from the midplane, whichproduces high resolution near the midplane.The unperturbed solution to the shearing-sheet wasadopted for the initial and outer-boundary conditions: v = (cid:18) − x − M hw (cid:19) e y , ρ = exp (cid:18) − z (cid:19) , (8)where the density has vertical stratification due to the stellargravity (the tidal force in the z -direction in our non-inertialframe). In the case of the lowest mass planet ( m = . )where we omitted the z -direction tidal force, we adopted ρ = ρ = . The reflective boundary condition was assumedfor the inner boundary. However, after conducting severaltest runs, we found that there is unphysical energy flow atthe inner boundary when gas accretes, which affected theentire calculated domain. This seems to be caused by thereflective boundary condition as it is precisely reflective onlyfor the (uniform) Cartesian grid. In order to prevent this un-physical energy flux at the inner boundary from affecting theentire flow, we introduced an artificial cooling at the innerboundary where a small value of β ( = − ) was adopted.A summary of simulations is given in Table 1. We as-sumed the planetary mass m = . as a nominal value fol-lowing Ormel et al. (2015b), which allows us to compare ourresults to theirs. This dimensionless mass corresponds to a0.12 M ⊕ planet orbiting a solar-mass star at 1 AU (Equa-tion 1). The z -direction tidal force and vertical stratificationwere omitted in the nominal case. For higher mass cases( m = . and m = , corresponding to 1.2 M ⊕ and 12 M ⊕ planet at 1 AU), the density stratification was considered.We note that, because our local simulations do not allow aplanet to open a ringed gap in the disc, low mass planets( m < ) would be preferred. Nevertheless, high mass simu-lations ( m = ) were performed in order to check that ourconclusions do not depend on the assumed planetary mass. MNRAS000
List of the simulations. The first column gives the name where shear , headw , and heado depict shear flow, shear flow withheadwind, and headwind only, respectively. The second to seventh columns show values of the cooling time, planetary mass, computationaldomain size, headwind speed, injection time of gravity of the planet, and termination time of computation, respectively. radiative-transfer calculations, and ii) we can obtain insightsinto the effect of radiative cooling by changing β values ar-tificially. The disadvantages are, i) this approximation be-comes unrealistic for optically thick and non-linear regimes,and ii) we assume one β value even though the relaxationtime may vary from place to place: the low-density disc gasand dense atmosphere have different relaxation time scales.We will change β values as a free parameter in Section 2in order to understand the influence. In subsection 4.2, wewill discuss the implications for gas accretion onto planetsby assuming that the β value represents the relaxation timeof disc gas descending into the Bondi radius.The forces in the non-inertial frame read as follows: theCoriolis force F cor = − e z × v , the tidal force F tid = x e x − z e z (where e i is the unit vector in the i -direction), the globalpressure force due to the sub-Keplerian motion of the gas F hw = M hw e x (where M hw is the Mach number of the head-wind), and the gravitational force due to the planet F p . Herewe assumed that the centre of the planet is located at theorigin of the coordinates. For the lowest mass simulations( m = . ), the tidal force in the z -direction was neglectedfollowing Ormel et al. (2015b) (subsection 2.2).The gravity term implementing a smoothing in spacewas gradually inserted over time: F p = −∇ (cid:18) m (cid:113) r + r (cid:19) · (cid:26) − exp (cid:20) − (cid:18) tt inj (cid:19) (cid:21)(cid:27) , (7)where m is the dimensionless mass of the planet, r is thedistance from the centre of the planet, r sm is the smooth-ing length, and t inj is the injection time. We assumed r sm = − × m . We note that Equation 7 differs from the formused in Ormel et al. (2015b) because they adopted a mod-ified gravitational potential of a planet to avoid numericaldifficulties found in their study. We solved the equations described in subsection 2.1 by usingthe hydrodynamical simulation code Athena++ (White et al. 2016, Stone et al. in prep.), an extended re-write of theAthena astrophysical magnetohydrodynamics code (Stoneet al. 2008). We adopted a spherical grid with coordinates( r , θ, φ ), centred on the planet. In order to resolve flow insidethe Bondi radius, a logarithmic grid for the radial dimensionwas employed. We adopted polar mesh spacing proportionalto ( ψ + ) where ψ is the angle from the midplane, whichproduces high resolution near the midplane.The unperturbed solution to the shearing-sheet wasadopted for the initial and outer-boundary conditions: v = (cid:18) − x − M hw (cid:19) e y , ρ = exp (cid:18) − z (cid:19) , (8)where the density has vertical stratification due to the stellargravity (the tidal force in the z -direction in our non-inertialframe). In the case of the lowest mass planet ( m = . )where we omitted the z -direction tidal force, we adopted ρ = ρ = . The reflective boundary condition was assumedfor the inner boundary. However, after conducting severaltest runs, we found that there is unphysical energy flow atthe inner boundary when gas accretes, which affected theentire calculated domain. This seems to be caused by thereflective boundary condition as it is precisely reflective onlyfor the (uniform) Cartesian grid. In order to prevent this un-physical energy flux at the inner boundary from affecting theentire flow, we introduced an artificial cooling at the innerboundary where a small value of β ( = − ) was adopted.A summary of simulations is given in Table 1. We as-sumed the planetary mass m = . as a nominal value fol-lowing Ormel et al. (2015b), which allows us to compare ourresults to theirs. This dimensionless mass corresponds to a0.12 M ⊕ planet orbiting a solar-mass star at 1 AU (Equa-tion 1). The z -direction tidal force and vertical stratificationwere omitted in the nominal case. For higher mass cases( m = . and m = , corresponding to 1.2 M ⊕ and 12 M ⊕ planet at 1 AU), the density stratification was considered.We note that, because our local simulations do not allow aplanet to open a ringed gap in the disc, low mass planets( m < ) would be preferred. Nevertheless, high mass simu-lations ( m = ) were performed in order to check that ourconclusions do not depend on the assumed planetary mass. MNRAS000 , 1–15 (2017)
H. Kurokawa and T. Tanigawa
The inner domain boundary was r inn = − . The numericalresolution was [ r , θ, φ ] = [ , , ] . Our model showed differences in flow patterns of recyclingbetween isothermal and non-isothermal cases. We comparedthe time sequence in isothermal and non-isothermal cases inFigures 1 and 2 for the nominal mass case ( m = . ). Gasstarted to accrete onto a planetary core as the gravity wasgradually inserted with t inj = . (Figures 1a and 2a). At t = , the density has reached its steady-state value (Figures1b and 2b). It means that a hydrostatic atmosphere has beenbuilt at this point. To the first approximation, a comparisonbetween t = (Figures 1c and 2c) and t = (Figures 1dand 2d) shows that the velocity field has reached a steadystate in terms of its basic properties (the presence/absenceof atmospheric recycling; see following paragraphs), thoughthe streamlines show some time-varying features. The time-varying velocity field has also been reported in Ormel et al.(2015b).The isothermal simulation (Figure 1) showed outflowleaving the planet in the midplane as part of the recyclingflow, which entered the Bondi sphere from the polar regionsas shown by Ormel et al. (2015b). In contrast, the non-isothermal run (Figure 2) showed more circular streamlinesaround the planetary core within r < ∼ . R Bondi (rememberthat the dimensionless planetary mass m is the measure ofthe Bondi radius), meaning that the recycling of envelopegas is limited and envelope gas is isolated in this region. Weobserved streamlines heading outwards outside the isolatedenvelope but still within the Bondi sphere. The differencein flow patterns has been maintained until the end of thecalculations ( t = , Figures 1d and 2d).The non-isothermal simulations showed recycling flowssimilar to what has been reported in the isothermal case(Ormel et al. 2015b), but the recycling was observed onlyoutside the isolated inner envelope. We showed 3D stream-lines in Figure 3a. Fluid descended in polar regions as itcirculated outside the isolated inner envelope. Once the de-scending gas reached the midplane, it flowed outwards in themidplane. In contrast, the isothermal simulations showed therecycling inflow and outflow being connected to the plane-tary core as reported previously (Ormel et al. 2015b). The size of the isolated inner envelope where we observedcircular streamlines is positively correlated to the coolingtime-scale β . We compared the midplane flow patterns ofthe isothermal ( β → ), β = − , β = − , β = − , andadiabatic ( β → ∞ ) simulations at the same time well afterthe hydrostatic structure was developed ( t = ) in Figure4. Similar to the isothermal run (Figure 4a), the outflowstreamlines were launched from the vicinity of the core in thecase where β = − was assumed (Figure 4b). We note thatthe flow patterns were not completely the same in the twocases. The longer the cooling time-scale β we assumed, the larger region of circular streamlines we observed (Figures 4b-d). There exists a dynamical boundary: the inner envelope isisolated from the shear and horseshoe flow of disc gas (Figure4e). As in the β = − case (Figure 3a), the recycling flowwas observed outside the isolated inner envelope (but stillwithin the Bondi sphere) in all cases.The isolated inner envelope disappears in the slow cool-ing limit ( β → ∞ ). The adiabatic simulation did not showthe dynamical boundary between the envelope and disc gasin its stream lines (Figure 4f). The absence of isolated en-velope agrees with the result of the isothermal case, thoughthe adiabatic case showed a rather time-varying flow pat-tern. This is possibly because the smaller density of the en-velope in the adiabatic simulation allowed the disc gas flowto perturb the envelope more easily. As we will explain insubsection 3.3, a comparison of these simulations suggeststhat the presence or absence of the dynamical boundary iscontrolled by the entropy difference between the descendinginflow gas and the envelope, which exists only in the non-isothermal (finite β ) simulations. We note that the natureof our adiabatic simulation is consistent with the previousstudy (Popovas et al. 2018).In contrast to the topology of the flow, the density pro-files were almost the same in the simulations assuming dif-ferent values of β except for the adiabatic case. Because theflow was subsonic, the pressure structure was determined tobe in hydrostatic equilibrium with the gravity of the planet.Temperature was nearly uniform due to the efficient radia-tive cooling. Therefore, the density profiles did not dependon β .In order to measure the efficiency of atmospheric recy-cling more quantitatively, Figure 5 shows a comparison ofthe shell-averaged recycling and perturbation flux, which isdefined by, F out ( r ) ≡ ∫ ρ · | v r | · r d Ω π r , (9)where the integral denotes a surface integral over the shell.Because the flux includes both actual recycling and pertur-bation within the bound atmosphere, the value representsthe upper limit of the recycling flux. Initially ( t < ), theplanet accreted the disc gas until the hydrostatic equilibriumwas established, so that the flux within the Bondi sphere( r < − ) was high due to the large influx (Figure 5a). Afterthe accretion period, the flux reached a nearly steady-statevalue. The larger β we assumed, the smaller the recyclingflux we obtained within the Bondi sphere ( r < − , Fig-ure 5b). The flux in the case of β = − is smaller thanthat in the isothermal simulations by one or two orders ofmagnitude. It means that the envelope became more sta-ble against the recycling as we assumed slower (inefficient)cooling of the descending disc gas (see subsection 3.3). Thechange in the position of dynamical boundary between theatmosphere and disc gas found in Figure 4 is also visiblein Figure 5b: the position where the flux has a minimummoved to outer regions as we assumed larger β . Because thedisc gas far from a planet has shear flow in all cases, the fluxvalues of isothermal and non-isothermal simulations outsidethe Bondi sphere matched each other. MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Figure 1.
Midplane time sequence of shear-is-m001 run at t = , , , and . The flow pattern is represented by arrows and streamlines.The length of the arrows scales to the speed. The density is shown by colour contour. A dashed circle denotes the Bondi Radius. The emergence of an isolated inner envelope around a planetin the non-thermal simulations is due to the buoyancy bar-rier to descending gas caused by the entropy gradient. Fig-ure 6 showed the entropy structure around a planet in themodels shown in Figure 4. The entropy structure was sim-ilar for isothermal and non-isothermal simulations (exceptfor the adiabatic case): there is a low-entropy region sur-rounding the planet, which is spherically symmetric, whilethe entropy is uniform in the outer region (Figures 6a-e).The pressure structure was determined to be in hydrostaticequilibrium with the gravity of the planet. The radiativecooling made the temperature nearly uniform and equal tothat of the background shear flow. As a result, the entropywas lower around the planet. In contrast, the adiabatic simu-lation has an isentropic structure (Figure 6f). A low-entropy shell exists around the inner boundary because of an artifi-cial cooling term introduced to the innermost cell of the grid(subsection 2.2), but its influence was limited to few cells inthe vicinity of the inner boundary.The positive entropy gradient inhibits the descendinggas from reaching deeper regions in the envelope in the non-isothermal cases. Let us suppose a parcel of descending gas.Because the flow is subsonic, the pressure of the gas parcelis equal to that of the surrounding gas. In contrast, the tem-perature would increase as the gas descends because of adi-abatic compression with some radiative loss. Consequently,the parcel has a higher temperature compared to that of itssurroundings and so the density is lower in the parcel. Thedensity difference exerts buoyancy, which prevents furtherdescent. The buoyancy caused by the entropy gradient iswell-known as one of the major factors governing stellar andplanetary structures (e.g. Kippenhahn et al. 2012).
MNRAS000
MNRAS000 , 1–15 (2017)
H. Kurokawa and T. Tanigawa
Figure 2.
The same as Figure 1, but the results of shear-B2-m001 run are shown.
The isothermal and adiabatic simulations overestimatethe efficiency of atmospheric recycling. Because the temper-ature of descending disc gas is immediately equilibrated withthat of the surroundings, no buoyant force is exerted on thegas in the isothermal limit. Thus, efficient atmospheric recy-cling is allowed. Whereas in the adiabatic limit, the descend-ing gas has the same entropy as its surroundings. Again, theabsence of buoyancy allows efficient recycling.The importance of the buoyancy barrier for the emer-gence of an isolated envelope can be demonstrated by sim-plified simulations. We compared the flow topology in theisothermal and non-isothermal simulations again in Figure 7,but here we removed the complexity – a planet is embeddedin uniform headwind without a shear component and iner-tial force: the Coriolis, tidal, and global pressure forces. Thedensity structure was the same because it is determined bythe hydrostatic equilibrium with the planet gravity (Figures 7a and 7c). The temperature was nearly uniform, and, con-sequently, the entropy structure is also the same and spheri-cally symmetric (Figures 7b and 7d). The flow patterns weredifferent between the isothermal and non-isothermal runs.The isothermal run showed that the flow focuses towardsthe gravitating body because of the mass conservation: asthe density increases around the planet, the speed of the flowhas to decrease in order to conserve the mass flux, leadingto the flow converging towards the centre (Figures 7a and7b). The convergent flow in the isothermal case has beenreported previously (Lee & Stahler 2011; Ormel 2013). Onthe contrary, the headwind diverged in front of the planetand an isolated envelope formed there in the non-isothermalcase (Figures 7c and 7d). As we have omitted all the com-plexity in these runs, the emergence of the isolated envelopecan be regarded as a consequence of the buoyancy barrierdiscussed above. This interpretation will hold in more real-
MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Figure 3.
3D streamlines near the Bondi sphere in (a) shear-B2-m001 and (b) headw-B2-m001 runs at t = . Colour representswhere the radial speed is negative (blue) or positive (red). Theinner black lines show the inner domain boundary. istic simulations performing radiative transfer calculations(see subsection 4.1). Planets are exposed to headwind when disc gas rotates sub-Keplerian. The buoyancy barrier prevented the headwindfrom penetrating the planetary envelope entirely, even inthese cases where both the headwind ( M hw = . , Table1) and shear flow were taken into account. The flow pat-terns in the isothermal and non-isothermal simulations werecompared in Figure 8. The density and entropy structures were again the same because of the hydrostatic equilibriumwith the planetary gravity and the efficient radiative cool-ing. The topography of the flow differs between the isother-mal and non-isothermal simulations in the headwind case aswell as the shear-flow case (Figure 4). Whereas the isother-mal run showed that the headwind penetrated deeply intothe Bondi sphere (Figures 8a and 8b), the non-isothermalrun showed an isolated envelope around the planet wherethe atmospheric recycling was limited (Figures 8c and 8d).This is because of the buoyancy preventing the headwindfrom penetrating the envelope.The non-isothermal simulation with the headwindshares a similarity with the isothermal one when it comesto the flow outside the isolated envelope. We showed 3Dstreamlines in Figure 3b. The headwind blowing against theenvelope in the mid plane slid down along the boundary ofthe isolated envelope and left above the midplane. A similarflow pattern has been reported in the isothermal simulation(Ormel et al. 2015b), while the headwind has reached theplanetary core in that case. The buoyancy barrier prevents the atmospheric recycling re-gardless of the planetary mass. Figure 9 compares the mid-plane flow patterns of the isothermal and non-isothermalsimulations for m = . and m = . The size of the circular-streamline region is larger in the non-isothermal simulationsthan in the isothermal simulations, as with the m = . cases. The recycling and perturbation flux (Equation 9)was smaller in the non-isothermal simulations than in theisothermal simulations by an order of magnitude (Figure10). It suggests that the buoyancy barrier induced by an en-tropy gradient prevented the inflow from descending onto aplanet in the non-isothermal simulations. We note that thosehigher-mass simulations have not reached the steady stateat the time of termination: the density in the vicinity of theplanet has still been increasing gradually. Their behaviour inmore long-term simulations will be addressed in our futurework.A difference from the lower mass simulation (Figure4a) is that the streamlines in the vicinity of the planet( ∼ . R Bondi ) is circular even in the case of the isother-mal simulations (Figures 9a and 9c). Because higher-massplanets accrete more disc gas in a more 2D-fashion, theyobtain more angular momentum than lower-mass planets(more precisely, because of the vortensity conservation in2D flow, Ormel et al. 2015a). As a consequence, higher-massplanets have more circular streamlines even in the isother-mal simulations. This picture is consistent with an inferencethat planets form circumplanetary discs in the higher-masslimit (e.g., Tanigawa et al. 2012).
We showed that the buoyancy barrier prevents the disc-gas inflow from recycling the gas in the planetary envelopeby comparing isothermal and non-isothermal calculations.
MNRAS000
MNRAS000 , 1–15 (2017)
H. Kurokawa and T. Tanigawa
Figure 4.
A comparison of flow patterns in the midplane at t = : (a) shear-is-m001 , (b) shear-B4-m001 , (c) shear-B3-m001 , (d) shear-B2-m001 , (e) the large-scale structure of shear-B2-m001 , and (f) shear-ad-m001 . The flow pattern is represented by arrows andstreamlines. The length of the arrows scales to the speed, but the scale is different for each figures. The density is shown by colourcontour. A dashed circle denotes the Bondi radius. MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Figure 5.
Shell-averaged recycling flux (Equation 9). (a) Timesequence of shear-B2-m001 run and (b) a comparison betweenisothermal ( shear-is-m001 ) and non-isothermal ( shear-B4-m001 , shear-B3-m001 , and shear-B2-m001 ) runs at t = are shown. Whereas the radiative cooling was implemented by the sim-ple β cooling approximation in our model, the results areuseful to understand the difference between previous 3D sim-ulations, some of which have carried out more sophisticated,radiative-transfer calculations.Isothermal simulations of flow past a planet in a pro-toplanetary disc have pointed to complete recycling of theenvelope – no closed streamlines have been found arounda planet. Ormel et al. (2015b) simulated inviscid isother-mal flow around an embedded planet (m = 0.01) in a lo-cal frame. They found that no clear boundary demarcatesbound atmospheric gas from disc material: gas enters theBondi sphere at high latitudes and leaves through the mid-plane regions or, vice versa when the rotation of disc gas issub-Keplerian. Fung et al. (2015) presented global simula-tions of the isothermal flow around an embedded planet (m= 0.56) with viscosity. They found a “transient” horseshoeflow along which high altitude gas descends rapidly into theplanet’s Bondi sphere, performs one horseshoe turn, and ex-its the Bondi sphere radially in the midplane. The flow pre-vented the planet from sustaining a hydrostatic envelope,that is, the gas is recycling. We note that they found a boundgas in the deep part of the envelope within ∼ . times thesize of the planetary core. Because the bound region was re- solved by only 3 grids, they argued that it is uncertain howmuch gas is truly bound to the planet.Adiabatic simulations also found no dynamical bound-ary between bound atmosphere and disc gas. Popovas et al.(2018) simulated inviscid adiabatic flow around an embed-ded planet in a local frame. Their numerical set up assumedMars- to Earth-sized planets at − . AU, which corre-sponds to m = . , . and . in the dimensionlessmass used in our study. The streamlines within the Bondisphere were non-circular and connected to disc gas.In contrast, radiation-hydrodynamics simulations havereported limited recycling – an inner part of the envelopeis isolated from the recycling flow of disc gas. D’Angelo& Bodenheimer (2013) performed global 3D radiation-hydrodynamics calculations of protoplanetary discs in whicha young planetary core is located with explicit viscosity. Re-alistic equation of state and gas opacity were used and theenergy released by the accretion of planetesimals was con-sidered in their simulations. They assumed 5–15 M ⊕ plan-ets located at 5–10 AU, which corresponds to m = . –0.25. They found that the interface between the materialorbiting the planet and material orbiting the star lies at ∼ . R Bondi . The global radiation-hydrodynamics simula-tions of Lambrechts & Lega (2017) also reported that theinterface is at ∼ . times the Hill radius of the planet. Theyassumed 5 M ⊕ and larger planets at 5 AU, which corre-sponds to m > . . Cimerman et al. (2017) conducted invis-cid radiation-hydrodynamics calculations of a flow around aplanet in a local frame. The assumed planetary masses were m = . , . , . , and . . They identified an “inner core”of the envelope – ∼ . R Bondi in size – where streamlines aremore circular and entropy is much lower than in the outerenvelope.To summarise, complete recycling has been found inisothermal and adiabatic calculations; a limited recyclinghas been observed in non-isothermal calculations, regardlessof many differences in these previous studies (local/global,viscid/inviscid, equation of state, opacity, planetesimal ac-cretion, and planetary masses). The clear difference betweenprevious isothermal and non-isothermal studies is consistentwith our explanation that the high-entropy disc-gas is pre-vented from descending into the low-entropy envelope by thebuoyant force in the non-isothermal cases.We note that the entropy gradient has also been pro-posed to influence the gas accreting onto young stars: high-entropy gas accreting onto stars cannot descend into thestellar interior with low entropy (Geroux et al. 2016; Kunit-omo et al. 2017). These studies seem to be consistent withour idea that high-entropy disc gas is prevented from de-scending into the low-entropy envelope in the case of gasaccretion onto planets.
The buoyancy barrier suppresses the atmospheric recycling,which might allow the cooling and accretion of the envelopesof super-Earths. Malygin et al. (2017) estimated the radia-tive relaxation time-scale of linear temperature perturba-tions in protoplanetary discs. The relaxation time-scale t relax MNRAS000
The buoyancy barrier suppresses the atmospheric recycling,which might allow the cooling and accretion of the envelopesof super-Earths. Malygin et al. (2017) estimated the radia-tive relaxation time-scale of linear temperature perturba-tions in protoplanetary discs. The relaxation time-scale t relax MNRAS000 , 1–15 (2017) H. Kurokawa and T. Tanigawa
Figure 6.
Data are the same as Figure 4, but colour contour represents p / ρ γ , which is the measure of entropy.MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Figure 7.
A comparison of flow patterns of (a,b) heado-is-m001 and (c,d) heado-B2-m001 (headwind only). The flow pattern is representedby arrows and streamlines. The length of the arrows scales to the speed. Colour contour represents (a,c) density and (b,d) p / ρ γ (entropy),respectively. A dashed circle denotes the Bondi Radius. in the optically thick regime was estimated by, t relax ∼ l D ∼ l κ R ρχ c , (10)with l being the length scale of perturbation, D the diffu-sion coefficient, κ R the Rosseland mean opacity, c the speedof light, and χ the flux limiter (Malygin et al. 2017). A par-cel of descending gas of inflow into the Bondi sphere of aplanet (subsection 3.3) would be heated by adiabatic com-pression. The typical scale of this temperature increase canbe given by the Bondi radius: R Bondi = H × m . The esti-mated relaxation time is ∼ − Ω − for the perturbationlength scale of λ ∼ . H at 1 AU (Malygin et al. 2017, intheir nominal model without dust depletion), which corre-sponds to the size of the Bondi radius of a m = . planet.As we presented in subsection 3.2, the simulation with thisrelaxation time ( β = − ) showed the emergence of an iso- lated inner envelope. The relaxation time becomes longer forlarger scales (namely, higher-mass planets) and smaller or-bital radius (Malygin et al. 2017). Because the atmosphericrecycling becomes less efficient for larger β (subsection 3.2),our results suggest that runaway gas accretion onto super-Earths may be inevitable once the envelopes have startedcooling. Fluid close to the core is bound to the planet andmodified 1D prescriptions may be suitable, where the cool-ing controls the accretion (D’Angelo & Bodenheimer 2013;Lambrechts & Lega 2017). We note that the relaxation timedepends on the dust opacity: the dust depletion would de-crease the relaxation time (Malygin et al. 2017) and, con-sequently, the size of the isolated inner envelope, althoughdust depletion would also decrease the time-scale of envelopecooling and runaway accretion (Lee et al. 2014).There is a caveat in our analysis by using the β ap-proximation: the presence of a planet is not a linear per- MNRAS000
A comparison of flow patterns of (a,b) heado-is-m001 and (c,d) heado-B2-m001 (headwind only). The flow pattern is representedby arrows and streamlines. The length of the arrows scales to the speed. Colour contour represents (a,c) density and (b,d) p / ρ γ (entropy),respectively. A dashed circle denotes the Bondi Radius. in the optically thick regime was estimated by, t relax ∼ l D ∼ l κ R ρχ c , (10)with l being the length scale of perturbation, D the diffu-sion coefficient, κ R the Rosseland mean opacity, c the speedof light, and χ the flux limiter (Malygin et al. 2017). A par-cel of descending gas of inflow into the Bondi sphere of aplanet (subsection 3.3) would be heated by adiabatic com-pression. The typical scale of this temperature increase canbe given by the Bondi radius: R Bondi = H × m . The esti-mated relaxation time is ∼ − Ω − for the perturbationlength scale of λ ∼ . H at 1 AU (Malygin et al. 2017, intheir nominal model without dust depletion), which corre-sponds to the size of the Bondi radius of a m = . planet.As we presented in subsection 3.2, the simulation with thisrelaxation time ( β = − ) showed the emergence of an iso- lated inner envelope. The relaxation time becomes longer forlarger scales (namely, higher-mass planets) and smaller or-bital radius (Malygin et al. 2017). Because the atmosphericrecycling becomes less efficient for larger β (subsection 3.2),our results suggest that runaway gas accretion onto super-Earths may be inevitable once the envelopes have startedcooling. Fluid close to the core is bound to the planet andmodified 1D prescriptions may be suitable, where the cool-ing controls the accretion (D’Angelo & Bodenheimer 2013;Lambrechts & Lega 2017). We note that the relaxation timedepends on the dust opacity: the dust depletion would de-crease the relaxation time (Malygin et al. 2017) and, con-sequently, the size of the isolated inner envelope, althoughdust depletion would also decrease the time-scale of envelopecooling and runaway accretion (Lee et al. 2014).There is a caveat in our analysis by using the β ap-proximation: the presence of a planet is not a linear per- MNRAS000 , 1–15 (2017) H. Kurokawa and T. Tanigawa
Figure 8.
A comparison of flow patterns of (a,b) headw-is-m001 and (c,d) headw-B2-m001 (shear flow with headwind). The flow patternis represented by arrows and streamlines. The length of the arrows scales to the speed. Colour contour represents (a,c) density and (b,d) p / ρ γ (entropy), respectively. A dashed circle denotes the Bondi radius. turbation to the disc in reality. Because the gas density inEquation 10 is largely modified in the vicinity of the planet,the linear approximation is no longer valid in this region.However, our simulation showed that the descending discgas does not reach the deep envelope and stays in the outerregion where gas density is not significantly modified by theplanet. Thus, we argue that the behaviour of the recyclingflow was captured in our simulations, though a comparisonto the radiation-hydrodynamics simulations (subsection 4.1)was necessary to confirm the validity.The presence of the isolated envelope does not necessar-ily mean that the recycling mechanism cannot explain theubiquity of super-Earths. Cimerman et al. (2017) has founda low-entropy outflow leaving the planet from the inner coreof the envelope. Advection flow of gas on the low-entropycore of the envelope seems to have removed a fraction ofthe gas. They stated that studies at higher resolutions are needed to assess whether this region can become hydrody-namically isolated on long time-scales.We propose that the structures of planetary envelopesembedded in a protoplanetary disc evolved as follows (Figure11). When the accretion rate of solid materials or luminos-ity of the core was high (Phase 1), the envelope was fullyconvective (Rafikov 2006; Lambrechts & Lega 2017; Popo-vas et al. 2018) and had an isentropic profile. Complete recy-cling would have been allowed because there was no buoyantforce on the descending disc gas (Figures 4f and 6f, subsec-tion 3.3). We note that Phase 1 may not be realised fortypical values of the solid accretion rate and opacity (e.g.,Ikoma et al. 2000; Lambrechts et al. 2014; Lambrechts &Lega 2017; Lee et al. 2017). Once the envelope started to becooled (Phase 2), a radiative layer characterised by a posi-tive entropy gradient appeared. The buoyancy barrier wouldbegin to prevent the disc gas from descending into the enve- MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier Figure 9.
A comparison of flow patterns in the midplane: (a) shear-is-m01 at t = , (b) shear-B2-m01 at t = , (c) shear-is-m1 at t = , and (d) shear-B0-m1 at t = . The flow pattern is represented by arrows and streamlines. The length of the arrows scales to thespeed, but the scale is different for each figure. The density is shown by colour contour. Dashed circles denote the Bondi (black) and Hill(red) radii. lope (subsection 3.3) so that the envelope would have threelayers: convective, radiative, and recycling. The three layerstructure has been proposed by Lambrechts & Lega (2017).The radiative-convective boundary is set by the convectivestability: the local entropy gradient determines the direc-tion of buoyancy force and, consequently, the stability. Thisis the same mechanism with the buoyancy barrier which wasproposed to set the boundary for the recycling flow in thisstudy. The innermost convective layer was missing in ourmodels (Figure 4) because we neglected heat sources: theaccretion of solid materials, luminosity of the core, and therelease of energy from gas. However, the convective layer islikely to exist in reality. The local flow past a planet investigated in this study mayhelp to explain the formation of super-Earths by delayingthe growth of their solid cores.In the pebble accretion theory (Ormel & Klahr 2010;Lambrechts & Johansen 2012), proto-cores can grow up to amass heavy enough to carve a gap in the pebble disc, whichis given by (Lambrechts et al. 2014), M pebbleiso ∼ (cid:18) H / a . (cid:19) M ⊕ . (11)This mass is large enough to allow disc gas to accrete in arunaway fashion within the lifetime of a protoplanetary disc(Lee et al. 2014).Outflow in the midplane region may reduce the peb- MNRAS000
A comparison of flow patterns in the midplane: (a) shear-is-m01 at t = , (b) shear-B2-m01 at t = , (c) shear-is-m1 at t = , and (d) shear-B0-m1 at t = . The flow pattern is represented by arrows and streamlines. The length of the arrows scales to thespeed, but the scale is different for each figure. The density is shown by colour contour. Dashed circles denote the Bondi (black) and Hill(red) radii. lope (subsection 3.3) so that the envelope would have threelayers: convective, radiative, and recycling. The three layerstructure has been proposed by Lambrechts & Lega (2017).The radiative-convective boundary is set by the convectivestability: the local entropy gradient determines the direc-tion of buoyancy force and, consequently, the stability. Thisis the same mechanism with the buoyancy barrier which wasproposed to set the boundary for the recycling flow in thisstudy. The innermost convective layer was missing in ourmodels (Figure 4) because we neglected heat sources: theaccretion of solid materials, luminosity of the core, and therelease of energy from gas. However, the convective layer islikely to exist in reality. The local flow past a planet investigated in this study mayhelp to explain the formation of super-Earths by delayingthe growth of their solid cores.In the pebble accretion theory (Ormel & Klahr 2010;Lambrechts & Johansen 2012), proto-cores can grow up to amass heavy enough to carve a gap in the pebble disc, whichis given by (Lambrechts et al. 2014), M pebbleiso ∼ (cid:18) H / a . (cid:19) M ⊕ . (11)This mass is large enough to allow disc gas to accrete in arunaway fashion within the lifetime of a protoplanetary disc(Lee et al. 2014).Outflow in the midplane region may reduce the peb- MNRAS000 , 1–15 (2017) H. Kurokawa and T. Tanigawa
Figure 10.
Shell-averaged recycling flux (Equation 9). (a) A com-parison between shear-is-m01 and shear-B2-m01 at t = . (b) Acomparison between shear-is-m1 , shear-B2-m1 , and shear-B0-m1 at t = . ble/dust accretion rate onto a planetary core and, conse-quently, delay the core growth. The flow induced by a planethas been shown to reduce the accretion rate of small dustparticles as they closely follow streamlines, while pebble-sized particles are only weakly affected (Ormel 2013; Popo-vas et al. 2018).Alibert (2017) has proposed that the outflow as “ad-vection wind” combined with the vaporisation of pebbles inthe envelopes prevents the growth of the planets at massessmaller or similar to Earth mass. Our results suggested thatthe advection wind might not reach the inner envelope wherepebbles vaporise. In addition, the advection wind containingthe vaporised pebbles as heavy elements would have a higherdensity than that of surrounding gas. The buoyancy barriermight prevent the advection wind from ascending from theplanet’s gravity potential, though analysis of convective sta-bility considering both composition and temperature gradi-ents (Kurokawa & Inutsuka 2015) is needed to understandthe consequences.In either case, further studies on the interaction ofthe planet-induced wind with pebbles/dusts are needed tounderstand the consequence on the formation scenarios ofsuper-Earths. Figure 11.
A schematic view of planetary envelopes embeddedin a protoplanetary disc proposed in this study. Phase 1: Whenthe accretion rate of solid materials or luminosity of the core washigh, the envelope was fully convective. Complete recycling waspossible at the time. Case 2: When the envelope started cooling,the radiative layer appeared and the buoyancy barrier began toprevent the disc-gas flow from recycling the envelope. The inner-most envelope was still convective.
The ubiquity of super-Earths in extrasolar planetary sys-tems poses a problem for the theory of planet formation asthey are expected to become gas giants through the run-away accretion of disc gas on a short time-scale. Ormel etal. (2015b) has proposed that rapid recycling of the envelopegas of planets embedded in a protoplanetary disc would de-lay the cooling and following accretion of disc gas. However,the topography of the recycling flow has been shown to dif-fer between previous simulations of the flow past a planet ina disc.We conducted a detailed comparison between isother-mal and non-isothermal 3D hydrodynamical simulations ofthe gas flow past a planet to investigate the effect of radia-tive cooling on the recycling flow pattern. Radiative coolingwas implemented by the β cooling model.We observed the pattern of recycling flow: fluid enteredthe Bondi sphere at high latitudes and leaves through themidplane regions. The topography of inflow and outflow isthe other way around when disc gas rotates in sub-Keplerian.Whereas the isothermal simulations showed efficient recy-cling where streamlines are connected to the inner part ofthe envelope gas, the non-isothermal runs showed limited re-cycling: at the inner part of the planetary envelope, stream-lines are circular. We demonstrated that the buoyancy forceinduced by the entropy difference between the atmosphere(low entropy) and disc gas (high entropy) suppresses the re-cycling. Though our models adopted the simple β cooling,the mechanism of the buoyancy barrier is useful to under-stand the difference between the previous studies.Our results suggested that, once the atmosphere startscooling, the buoyancy barrier prevents the high-entropy disc MNRAS , 1–15 (2017) uppression of recycling by buoyancy barrier gas from intruding the cooled atmosphere, which may leadfurther cooling of the atmosphere and runaway gas accretiononto the core. Nevertheless, the interaction of the recyclingflow with accreting solid materials may delay the growth ofsuper-Earth cores and help to explain the ubiquity of super-Earths in exoplanetary systems. ACKNOWLEDGEMENTS
We would like to thank the reviewer, Michiel Lambrechtsfor thoughtful and helpful comments that have substan-tially improved the quality of this manuscript. We thankAthena++ developers: James M. Stone, Kengo Tomida,and Christopher White. HK particularly acknowledges as-sistance by Kengo Tomida to start using the code. Thisstudy has greatly benefited from fruitful discussion withChris W. Ormel, Shu-ichiro Inutsuka, Masanobu Kunit-omo, Shigeru Ida, Ayumu Kuwahara, Takayuki Saito, andMasahiro Ikoma. HK was supported by JSPS KAKENHIGrant number 15J09448 and 18K13602. TT was supportedby JSPS KAKENHI Grant number 26800229 and 15H02065.Numerical computations were in part carried out on CrayXC30 at the Center for Computational Astrophysics, Na-tional Astronomical Observatory of Japan.
REFERENCES
Alibert, Y. 2017, A&A, 606, A69Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS,471, 4662D’Angelo, G., & Bodenheimer, P. 2013, ApJ, 778, 77Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766,81Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101Gammie, C. F. 2001, ApJ, 553, 174Geroux, C., Baraffe, I., Viallet, M., et al. 2016, A&A, 588, A85Ginzburg, S., & Sari, R. 2017, MNRAS, 464, 3937Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, Protostars andPlanets II, 1100Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66Inamdar, N. K., & Schlichting, H. E. 2015, MNRAS, 448, 1751Kokubo, E., & Ida, S. 2000, Icarus, 143, 15Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structureand Evolution, Astronomy and Astrophysics Library. ISBN978-3-642-30255-8. Springer-Verlag Berlin Heidelberg, 2012Kominami, J., & Ida, S. 2002, Icarus, 157, 43Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A,599, A49Kurokawa, H., & Inutsuka, S.-i. 2015, ApJ, 815, 78Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572,A35Lambrechts, M., & Lega, E. 2017, A&A, 606, A146.Lee, A. T., & Stahler, S. W. 2011, MNRAS, 416, 3177Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95Lee, E. J., Chiang, E., & Ferguson, J. W. 2017, arXiv:1710.02604Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dulle-mond, C. P. 2017, A&A, 605, A30Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka,H. 2016, ApJ, 821, 82 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43Ormel, C. W. 2013, MNRAS, 428, 3526Ormel, C. W., Kuiper, R., & Shi, J.-M. 2015, MNRAS, 446, 1026Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512Ormel, C. W. 2017, Formation, Evolution, and Dynamics ofYoung Solar Systems, Proceedings of the Sant Cugat Forumon Astrophysics. Springer.Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus,124, 62Popovas, A., Nordlund, ˚A., Ramsey, J. P., & Ormel, C. W. 2018,arXiv:1801.07707Rafikov, R. R. 2006, ApJ, 648, 666Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon,J. B. 2008, ApJS, 178, 137-177Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47Venturini, J., Alibert, Y., Benz, W., & Ikoma, M. 2015, A&A,576, A114Weidenschilling, S. J. 1977, Ap&SS, 51, 153Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22Yu, C. 2017, arXiv:1711.00594This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000