Planets embedded in their natal discs harbour hot envelopes. When pebbles are accreted by these planets the contained volatile components may sublimate enriching the envelope and potentially changing its thermodynamical properties. However the envelopes of embedded planets actively exchange material with the disc which would limit the buildup of a vapourrich atmosphere. To properly investigate these processes we have developed a new phase change module to treat the sublimation process with hydrodynamical simulations. Combined with the recently developed multidust fluid approach we conduct 2D selfconsistent hydrodynamic simulations to study how pebble sublimation influences the water content of superEarths and subNeptunes. We find the extent and the amount of vapour that a planet is able to hold on to is determined by the relative size of the sublimation front and the atmosphere. When the sublimation front lies far inside the atmosphere vapour tends to be locked deep in the atmosphere and keeps accumulating through a positive feedback mechanism. On the other hand when the sublimation front exceeds the bound atmosphere the ice component of incoming pebbles can be fully recycled and the vapour content reaches a low steady value. Low disc temperature small planet mass and high pebble flux omitting accretion heating by pebbles render the planet atmosphere vapourrich while the reverse changes render it vapourpoor. The phase change module introduced here can in future studies also be employed to model the chemical composition of the gas in the vicinity of accreting planets and around snowlines.
['hydrodynamics', 'methods: numerical', 'planets and satellites: atmospheres', 'planets and satellites: composition', 'Astrophysics - Earth and Planetary Astrophysics']
Atmospheric recycling of volatiles by pebbleaccreting planets
{'Atmospheric Recyling of Volatiles by Pebble-Accreting Planets': 'Yu Wang ( 王 雨 ) 1 , Chris W. Ormel 1 ★ , Pinghui Huang ( 黄 平 辉 ) 2 , 3 , Rolf Kuiper 4  \n1 Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, Haidian DS 100084, Beĳing, China  \n- 2 University of Victoria, Victoria, British Columbia, Canada\n- 3 Institute for Advanced Study, Tsinghua University, 30 Shuangqing Rd, Haidian DS 100084, Beĳing, China\n- 4 Faculty of Physics, University of Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'Planets, embedded in their natal discs, harbour hot envelopes. When pebbles are accreted by these planets, the contained volatile components may sublimate, enriching the envelope and potentially changing its thermodynamical properties. However, the envelopes of embedded planets actively exchange material with the disc, which would limit the buildup of a vapour-rich atmosphere. To properly investigate these processes, we have developed a new phase change module to treat the sublimation process with hydrodynamical simultions. Combined with the recently developed multi-dust fluid approach, we conduct 2D self-consistent hydrodynamic simulations to study how pebble sublimation influences the water content of super-Earths and sub-Neptunes. We find the extent and the amount of vapour that a planet is able to hold on to is determined by the relative size of the sublimation front and the atmosphere. When the sublimation front lies far inside the atmosphere, vapour tends to be locked deep in the atmosphere and keeps accumulating through a positive feedback mechanism. On the other hand, when the sublimation front exceeds the (bound) atmosphere, the ice component of incoming pebbles can be fully recycled and the vapour content reaches a low, steady value. Low disc temperature, small planet mass and high pebble flux (omitting accretion heating by pebbles) render the planet atmosphere vapour-rich while the reverse changes render it vapour-poor. The phase change module introduced here can in future studies also be employed to model the chemical composition of the gas in the vicinity of accreting planets and around snowlines.  \nKey words: planets and satellites: atmospheres - hydrodynamics - planets and satellites: composition - methods: numerical', '1 INTRODUCTION': "Planets are born in gaseous discs. The close-in, low-mass planet population revealed by the Kepler and the Transiting Exoplanet Survey Satellite ( TESS ) missions, known as super-Earths and sub-Neptunes (Fressin et al. 2013; Winn & Fabrycky 2015; Zhu et al. 2018; Zhu & Dong2021), are believed to hold - or have held in the past - moderate amounts of H and He gas. This strongly suggests these planets have their genesis in the gas-rich disc (Fulton et al. 2017; Owen & Wu 2017; Jin & Mordasini 2018).  \nIn the core accretion model, planets need to grow massive enough before they start to bind an atmosphere (Mizuno 1980; Pollack et al. 1996). Low-mass planets can be defined to have their atmosphere in hydrostatic balance with the disc. Their imprints on discs remain limited, in contrast to high-mass planets which open gaps and that have been proposed to be responsible for the substructure seen in young discs (Benisty et al. 2017; Andrews et al. 2018; Long et al. 2018; Zhang et al. 2018; Avenhaus et al. 2018; van der Marel et al. 2019). To understand the transition from low-mass planets to gas giants, one must understand how efficient atmospheres cool (Lee et al. 2014; Ginzburg & Sari 2017; Yu 2017; Ormel et al. 2021). For the post-disc phase, atmospheres are further affected by the stellar insolation (Lopez & Fortney 2013; Owen & Wu 2017; Jin & Mor-  \ndasini 2018) or the heat release from the core (Ginzburg et al. 2018; Vazan et al. 2018a,b; Gupta & Schlichting 2019). Importantly, all these works consider a fixed composition for the gas (the standard Hand He-dominated gas) and also assume that the atmospheres stay thermodynamically isolated in the disc phase.  \nThese assumptions can be challenged, in particular when planet assembly occurs by accreting pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012). As atmospheres around low-mass planets heat up, the more volatile material component contained in the pebbles sublimates (Alibert 2017; Brouwers et al. 2018). Pebble sublimation carries major implications. First, the atmosphere will be polluted by the high-Z vapour and has a larger mean molecular weight (Iaroslavitz & Podolak 2007; Bodenheimer et al. 2018). The polluted atmosphere becomes significantly heavier than the pure H-He atmosphere, thus triggering runaway gas accretion at a lower mass (Wuchterl 1993; Venturini et al. 2015, 2016), which is renamed the 'critical metal mass' by Brouwers & Ormel (2020).  \nSecond, pebble sublimation in atmosphere can affect the volatile delivery by setting an effective iceline inside the planet atmosphere. This contrasts the standard assumption that the chemical inventory of a planet atmosphere is set by the local disc properties. That is, a planet accretes solids (ices) if it is situated outside the disc iceline, while it accretes vapour only interior to it (Öberg et al. 2011; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017; Ida et al. 2019; Booth & Ilee 2019; Krĳt et al. 2022; Mollière et al. 2022). The super-  \n𝑁 2 abundance of Jupiter, for example, has been suggested as an imprint for its core to be assembled outside the N 2 snowline (Öberg & Wordsworth 2019; Bosman et al. 2019). However, Barnett & Ciesla (2022) pointed out that an accreting planet would greatly elevate the ambient temperature and halt N 2 delivery by sublimation, putting even more stringent constraints on Jupiter's core birth locations. Similarly, Johansen et al. (2021) proposed that Earth-like planets' atmosphere can be largely refreshed by the deep recycling flows from the disc, which potentially remove the sublimated vapour and prevent Earth from becoming abundant in water and carbon, even though the respective disc snowlines move interior to the Earth's location (Johansen et al. 2023a,b,c). Knowing the efficacy of the recycling process is therefore instrumental to understand the composition of planets and to guide characterization efforts of exoplanets.  \nThe recycling mechanism has been studied with 2D and 3D hydrodynamical simulations, which solve for the flow patterns in the vicinity of small planets embedded in the disc. In 2D, the gas comes in and leaves through the horseshoe orbit and there is an inner bound region supported by rotation (Ormel et al. 2015a; Fung et al. 2015; Béthune & Rafikov 2019). On the other hand, in 3D gas tends to be accreted from the polar direction and flows back to the disc in the equatorial plane without a clear boundary between the disc and the atmosphere (Ormel et al. 2015b; Cimerman et al. 2017; Lambrechts &Lega2017;Kurokawa&Tanigawa2018;Popovasetal.2018;Fung et al. 2019; Kuwahara et al. 2019; Moldenhauer et al. 2021, 2022). Accounting for radiation transport in 3D, Moldenhauer et al. (2021) finds a fully-recycled atmosphere for an Earth-mass planet, which halts cooling and subsequent (runaway) gas accretion by continuously refreshing the low-entropy atmosphere with the higher entropy gas from the disc. Still, these simulations assume a single component H-He gas and it is unclear whether the recycling mechanism also operates in polluted atmosphere. That is, what happens to the (volatile) components that pebbles contain: do they stay in the atmosphere or do they flow back to the disc (recycle)?  \nTo investigate the fate of pebble sublimation for low-mass planets (super-Earths), we have designed a new phase change module for Athena++, on top of the recently developed multi-dust fluid module by Huang & Bai (2022). The phase change module treats the mass transfer, energy exchange and momentum conservation processes during sublimation and condensation (freezing-out) self-consistently (Li & Chen 2019), enabling us to appropriately study the coupled thermodynamic and hydrodynamical effect of phase change processes in numerical simulations for the first time. In this pioneering work, our focus lies on describing phase change processes: the sublimation of the volatile components of pebbles upon entering hot atmospheres, the thermodynamical consequences of the released vapour, and the freeze-out of vapour in the form of ice grains, which can flow back to the disc. Accordingly, we include only the minimum amount of physical processes (e.g., adiabatic equation of state, single ice species), omitting for the moment processes as feedback of pebbles on gas, grain growth, deposition of vapour on incoming pebbles, as well as radiation transport. To further reduce the simulation runtime and explore the parameter space, simulations are carried out in 2D. Our over-arching goal is to identify trends of the pebble recycling process with changing disc conditions and planet mass - trends that can be verified in the future with more detailed approaches. Simulation caveats are addressed in Sect. 4.2.  \nIn this work, we focus on a single ice species - water. We identify a dichotomy in the amount of water vapour a planet can hold on to. For recycling-dominated planets, the ice component of the incoming pebbles is fully recycled and the vapour content of the atmosphere reaches a steady, limited value. On the other hand, in the vapour-  \ndominated regime, sublimation occurs deep in the atmosphere, out of reach for recycling. A few simulations, falling in between these limits, delineate these two outcomes. In general, low ambient disc temperatures, small planets, high pebble fluxes and large Stokes numbers push the planet towards the vapour-dominated limit. In Sect. 4 we compile the simulation results to obtain a predictive estimate of which regime applies to the combination of disc condition and pebble material properties - an estimator that is imprecise in light of the above-mentioned limitations, but that we expect to convey the trends well.  \nThe structure of the paper is as follows. In Sect. 2 we describe the basic work principles of the phase change module and our implementation in different hydrodynamic codes, which we conduct a few tests to verify the robustness of our model (see also App. A). Section 3 begins with the simulation set-up of our fiducial run and then introduce the main results of this paper: various flow patterns and vapour content gained in the atmosphere with different disc parameters and planet mass. We identifiy three qualitative regimes in this section and explain the underlying mechanisms. In Sect. 4 we first summarize the results and then propose an estimator to determine the evolving trend of vapour content. Then we discuss the caveats and future improvements and applications of our phase change module. Section 5 lists our conclusions.", '2 MODEL': 'This section describes the phase change extension of the multi-fluid module. In Sect. 2.1 we outline the problem and introduce our modelling approaches. In Sect. 2.2 we present the governing equations for both gas and the multi-particle fluid hydrodynamics, building upon the work of Huang & Bai (2022), and introduce our choices of parameters and source terms regarding the problem we focus on. In Sect. 2.3 we introduce the phase change model with which we simulate evaporation fronts in hydrodynamical codes. Finally, in Sect. 2.4, we describe how we implement the model in both PLUTO and ATHENA++.', '2.1 Overall model design': 'We illustrate the problem in Fig. 1. In this figure, a planetary core (black dot at the centre) is embedded into the disc and surrounded by a denser atmosphere. Depicted in the co-moving frame, this results in the familiar horeshoe, atmosphere and circulating streamlines (grey lines; (Ormel et al. 2015a; Béthune & Rafikov 2019; Moldenhauer et al. 2021). In addition, there are pebbles, embedded in the flow, which are accreted by the planet (pebble accretion). When they approach the inner atmosphere, the temperature of the ambient gas is hot enough to sublimate the ice from the pebble, during which latent heat is absorbed from the environment. The released vapour can also flow back to the disc and recondense in the form of ice grains, which are carried away with the gas flow. On the other hand, the refractory components of the pebble continue their trajectory towards the planet surface.  \nWe simulate these processes with a multi-fluid modelling approach. In this setup, we distinguish between gas, tracer and particle fluids. The \'gas fluid\' is the mixture of all non-condensable gaseous species (mostly H and He) and all vaporized ices, which are characterized by the standard gas properties as pressure, temperature and mean molecular weight. There is only one gas fluid. In addition, the density of each vapour species is followed by a \'tracer fluid\'. The tracer fluid inherits the hydrodynamical properties (i.e., velocity  \nFigure 1. Schematic displaying the processes occurring at the sublimation fronts. Grey curves are gas streamlines as seen in the co-moving frame. Around the planetary core (black dot) a hot atmosphere forms. From top to bottom: pebbles consisting of ice and refractory materials approach the atmosphere of an embedded planet. In the hot atmosphere the ice sublimates off the pebbles, while the refractory remnant of the pebble proceeds towards the core. vapour can recondense to ice grains in colder regions in the atmosphere, which may be advected back to the disc.  \n<!-- image -->  \nfield) of the gas fluid. Next, there are the particle fluids, with which the particle species are simulated. In general, each particle species is a compound consisting of various condensates and refractory materials. For example, cm-sized \'pebbles\' are a compound whose components are refractory silicates and water ice, while micron-sized water ice grains are another single-component compound. Each of these material components is described by a pressureless \'particle fluid" which inherits the aerodynamic properties of their parent compound. In Fig. 2, we present an example of the multi-fluid setup. We assume there are two particle species A and B (varying in size e.g.) and they each contain one refractory component and one ice component, but in generally not in the same proportions. During the phase change, the ice component in the particles will be transferred to vapour, increasing the local concentration of vapour over non-condensible gas (here: hydrogen and helium). For the entire setup, we need four particle fluids to describe the refractory material and the ice of each particle compound, a single gas fluid for the H/He and vapour, and one tracer fluid to follow the vapour component corresponding to the single ice component.  \nThe simulation is separated into a hydrodynamic and a thermodynamic step. Figure 3 describes the workflow formally. Initially 𝑃 𝑖 , the pressure of vapour component 𝑖 , meets the vapour saturation pressure equation ( 𝑃 𝑖 = 𝑃 eq ,𝑖 ( 𝑇 ) , Eq. (20)) and serves as input of the simulation. The hydrodynamic step solves a set of multi-fluid hydrodynamical equations to update the density and velocities of the gas and particle. After this step, 𝑃 𝑖 = 𝑃 eq ,𝑖 ( 𝑇 ) is no longer held. For example, the saturation pressure will exceed the partial pressure of vapour due to temperature increase, which leads to sublimation. Then, the thermodynamic step is needed to update the properties of the gas (temperature, density and vapour fraction) with the phase change model, which represents the major new addition of this work.  \nFigure 2. An example of the multi-fluid setup. There are two particle species A and B containing two material components: refractory and ice, which in total need four particle fluids to follow. The gas contains one non-condensable component (H, He mixture here) and one vapour component, which are described by one gas fluid in total and one additional tracer fluid.  \n<!-- image -->  \nFigure 3. Flowchart representing the workflow of the simulation. It starts with the vapour partial pressure 𝑃 𝑖 for material 𝑖 equaling the saturation pressure 𝑃 eq ,𝑖 . And after the hydrodynamic step, the output doesn\'t satisfy it anymore. Then the thermodynamic step solves for 𝑃 𝑖 = 𝑃 eq ,𝑖 ( 𝑇 ) with mass transfer, energy exchange and momentum conservation. On the right, we list the equations related to both the hydrodynamic and thermodynamic steps.  \n<!-- image -->  \nThis model solves the mass transfer, momentum and energy exchange between ice and gas as determined by the vapour saturation pressure equation. It connects the solid and gas phases of volatiles. As a result, weare able to investigate both the hydrodynamics and thermodynamics of volatiles in a self-consistent way. We will describe both steps in detail in the next two sections.', '2.2 Governing equations': 'The hydrodynamic equations of gas and particle fluids are presented in conservative form. We use subscripts \'g" \'p" and \'tr" to denote gas, particles and the tracers respectively. We also use \'z\' to denote the material components. Furthermore, let 𝑁 p and 𝑁 z describe the total number of particle species and material components (both refractory and volatile). Thus, in total we have one gas fluid, 𝑁 z × 𝑁 p particle fluids and at most 𝑁 z tracer fluids (there are no refractory tracer fluids.).  \nThe governing equations read:  \n𝜕𝜌 g 𝜕𝑡 + ∇ · ( 𝜌 g 𝒗 g ) = 0 , (1)  \n𝜕 ( 𝜌 g 𝒗 g ) 𝜕𝑡 + ∇ · ( 𝜌 g 𝒗 g 𝒗 g + 𝑃 g 𝑰 ) = 𝜌 g 𝒇 g , src , (2)  \n𝜕𝐸 g 𝜕𝑡 + ∇ · GLYPH<2> GLYPH<0> 𝐸 g + 𝑃 g GLYPH<1> 𝒗 g GLYPH<3> = 𝜌 g 𝒇 g , src · 𝒗 g ,  \nFollowing Moldenhauer et al. (2021), we omit smoothing of the planet gravity since we aim to reveal flow patterns deep in the atmosphere. Also, we omit self-gravity. As our simulations are conducted in a local frame centred on the planet, two additional forces should be included in this non-inertial frame. The first is the Coriolis forces,  \n𝒇 cor = -2 𝛀 𝐾 × 𝒗 . (11)  \nThe second is the tidal force induced by the summation of stellar gravity and centrifugal force. Since we focus on a small domain near the planet, the tidal force can be expanded and approximated to the first order as,  \n𝒇 tid = (-𝐺𝑀 ★ / 𝑎 2 semi + Ω 2 K 𝑎 semi ) · 𝒆 𝑥 ≈ 3 Ω 2 K 𝑥 · 𝒆 𝑥 (12)  \nwhere 𝒆 𝑥 denotes the direction pointing from the star to the planet and 𝑥 is the corresponding coordinate in this comoving frame. Here, 𝑎 semi is the semi-major axis of the planet (in circular orbit). Finally we have,  \n(3) 𝒇 src = 𝒇 cor + 𝒇 tid + 𝒇 planet . (13)  \n𝜕𝜌 p ,𝑛𝑖 𝜕𝑡 + ∇ · GLYPH<16> 𝜌 p ,𝑛𝑖 𝒗 p ,𝑛 + ℱ p , dif ,𝑛𝑖 GLYPH<17> = 0 , (4)  \n𝜕𝜌 p ,𝑛𝑖 GLYPH<16> 𝒗 p ,𝑛 + 𝒗 p , dif ,𝑛𝑖 GLYPH<17> 𝜕𝑡 + ∇ · GLYPH<0> 𝜌 p ,𝑛𝑖 𝒗 p ,𝑛 𝒗 p ,𝑛 + 𝚷 dif ,𝑛𝑖 GLYPH<1> = 𝜌 p ,𝑛𝑖 𝒇 p , src ,𝑛𝑖 + 𝜌 p ,𝑛𝑖 𝒗 g -𝒗 p ,𝑛 𝑡 s ,𝑛 , (5)  \n𝜕𝜌 tr ,𝑖 𝜕𝑡 + ∇ · GLYPH<0> 𝜌 tr ,𝑖 𝒗 g + ℱ tr , dif ,𝑖 GLYPH<1> = 0 . (6)  \nIn these equations the index 𝑛 refers to the particle species and 𝑖 to the material component. The particle concentration diffusion flux is given by,  \nℱ p , dif ,𝑛𝑖 ≡ -𝜌 g 𝐷 p ,𝑛 ∇ GLYPH<18> 𝜌 p ,𝑛𝑖 𝜌 g GLYPH<19> ≡ 𝜌 p ,𝑛𝑖 𝒗 p , dif ,𝑛𝑖 (7)  \nwhich defines the effective diffusive velocity 𝒗 p,diff ,𝑛𝑖 . The particle diffusivity 𝐷 p ,𝑛 is a free parameter in our simulation related to the gas diffusivity 𝐷 g:  \n𝐷 p ,𝑛 = 𝐷 g 1 + St 2 𝑛 (8)  \nwhere St 𝑛 = 𝑡 s ,𝑛 Ω K is the dimensionless stopping time, i.e, the Stokes number (Youdin & Lithwick 2007). 𝚷 dif ,𝑛𝑖 is the momentum diffusion tensor, which describes the momentum flux caused by diffusion. Combined with 𝒗 p,diff ,𝑛𝑖 , Eq. (5) ensures the Galilean invariance (Huang & Bai 2022). And the last term of Eq. (5) denotes the aerodynamic drag felt by particles, while we omit the backreaction on gas (see Huang & Bai 2022 for a complete form of Eqs. (2) and (3) including dust feedback). Eq. (6) is the continuity equation of the vapour tracers whose velocity is the same as gas. Similar to Eq. (7), the tracer concentration diffusion flux is defined as,  \nℱ tr , dif ,𝑖 ≡ -𝜌 g 𝐷 tr ∇ GLYPH<18> 𝜌 tr ,𝑖 𝜌 g GLYPH<19> (9)  \nwhere 𝐷 tr ≡ 𝐷 g is the tracer diffusivity.  \nIn Eqs. (2) and (5) 𝒇 src denotes the source term from external forces. Planet gravity is one of the external forces:  \n𝒇 planet = -𝐺𝑀 p 𝑟 2 . (10)', '2.3 Phase change model': 'We adopt the ideal equation of state (EOS) for a multi-component fluid:  \n𝑃 = 𝑁 g ∑︁ 𝑖 𝜌 g ,𝑖 𝑘 B 𝑇 𝜇 𝑖 𝑚 p = 𝜌 g 𝑘 B 𝑇 𝜇 g 𝑚 p (14)  \nwhere 𝜌 g ,𝑖 is the density of gas component 𝑖 , 𝜇 𝑖 is the mean molecular weight of gas component 𝑖 , 𝜌 g is total density and 𝜇 g is the mean molecular weight of the gas:  \n1 𝜇 g = 𝑁 g ∑︁ 𝑖 𝑓 g ,𝑖 𝜇 𝑖 (15)  \nwhere 𝑓 g ,𝑖 = 𝜌 g ,𝑖 / 𝜌 g are the mass fraction of the gas phase components. These include the vaporized ices as well as the noncondensable materials and we denoted with 𝑁 g the total number of them. In our case, we combine H and He into one non-condensable component so that 𝑁 g = 𝑁 z + 1.  \nTo link the EOS to the specific energy of the gas, one needs to supply another caloric EOS, which is written as,  \n𝜌 g 𝑒 g = 𝑃 𝛾 -1 (16)  \nwhere for an ideal fluid 𝛾 is the heat capacity ratio and 𝑒 is the specific gas internal energy. For an ideal gas, the internal energy can be written as a summation of the specific heat of the different gas components plus the chemical potential of the vaporized ices,  \n𝜌 𝑔 𝑒 g = 𝑁 g ∑︁ 𝑖 𝐶 𝑉 g ,𝑖 𝜌 g ,𝑖 𝑇 + 𝑁 z ∑︁ 𝑖 𝜌 𝑧,𝑖 𝜙 𝑧,𝑖 (17)  \nSimilarly, for particles we have,  \n𝜌 p 𝑒 p = 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 𝐶 p ,𝑖 𝜌 p ,𝑛𝑖 𝑇 + 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 𝜌 p ,𝑛𝑖 𝜙 p ,𝑖 (18)  \nHere 𝐶 𝑉 g ,𝑖 is ideal gas heat capacity at constant volume and 𝐶 p ,𝑖 is the heat capacity of solids. In practice, only the difference in chemical potential matters. It is convenient to set the chemical potential of vapour to be zero, 𝜙 𝑧,𝑖 = 0, such that the caloric EOS stays ideal.  \nAlso, we omit the chemical potential for non-condensable components. Then,  \n1 𝛾 -1 = 𝜇 g 𝑁 g ∑︁ 𝑖 𝑓 g ,𝑖 𝐶 𝑉 g ,𝑖 𝑘 B / 𝑚 p (19)  \nwhere we combined Eqs. (14)-(17). Given the gas mass fractions 𝑓 g ,𝑖 and the heat capacity values, we hence obtain 𝛾 . Furthermore, under the ideal gas assumption, 𝛾 equals the adiabatic index, which is needed in the Riemann solver.  \nThe next step is to determine the fraction of each vapour component. During accretion, the planet atmosphere becomes much hotter than the ambient disc. For simplicity, we assume that phase change processes happen instantaneously and that all vapour components follow their saturation pressure in equilibrium. This assumption follows previous studies (Brouwers & Ormel 2020; Ormel et al. 2021), but the sublimation process can alternatively be described by rate equations (Ros & Johansen 2013; Schoonenberg & Ormel 2017). According to the Clausius-Clapeyron equation, the saturation vapour pressure of species 𝑖 reads,  \n𝑃 eq ,𝑖 = 𝑃 eq0 ,𝑖 exp (-𝑇 a ,𝑖 / 𝑇 ) (20)  \nwhere 𝑃 eq0 ,𝑖 and 𝑇 a ,𝑖 are constants specific to different materials.  \nTo consistently deal with mass and energy transfer during phase change, we use the energy and mass conservation relationship in a local grid cell following Li & Chen (2019). For a gas parcel containing various material components and particles, the total energy can be written as:  \n𝜌𝐸 = 𝜌 g 𝑒 g + 𝜌 g 𝑒 𝐾 g + 𝜌 p 𝑒 p + 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 𝜌 p ,𝑛𝑖 𝑒 𝐾 p ,𝑛 (21)  \nHere 𝑒 𝐾 = 1 2 𝒗 2 denotes the specific kinetic energy. Note that all gaseous species have the same specific kinetic energy 𝑒 𝐾 g . Moreover, latent heat is defined as the enthalpy difference between the vapour and the condensate:  \n𝐿 heat ,𝑖 ( 𝑇 ) = 𝐶 𝑃 g ,𝑖 𝑇 + 𝜙 𝑧,𝑖 - ( 𝐶 p ,𝑖 𝑇 + 𝜙 p ,𝑖 ) (22)  \nwhere 𝐶 𝑃 g ,𝑖 is heat capacity at constant pressure and 𝐿 heat ,𝑖 is the specific latent heat for vapour species 𝑖 . Therefore, the difference of the chemical potential is,  \nΔ 𝜙 𝑖 = 𝜙 𝑧,𝑖 -𝜙 p ,𝑖 = 𝐿 heat ,𝑖 ( 𝑇 ) -Δ 𝐶 𝑃,𝑖 𝑇 (23)  \nwhere Δ 𝐶 𝑃,𝑖 = 𝐶 𝑃 g ,𝑖 -𝐶 p ,𝑖 . Similarly, we can also define the specific kinetic energy difference between the solid and the gaseous phase Δ 𝑒 𝐾,𝑛 = GLYPH<0> 𝑒 𝐾 g ,𝑛 -𝑒 𝐾 p GLYPH<1> .  \nLet the mass transfer of material component 𝑖 from the solid phase to the gas phase be denoted 𝛿𝜌 𝑖 ( 𝛿𝜌 𝑖 > 0 for sublimation; also, quantities involving 𝛿 are unknown and solved for by the vapour module, while Δ -labelled quantities are constants or follow from the hydro step.). Furthermore, let the corresponding temperature change be 𝛿𝑇 . We need to partition 𝛿𝜌 𝑖 among the different particle species. For simplicity, we assume that 𝛿𝜌 𝑖 is distributed proportional to the mass fractions of the different particle species among material component 𝑖 . That is,  \n𝛿𝜌 p ,𝑛𝑖 = -𝛿𝜌 𝑖 𝜌 p ,𝑛𝑖 " 𝑁 p 𝑛 𝜌 p ,𝑛𝑖 (24)  \nThen energy conservation yields,  \n𝑁  \ng  \n∑︁  \n\'  \n›  \n«  \n\'  \nfi  \n‹  \n𝑁 z ∑︁ 𝑖 Δ 𝜙 𝑖 𝛿𝜌 𝑖 + 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 Δ 𝑒 𝐾,𝑛 𝛿𝜌 p ,𝑛𝑖 + ( 𝑇 + 𝛿𝑇 ) × (25)  \n     𝑁 g ∑︁ 𝑖 𝐶 𝑉 g ,𝑖 𝜌 g ,𝑖 + 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 𝐶 p ,𝑖 𝜌 p ,𝑛𝑖 + 𝑁 z ∑︁ 𝑖 𝛿𝜌 𝑖 GLYPH<0> 𝐶 𝑉 g ,𝑖 -𝐶 p , 𝑗 GLYPH<1>       \nIn this equation, the LHS is the total gas and particle energy before sublimation. These terms also appear on the RHS. The additional terms on the RHS represent the energy exchange due to sublimation. The first term on the RHS is caused by the chemical potential difference, which mainly accounts for latent heat absorption ( Δ 𝜙 𝑖 is usually positive). The second term on the RHS accounts for kinetic energy exchange due to the mass transfer. The last term on the RHS is due to both the temperature change and the heat capacity difference of the material components between gas and ice states.  \nEquation (25) gives a relation between 𝛿𝑇 and 𝛿𝜌 𝑖 . In addition, each 𝛿𝜌 𝑖 should fulfill the material supply limit. That is, the changes are limited to:  \n𝛿𝜌 𝑖 ≥ -𝜌 𝑧,𝑖 and 𝛿𝜌 p ,𝑛𝑖 ≥ -𝜌 p ,𝑛𝑖 . (26)  \nThe first expression stipulates that we cannot freeze out more than the amount of available vapour, the second demands that we cannot sublimate all the available ice mass.  \nWe can also include momentum conservation by assuming that sublimation happens isotropically so that the velocity of pebbles does not change. Therefore, the specific momentum of pebbles is conserved during phase change while the released vapour alters the momentum of the gas parcel during its mixing. This process can be described as,  \n( 𝜌 g + 𝑁 z ∑︁ 𝑖 𝛿𝜌 𝑖 )( 𝒗 g + 𝛿 𝒗 g ) -𝜌 g 𝒗 g + 𝑁 p ∑︁ 𝑛 𝑁 z ∑︁ 𝑖 𝛿𝜌 p ,𝑛𝑖 𝒗 p ,𝑛 = 0 (27)  \nwhere 𝛿 𝒗 g is the corresponding gas velocity change during the phase change process. This equation can be solved for 𝛿 𝒗 g and the gas specific kinetic energy 𝑒 𝐾 g is accordingly modified. Since in our simulation setup we only consider vapour freeze-out on tiny, well-coupled ice grains, where 𝒗 p ,𝑛 = 𝒗 g and 𝛿 𝒗 g = 0, momentum conservation is guaranteed. In general, however, e.g., with vapour freezing out on pebbles that drift with respect to the gas, an equation similar to Eq. (27) applicable to the condensation process, needs to be invoked.  \nIn Eq. (25), there are ( 1 + 𝑖 ) unknowns ( 𝛿𝑇 and 𝛿𝜌 𝑖 ). However, recall that we also have 𝑖 vapour saturation pressure equations for all ice components (Eq. (20)). Thus in total there are ( 1 + 𝑖 ) unknowns and ( 1 + 𝑖 ) equations. We numerically solve for the unknowns using the standard secant method (Press et al. 1986). In this procedure, the root is initialized with the values from the hydrodynamical step. By varying 𝛿𝑇 , 𝛿𝜌 𝑖 is acquired from the saturation curve and the root is found when the energy conservation is met simultaneously (Eq. (25)).', '2.4 Implementation in hydrodyanmical codes': 'Unless specified otherwise, our default setup contains one noncondensable H-He mixture, one water ice and one refractory (silicate). Their material properties are listed in Table 2. In this setup, we need one gas fluid, one tracer fluid and two dust fluid.  \n𝑁  \np  \n∑︁  \n𝑛  \n𝑁  \nz  \n∑︁  \n𝑖  \n𝑖  \n𝐶  \n𝑉  \ng  \n,𝑖  \n𝜌  \ng  \n,𝑖  \n+  \n𝐶  \np  \n,𝑖  \n𝜌  \np  \n,𝑛𝑖  \n𝑇  \n=', '2.4.1 PLUTO': 'The public version of the PLUTO code (Mignone et al. 2007, 2012) does not accommodate a dust-fluid module so that we cannot follow the hydrodynamics of pebbles directly. However, it contains a built-in chemical module. In this module the material component fractions are specified and can be updated for each timestep. For example, H-He and water vapour are two chemical components with different mean molecular weight. Thus for well-coupled ice grains (St -→ 0) which exactly follow the gas, we can also treat them as another chemical components with their mean molecular weight infinity. With this setup, PLUTO is capable to simulate special cases with well-coupled ice grains. Due to this limitation, we mainly run PLUTO to compare with the results from Athena++ and verify the design of our phase change model.', '2.4.2 ATHENA++': 'To fully capture the dynamics of pebbles as well as hydrodynamic behavior of gas simultaneously, we employ the recently developed multi-fluid dust module in Athena++ (Huang & Bai 2022). In this module, dust is modelled as a pressureless fluid. An arbitrary number of dust species, interacting with the gas via aerodynamic drag regulated by their stopping time, can be specified. Apart from the dust fluids, this module accommodates only one single gas fluid. Thus, to follow the vapour, we add several tracer fluids (see Sect. 2.1 for details). In addition, we have modified the module to allow mass transfer among the components, in accordance to the phase change model Sect. 2.3.  \nThe hydrodynamic part is calculated by the multi-fluid dust module, which solves Eqs. (1)-(5). The thermodynamic part is calculated according to the phase change model described in Sect. 2.3. Since the thermodynamic calculations also update the local sound speed, the timestep is affected by the phase change model via the CourantFriedrichs-Lewy (CFL) condition for numerical stability.', '2.4.3 Specification of thermodynamic quantities': 'In principle, 𝐶 𝑃 g ,𝑖 , 𝐶 𝑉 g ,𝑖 , 𝐶 p ,𝑖 and 𝐿 heat ,𝑖 all depend on temperature. Their precise value can be obtained from experiments combining with fittings (Fray & Schmitt 2009), where higher order temperature terms are also extended based on the classical Clausius-Clapeyron equation. Hydrogen and helium as simple molecules have their heat capacity well approximated with the ideal gas:  \n𝐶 𝑉 g ,𝑖 = 𝑛 𝑖 2 𝑘 B 𝜇 𝑖 𝑚 p (28)  \nwhere 𝑛 𝑖 is the degree of freedom of motion. On the other hand, water deviates more from this heat capacity relation (Fray & Schmitt 2009). However, in this paper we focus on demonstrating hydrodynamic effect of phase change rather than digging into precise thermodynamic properties, thus we apply simple heat capacity relationship that writes as,  \n𝐶 𝑃 g ,𝑖 = 𝑘 B 𝜇 𝑖 𝑚 p + 𝐶 𝑉 g ,𝑖 𝐶 p ,𝑖 = 3 𝑘 B 𝜇 𝑖 𝑚 p (29)  \nHere for solid heat capacity 𝐶 p ,𝑖 , we have assigned six degrees of freedom of motion (three transverse and three rotational, 𝑛 𝑖 = 6).  \nThen Eq. (19) simplifies as,  \n1 𝛾 -1 = 𝜇 g 𝑁 g ∑︁ 𝑖 𝑓 g ,𝑖 2 𝑛 𝑖 𝜇 𝑖 (30)  \nFor an ideal gas, the latent heat as a function of temperature can be expressed as:  \n𝐿 heat ,𝑖 ( 𝑇 ) = ∫ 𝑇 𝑇 𝑟 -Δ 𝐶 𝑃,𝑖 ( 𝑇 ) 𝑑𝑇 + 𝐿 heat ,𝑖 ( 𝑇 𝑟 ) ≈ 𝐿 heat ,𝑖 ( 𝑇 𝑟 ) -Δ 𝐶 𝑃,𝑖 ( 𝑇 )( 𝑇 -𝑇 𝑟 ) (31)  \nwhere 𝑇 𝑟 is the reference temperature. In the temperature range we are interested in, the experimental results and analytic solutions show that 𝐿 heat ,𝑖 ( 𝑇 ) varies little but has large systematic uncertainties (Feistel &Wagner 2007). Considering this and our focusing on hydrodynamics, we decide to use a constant latent heat value in simulation. Thus Eq. (23) reduces to:  \nΔ 𝜙 𝑖 = 𝜙 p ,𝑖 -𝜙 𝑧,𝑖 = -𝐿 heat ,𝑖 ( 𝑇 𝑟 ) (32)  \nHere after we omit 𝑇 𝑟 and list the value used in our simulation with 𝐿 heat in Table 2.', '2.5 disc model': "Throughout the simulation runs, we adopt the Minimum Mass Solar Nebula (MMSN) disc model (Weidenschilling 1977; Hayashi 1981). The surface density and temperature profile read:  \nΣ ( 𝑎 semi ) = 1 . 7 × 10 3 GLYPH<16> 𝑎 semi 1 au GLYPH<17> -3 / 2 g cm -2 (33)  \n𝑇 ( 𝑎 semi ) = 270 GLYPH<16> 𝑎 semi 1 au GLYPH<17> -1 / 2 K (34)  \nThe mass of the planetary core is characterized by a dimensionless parameter 𝑚 = 𝑀 p / 𝑀 th , where 𝑀 th = 𝑐 3 s , iso /( 𝐺 Ω K ) is the thermal mass. Here 𝑐 s , iso is the isothermal sound speed corresponding to the disc temperature. In the MMSN disc, the dimensionless mass 𝑚 and the planetary core mass are related as follows,  \n𝑀 p 𝑀 ⊕ ≈ 11 . 8 𝑚 GLYPH<16> 𝑎 semi 1 au GLYPH<17> 3 / 4 . (35)  \nAlso, the Bondi radius of the planet is defined as,  \n𝑅 Bondi = 𝐺𝑀 p 𝑐 2 s , iso ≡ 𝑚𝐻 g (36)  \nwhere 𝐻 g = 𝑐 s , iso / Ω is the local scale height at the planet's orbital radius. Unless further specified, we assume an inviscid gas flow following previous studies (Ormel et al. 2015a; Béthune & Rafikov 2019; Kurokawa & Tanigawa 2018). However, in order to speedup the simulation convergence we do adopt a small but non-zero diffusivity for the vapour and ice grains.  \n𝐷 g = 𝛽𝑐 s , iso 𝐻 g (37)  \nHere we take 𝛽 = 10 -5 for our simulation (see Table 1).  \nUnless specified, we use a set of dimensionless units to show the results, where the length is in unit of local scale height 𝐻 g, time is in unit of Ω -1 k and density is in units of the disc midplane density at planet position: 𝜌 g , 0 = Σ /( √ 2 𝜋𝐻 g ) . The temperature is still expressed in Kelvin.  \nTable 1. Parameters of all simulations runs. Columns denote 𝑎 semi :semi-major axis of planets; 𝑚 : dimensionless planet mass (see Sect. 2.5); 𝑁 p: number of particle species; material: material components included in pebbles; 𝑓 𝑧 : initial material component fractions for the particle species with non-zero St; St: dimensionless stopping time of each particle species with respect to the background disc; 𝑓 p2g , ini : initial pebble-to-gas ratio; 𝐷 g: diffusivity of gas and vapour in terms of 𝐻 2 Ω -1 ; 𝑡 soft : planet gravity softening time; RES: resolution of simulation ( 𝑁 𝑟 × 𝑁 𝜙 ); domain: simulation domain ( 𝑟 min , 𝑟 max ) × ( 𝜙 min , 𝜙 max ) in unit of local gas scale height ( 𝐻 g) and radian.", '2.6 Pebble dynamics': 'We include two drag regimes: the Epstein regime and the Stokes regime drag. These drag regimes are determined by the size of the pebble in relation to the mean free path of the gas molecules 𝑙 mfp ,  \n𝑙 mfp = 𝜇 g 𝑚 p √ 2 𝜌 g 𝜎 mol . (38)  \nHere, 𝜎 mol is the molecular collision cross-section, which we approximate by the collision cross-section of molecular hydrogen 𝜎 mol = 2 × 10 -15 cm 2 (Chapman & Cowling 1991). The stopping time is given by (we drop the subscript 𝑛 ),  \n𝑡 s =              𝜌 · , p 𝑠 p 𝑣 th 𝜌 g , ( Epstein : 𝑠 p < 9 4 𝑙 mfp ) , 4 𝜌 · , p 𝑠 2 p 9 𝑣 th 𝜌 g 𝑙 mfp , ( Stokes : 𝑠 p > 9 4 𝑙 mfp ) . (39)  \nwhere 𝑠 p is the particle size, 𝜌 · , p is the particle internal density and 𝑣 th = √︁ 8 / 𝜋𝑐 s is the thermal velocity of the gas molecules. Here 𝑐 s is the local sound speed. For a particle compound, its internal density 𝜌 · , p is given by,  \n𝜌 · , p = \' › 𝑁 𝑧 ∑︁ 𝑖 𝑓 p ,𝑖 𝜌 · , p ,𝑖 \' fi -1 (40)  \n«  \n‹  \nwhere 𝑓 p ,𝑖 = 𝜌 p ,𝑖 /( " 𝑁 𝑧 𝑖 𝜌 p ,𝑖 ) is the fraction of material component 𝑖 and 𝜌 · , p ,𝑖 is the internal density of material component 𝑖 (see Table 2 for detailed value).  \nIn order to calculate stopping time, we need to compute the particle size 𝑠 p. Since the mass of the refractory component of each particle species is unaffected by the phase change process, we obtain the particle number density 𝑛 p accordingly,  \n𝑛 p = 𝜌 p , refrac 𝑚 p , refrac , 0 (41)  \nwhere 𝑚 p , refrac , 0 is the mass of refractories in this particle species and 𝜌 p , refrac the corresponding internal density. From the number  \ndensity, the particle mass is obtained, 𝑚 p = " 𝑁𝑧 𝑖 𝜌 p ,𝑖 / 𝑛 p, which, combined with its internal density, gives 𝑠 p:  \n4 𝜋 3 𝑠 3 p = 1 𝑛 p " 𝑁 𝑧 𝑖 𝜌 p ,𝑖 𝜌 · , p (42)', '2.7 Model validation': 'To verify the hydrodyanmical part of our model, we reproduce the case that inviscid gas flow composed of pure H and He passes through a planetary atmosphere, a problem that has already been studied before (Ormel et al. 2015a; Fung et al. 2015; Béthune & Rafikov 2019). By design, the gas is adiabatic and 𝛾 = 1 . 43 (assume 71% H 2 and 29% He). We also run cases with 𝛾 = 1 . 3 , 1 . 1. Overall the flow patterns show the typical horseshoe orbit and inner rotational atmosphere that are consistent with previous studies (Fig. A1). We also reproduce the trend that the atmosphere expands with 𝛾 decreasing. However, there emerges unphysical entropy violation possibly caused by numerical viscosity for the run 𝛾 = 1 . 1 (Fig. A2). Please refer to Sect. A for detailed results and discussions.  \nTo test the robustness of the phase change module, we apply it to the 1D snowline model where Schoonenberg & Ormel (2017) investigated the feasibility of planetesimal formation by streaming instability near the water snowline. They found that the solid surface density would be enhanced through vapour diffusing outwards across the snowline to recondense back to ices. We have also successfully reproduced the ice density enhancement just outside the snowline seen in their fiducial model.', '3 RESULTS': "In this section, we conduct simulations with either well-coupled grains or pebbles of non-zero Stokes number. All runs are listed in Table 1. For non-zero Stokes number runs, we name them with the prefix ' pebble ' and the text following the dash represents the exclusive parameters that are varied in this run. For example, pebble-fiducial represents the fiducial run where the semi-major axis of the planet  \nTable 2. Thermodynamic constants of different material components. The given ID refers to the 'material' column in Table 1.  \n1  \nFray & Schmitt (2009)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 4. Comparison between well-coupled cases for the Athena++ (upper panels) and PLUTO (lower) runs ( coupled and coupled-Pluto ). The color in three columns denotes ice grain fraction ( 𝑓 grain , ice ), vapour fraction ( 𝑓 vap) and gas density in logarithmic scale, respectively. The text in the title indicates the simulation name (Table 1). The critical streamline for gas (blue) is plotted in the middle and right columns (Sect. 3.2).  \n<!-- image -->  \n( 𝑎 semi ) is 4.0 au. Then pebble-3.5au means the planet is put at 3.5 au while other parameters are kept the same as the fiducial run. The number of particle species included in this run is given by 𝑁 p. For each particle species we list their Stokes numbers (St) and initial pebble-to-gas mass ratio ( 𝑓 p2g , ini ). The Material column represents the material components contained in each particle species and their properties and corresponding ID are listed in Table 2. Specifically, (b) represents water and ETC. For each material component, the initial mass fraction is given by 𝑓 𝑧 . For example, for the fiducial run, only a single particle species with 𝑓 p2g , ini = 0 . 02 is initialized in the domain and water ice accounts for 50% of it.", '3.1 Simulation domain and setup': "To better capture the hydrodynamics near the core, a polar coordinate with a logarithmically spaced grid in the radial direction is adopted (Ormel et al. 2015a), while in the azimuthal direction the grid is evenly spaced. With this design the outer boundary can be placed far enough to not feel the perturbations by the planet. Therefore, we  \nuse a fixed outer boundary with constant density and a Keplerian shear velocity (both for pebbles and gas). As a consequence, gas and pebbles are replenished continuously from the outer boundary. For the inner boundary, if not specified otherwise, we use reflective boundary condition for gas and outflow for pebbles. We list our simulation resolution for all simulation runs in Table 1.  \nIn order to have a smooth start without shocks near the planet surface, a gravitational softening coefficient 𝑓 soft multiplying the planet's gravity is introduced. Here, we choose  \n𝑓 soft = 1 -exp [-0 . 5 ( 𝑡 / 𝑡 soft ) 2 ] (43)  \nfollowing Ormel et al. (2015a); Moldenhauer et al. (2021), where 𝑡 soft is the gravitational softening time.  \nIn the St = 0 simulations, we only consider one particle species: the well-coupled grains which contain only one material component, water ice. Thus in total the simulation contains one dust fluid. In addition, one tracer fluid is needed to follow the vapour density. Initially, well-coupled grains (St = 0) are present in the domain at a pebble-to-gas ratio 𝑓 p2g , ini (Table 1). Also grains with same pebbleto-gas ratio 𝑓 p2g , ini continuously replenish from the outer boundary,  \n2  \nCosta et al. (2017)  \n3  \nGoldfarb & Suuberg (2008)  \nTable 3. Summary of several characteristic outputs of all runs. Columns denote: 𝑅 subl : radius of the sublimation front; 𝑅 atm: atmosphere radius; 𝑓 vap , peak : peak vapour fraction; 𝑀 vap: total vapour mass in the whole domain; 𝑀 vap , bound : bound vapour mass inside 𝑅 atm. In status, we summary the runs's behavior by three regimes defined in the main text (recycling-dominated as 'R', vapour-dominated as 'V' and the intermediate regime as 'I'). For the runs that have already reached a steady state, we list their steady value of all output quantities (take the average of last 10 Ω -1 ) while for the other runs (either quasi-steady or not steady), we list their value for the final simulation time.  \nFigure 5. Comparison of azimuthally averaged quantities of runs gasOnly and coupled . From top to bottom, figures show profiles of midplane gas density, vapour fraction and temperature.  \n<!-- image -->  \nrendering a steady state solution possible. Therefore, it is expected that the vapour fraction in the inner region (near the core) should level off at ( 𝑓 p2g , ini )/( 1 . 0 + 𝑓 p2g , ini ) .  \nIn the St > 0 simulations, we consider two material components: (volatile) water ice and (refractory) silicates, along with the noncondensible H/He gas mixture (assume 71% H 2 and 29% He). For simplicity, we assume that these materials are contained into pebbles of a single size (or Stokes number), of which only the H 2 O ice can  \nsublimate. We further assume that the water vapour will recondense onto grains that are essentially massless, St = 0. Hence, there are two particle species - pebbles and ice grains - represented by in total three dust fluids. Also, one tracer fluid is added to follow the vapour density. For simplicity, we omit condensation of water vapour on pebbles and other processes as grain growth and fragmentation.  \nIn contrast to the gas-only and perfectly-coupled simulations, we found that the St > 0simulations must be carefully initialized to avoid any spurious instabilities or vortices triggered by the presence of decoupled pebbles near the planetary core. Thus an 'relaxing' boundary condition was introduced, which starts from an inflow/outflow boundary but gradually relaxes to reflective. Specifically, the velocity at the ghost cell is written as,  \n𝑣 ' g ,𝑟 = -2 𝑣 g ,𝑟 𝑓 soft + 𝑣 g ,𝑟 (44)  \nwhere 𝑣 ' g ,𝑟 is the radial gas velocity of the ghost zone and 𝑓 soft = 1 -exp [-0 . 5 ( 𝑡 / 𝑡 soft ) 2 ] , which is the same as used for the gravitational softening factor (Eq. (43)). In our runs, we choose 𝑡 soft = 5 Ω -1 K and larger 𝑡 soft only bring slight difference on the steady state vapour fraction, which will not influence our conclusions.  \nIn App. A we report an entropy violation issue, possibly induced by the high azimuthal velocity . For this reason, we increase the size of the inner radius with respect to the gasOnly and coupled runs and lower the resolution accordingly in the St ≠ 0 simulations (Table 1). Moreover, we conduct only half-domain simulations where only 𝜙 = [ 0 , 𝜋 ] is simulated considering that all the forces and initial conditions are symmetric about the central planet. In runs that nevertheless yielded unacceptably high numerical viscosity, we have increased the resolution (for example, pebble-m0.2 ).", '3.2 Strongly coupled dust runs': 'Since well-coupled grains follow the gas, we can implement them both in PLUTO with Chemical module (Sect. 2.4.1) and Athena++ with the multifluid dust module. By comparing their outputs, we are able to verify the robustness of the phase change module presented in Sect. 2.3.  \nIn Fig. 4, we show the results from both PLUTO and Athena++ after a steady state has emerged for the run with 𝑓 p2g , ini = 0 . 25 (Table 1, coupled and coupled-Pluto ). The flow pattern for the PLUTO and Athena++ runs are nearly identical. They both show symmetric horseshoe orbits and central rotational atmospheres. Near  \n<!-- image -->', 'pebble-fiducial': "<!-- image -->  \nFigure 6. Flow patterns of run pebble-fiducial . From left to right, panels depict ice density contained in pebbles, the H2O vapour fraction, and the ice density contained in grains. The black dashed line in each panel denotes the Bondi radius and the orange circle denotes the sublimation front (Sect. 3.3). The critical streamline (Sect. 3.2) for gas (blue) is plotted in the middle panel (vapour) while the 'ice sublimation band' (Sect. 3.3) where pebbles sublimate is denoted in the left panel by the red streamlines. A movie showing the temporal evolution is available in the online supplementary material and on github (pebble-fiducial).  \n<!-- image -->  \nthe core where the environment is hot, ice grains (left) sublimate and convert to vapour (centre). Further from the core, no vapour exists. There is a sharp increase in the vapour fraction as well as a sharp decrease in the ice grain fraction, which are due to the exponential form of the saturation profile (Eq. (20)). Following the horseshoe orbit, some vapour outside the blue streamline can be transported to cooler region and freeze out immediately as ice grains, which flow back to the disc (right panel). In contrast, the vapour inside the blue streamline rotates around the planet on closed streamlines in a bound atmosphere. We then define the blue streamline as the critical streamline that distinguish the recycling efficacy of H 2 O. Most incoming pebbles already sublimate before they enter the critical streamline and do not contribute to the vapour content in the planetary atmosphere.  \nTo further investigate the result, in Fig. 5 we present azimuthallyaveraged profiles of gas density, vapour fraction and temperature for runs gasOnly and coupled . The middle panel shows the radial profile of the vapour fraction. As expected, it increases sharply from the outer region and reaches the value of 0.2 for the inner atmosphere. The sharp increase denotes the sublimation front, where the ice grains sublimate and the vapour freezes out. Due to the injection of vapour the density of the coupled run is higher than the gasOnly run. There emerges a temperature plateau at the sublimation front where the temperature profile becomes flatter and lower than that in the gasOnly run. This plateau is due to the latent absorption associated with the phase change process, which suppress the temperature increase. Interior of the plateau where all ice is sublimated, the temperature begins to increase again at a steeper rate than the gasOnly run. This is because of a lower 𝛾 value caused by vapour injection (see Eq. (30), Brouwers & Ormel (2020)).  \nThese simulations demonstrate that the phase change module work for different hydrodynamic code and that it models thermodynamic processes like latent heat energy exchange as expected.", '3.3 Fiducial run': 'For the fiducial run ( pebble-fiducial in Table 1), the planet position is fixed at 4 au and the dimensionless mass 𝑚 = 0 . 1 ( ∼ 3 . 3 𝑀 ⊕ at 4 au). The Stokes number of the pebbles is 0.01 and the pebble-to-gas ratio is 0.01.  \nThe flow patterns corresponding to this run at time 𝑡 = 120 Ω -1 are shown in Fig. 6, when a steady state pattern has emerged. The left panel shows that pebbles enter the atmosphere where they sublimate to vapour. The streamlines clearly show the inspiralling trajectories of the accreted pebbles. To indicate how much ice is sublimated, we denote the \'ice sublimation\' band with red lines. Inside the red line, pebbles get accreted inside the sublimation front and the ice contained sublimates to vapour while outside it pebbles leave directly via the horseshoe. We note that our definition is different from the typical pebble accretion band, which denotes the accretion of refractory pebbles to the core. Because of adiabatic compression, the temperature of the ambient gas starts to increase around the planet Bondi radius. For the parameters of this simulation, the sublimation front 𝑅 subl - the radius where the ice has fully sublimated from the pebbles - happens to be around this point. Within 𝑅 subl all H 2 O is in vapour form (middle panel).  \nSpecifically, the radius of the sublimation front 𝑅 subl is defined as  \n𝜌 ice ( 𝑟 < 𝑅 subl ) = 0 (45)  \nwhere 𝜌 ice = 𝜌 peb , ice + 𝜌 grain , ice is the total ice density and 𝜌 ice is its azimuthal average. In addition, we define an atmosphere radius 𝑅 atm as the minimum distance from the critical streamline to the planetary core 1 . Approximately, vapour inside 𝑅 atm is harder to be recycled back to the disc than vapour outside it. Clearly, 𝑅 subl > 𝑅 atm in the fiducial run, which implies that the vapour can readily flow back to the disc. We list the value of 𝑅 subl and 𝑅 atm in Table 3.  \nTo examine whether the atmosphere vapour content has reached a steady value, we measure the mass influx of ice contained in pebbles, F peb , ice , and the outflux of recycled ice grains, F grain , ice . In Fig. 7 both fluxes are measured at 3 𝑅 Bondi and the sum of them gives the net ice flux. In addition, Fig. 7 presents the total vapour mass 𝑀 vap inside the whole domain (the red thick line). Not all the vapour is bounded to the planet; the vapour outside the critical streamline is on its way to freeze out as ice grains. For reference, we also give the amount of vapour inside 𝑅 atm ("star" symbol), which we define as the bound vapour mass 𝑀 vap , bound (see Fig. 7). Initially, the net flux  \nFigure 8. Azimuthally-averaged radial profiles of run pebble-fiducial . The 𝑦 -axis denotes the material fraction normalized by the gas density (e.g., 𝑓 peb , ice = 𝜌 peb , ice / 𝜌 g). The black dash-dotted line represents the saturation profile and 𝑓 ice = 𝑓 peb , ice + 𝑓 grain , ice denotes the total ice fraction. The silicate fraction has a small bump near the core, which is caused by the Keplerian rotation of the gas that slows down pebble infall.  \n<!-- image -->  \nFigure 7. Mass flux and total vapour mass of run pebble-fiducial . The 𝑥 -axis denotes the simulation time. The left 𝑦 -axis represents indicates the pebble, grain, and net mass flux and the right y-axis (red) gives the total vapour mass inside the atmosphere. The final total vapour mass 𝑀 vap is marked with a red circle while the bounded vapour mass 𝑀 vap , bound (the vapour mass inside 𝑅 atm) is marked with a star.  \n<!-- image -->  \nis negative, indicating net vapour injection. Quickly, after ∼ 4 𝑡 soft , the net flux converges to zero and the total vapour mass reaches a steady value. The incoming pebbles get completely recycled back to the disc after sublimation and then freeze-out to grains.  \nThe steady state result is better illustrated in Fig. 8, where the material fractions are shown (e.g., 𝑓 peb , ice = 𝜌 peb , ice / 𝜌 g). 𝑓 satur represents the saturation profile calculated from the /local temperature and Eq. (20). The position of the sublimation radius as defined above is indicated by the vertical grey line. The vapour fraction 𝑓 vap reaches a steady value of ≈ 0 . 013 inside the sublimation front, which we defineas 𝑓 vap , peak . Here 𝑓 vap , peak > 𝑓 p2g , ini /( 1 . 0 + 𝑓 p2g , ini ) , higher than the coupled run. This increase is a consequence of the gravitational focusing experienced by the pebbles.  \nIn summary, since 𝑅 subl > 𝑅 atm, incoming pebbles already sublimate before they enter the atmosphere. A steady-state emerges where vapour injection via incoming pebble sublimation is fully balanced by freeze-out of ice grains that are recycled back to the disc. The peak vapour fraction in the atmosphere 𝑓 vap , peak is higher than what is expected for well-coupled grains but remains low.', '3.4 Varying the planet position': "In this section, we vary the semi-major axis of the planets while keeping other parameters the same as in the fiducial run (see pebble-3.5au , pebble-5au and pebble-6au in Table 1). Figure 9 is a collection of flow patterns for several parameter variations with the centre panel that of the fiducial run. In these panels we combine the vapour and ice grains fractions, which share the same velocity field. Inside the sublimation front the intensity of the color denotes the vapour fraction while outside it the color denotes the ice grain fraction. Along each of the four colored lines, one parameters is varied with respect to the fiducial run. In this section, we are concerned with the middle panels. In addition, Fig. 10 presents several output quantities as function of time to quantitatively investigate the evolutionary trend.  \nThe 3.5 au run is qualitatively similar to the fiducial run. The time evolution of the characteristic outputs, shown in Fig. 10, reach a steady state quickly, which is similar to the default run (4 au). For both runs, 𝑅 subl > 𝑅 atm and the ice contained in incoming pebbles  \nis fully recycled by freezing-out ice grains. At 3.5 au, however, the sublimation front is significantly larger than the Bondi radius, which is a result of the increasing ambient temperature at this planet location (Eq. (34)). Also, further from the planet, pebbles experience less significant gravitational focusing by the planet and remain more tightly-coupled to the gas. This results in the drop of the peak vapour fraction for the 3.5 au run ( ≈ 0 . 01), which is only slightly above the value of the perfectly-coupled runs ( ≈ 0 . 0091). Although the 3.5 au run contains a higher total vapour mass (dashed line) than the 4 au run, less vapour is bound to the atmosphere (star, Sect. 3.3).  \nIn contrast, the run pebble-5au does not reach a steady state within the simulation time ( ∼ 900 Ω -1 ) and the snapshot at 𝑡 = 300 Ω -1 is shown in Fig. 9. Compared to the fiducial run, 𝑅 atm expands while 𝑅 subl shrinks significantly. Consequently, the sublimation front lies deep inside the atmosphere. Ice keeps being accreted (purple lines) while the vapour cannot (easily) escape the atmosphere, leading to a sustained growth of the vapour content (see Fig. 10). Although vapour does freeze-out, the ice grains similarly remain locked in the atmosphere and keep accumulating. The pebble-6au run (red lines) results in a picture similar to the 5 au run with an even more rapid vapour growth and higher vapour fraction. A positive feedback loop is at work to account for these features. Initially, the sublimation fronts of the 5 and 6 au runs are already smaller than their atmosphere radii due to the lower ambient temperature (see Fig. 10). Cooling caused by ice sublimation further shrinks the sublimation front. At the same time the injection of vapour grows the atmosphere in terms of its mass, which also causes it to expands (Ormel et al. 2015a; Moldenhauer et al. 2022). The divergence of these two critical radii suppresses the efficacy of the ice recycling. Hence, we see a continuous growth of the peak vapour fraction and vapour mass. Because of the significant expansion of the atmosphere, nearly all vapour is bound, see Fig. 10.  \nTo verify this positive feedback, we plot the azimuthally-averaged radial profiles of pebble-5au in Fig. 11. The color of the lines indicates time in units of Ω -1 . First, the temperature profiles (f) keep evolving, extending the plateau (Sect. 3.2) towards the core because of the continuous cooling caused by ice sublimation. This implies significant shrinking of 𝑅 subl . Second, the velocity profile  \nFigure 9. Flow patterns of 'pebble' runs (listed in Table 1) following the simulation setup of Sect. 3.1. The fiducial run is shown in the centre while other panels show runs that have different parameters varied. Again the yellow circle denotes the sublimation front, the black dashed circle denotes the Bondi radius, the blue line denotes the critical streamline and the red line indicates the ice sublimation band. In each plot, the magenta-to-brown color inside the sublimation front denotes vapour fraction ( 𝑓 vap) while the intensity of the blue color outside it denotes ice grain fraction ( 𝑓 ice , grain ). Movies corresponding to pebble-3.5au and pebble-5au are available as online supplementary material and on github (pebble-3.5au and pebble-5au).  \n<!-- image -->  \nFigure 10. Characteristic outputs of runs of different location evolving with the simulation time. Left panel shows the atmosphere radius 𝑅 atm (transparent thick line) and the sublimation front 𝑅 subl (opaque thin line) together in unit of disc scale height. The Bondi radius 𝑅 Bondi and inner boundary size 𝑅 min are also indicated with black dashed lines; On the right, the upper panels shows the peak vapour fraction 𝑓 vap , peak while the lower panel vapour mass 𝑀 vap. The left y-axis (black) represents the value of pebble-3.5au and pebble-fiducial while the right y-axis (red) represents pebble-5au and pebble-6au to avoid showing data that differs too much with the same axis. The dots in the last panel denote the total vapour mass contained in the domain at the final time while the stars denote the bound vapour mass.  \n<!-- image -->  \nFigure 11. Time evolution of azimuthally-averaged radial profiles of run pebble-5au . From left to right, top to bottom, the 𝑦 -axes show the norm of the velocity vector (normalized by the planet's Keplerian velocity), gas pressure, gas density, vapour fraction, total ice density and temperature, respectively. The 𝑥 -axes denote distance from the planet's centre. Line color indicates simulation times, as given by the colorbar on the right. The black dashed line in the first panel represents the Keplerian limit, 𝑣 / 𝑣 k = 1 . 0.  \n<!-- image -->  \n(a) also evolves. The radius where the azimuthally-averaged velocity displays its minimum, indicating the transition from a Keplerian shear-dominated velocity to a rotationally-dominated velocity in the planetary atmosphere, shifts outwards in time (from 𝑟 ∼ 0 . 1 to 0 . 3 in Fig. 11 a), indicating that 𝑅 atm expands. The expansion of the atmosphere is exacerbated by the decrease of the adiabatic index 𝛾 due to the increasing vapour fraction (Eq. (30)), which leads to more material falling in, as can also be seen in the Gas-only runs (Sect. A). Finally, the gas density (c) and pressure (b) also increase  \nsignificantly. Near the core, the density increase by a factor of ∼ 15 while the vapour fraction (d) only become ∼ 4 times larger, meaning that the entire atmosphere becomes significantly heavier and the additional accretion of H and He contributes more than the vapour.  \nIn summary, the respective values of 𝑅 atm and 𝑅 subl are the key parameters that determine the evolutionary trend of vapour content. When the sublimation front exceeds the (bound) atmosphere, incoming pebbles can be fully recycled and the vapour content reaches a steady value (recycling-dominated). On the other hand, when the  \nFigure 12. Characteristic outputs of runs of different planet masses evolving with the simulation time. Left and right panels are similar to Fig. 10 except 𝑅 atm and 𝑅 subl are in unit of Bondi radius. The vertical dashed line in each panel denotes 𝑡 = 300 Ω -1 , which is the starting point of the run with planet mass growth from 𝑚 = 0 . 05 to 0.2.  \n<!-- image -->  \nsublimation front is inside the atmosphere, vapour tends to be locked deep in the atmosphere and keeps accumulating. The positive feedback works effectively to shrink the sublimation front and expands the atmosphere, leading to continuous growth of the vapour content (vapour-dominated) and the entire atmosphere. Consequently, the extent and the amount of vapour that a planet is able to hold on to strongly depend on its location in the disc (ambient temperature).", '3.5 Varying the planet mass': "In this section, we vary the planet mass while keeping other parameters the same as in the fiducial case (see pebble-m0.05 , pebble-m0.2 ) in Table 1. The corresponding flow patterns are shown in Fig. 9 along the upper diagonal (red line). The run pebble-m0.2 has already reached a steady state while the pattern of the pebble-m0.05 run is still evolving and we show the snapshot at 𝑡 = 500 Ω -1 . In Fig. 12 the characteristic outputs are plotted where 𝑅 atm and 𝑅 subl are in units of their respective Bondi radius.  \nThe vapour fraction of the high-mass run pebble-m0.2 quickly reaches a steady value that is slightly lower than that of the fiducial case. Initially, just after the gas infall time (at 𝑡 ∼ 4 𝑡 soft = 20 Ω -1 ), its sublimation front lies already significantly further away from the atmosphere radius 𝑅 atm, leading to the recycling-dominated regime. We now explain the mass-dependence of 𝑅 subl / 𝑅 atm.  \nFirst, 𝑅 subl can be obtained by assuming that at the sublimation front the saturatd vapour fraction reaches 𝑓 vap , peak :  \n𝑓 vap , peak = 𝑃 eq , 0 exp GLYPH<16> -𝑇 a ,𝑖 𝑇 subl GLYPH<17> 𝑃 subl 𝜇 H2O 𝜇 g (46)  \nwhere 𝑃 subl and 𝑇 subl are the gas pressure and temperature at the sublimation front respectively, and 𝜇 H2O is the mean molecule weight of water. Given that the adiabatic index 𝛾 changes little for the recyclingdominated runs, the temperature and pressure profile follows the 1D adiabatic solution:  \n𝑇 subl = GLYPH<20> 1 + 𝛾 -1 𝛾 𝑅 Bondi GLYPH<18> 1 𝑅 subl -1 𝑟 max GLYPH<19> GLYPH<21> 𝑇 disc 𝑃 subl = GLYPH<20> 1 + 𝛾 -1 𝛾 𝑅 Bondi GLYPH<18> 1 𝑅 subl -1 𝑟 max GLYPH<19> GLYPH<21> 𝛾 𝛾 -1 𝑃 disc (47)  \nwhere 𝑇 disc and 𝑃 disc are the disc background temperature and pressure and 𝑟 max denotes outer boundary of the simulation domain. From Eqs. (46) and (47), 𝑅 subl / 𝑅 Bondi is solely determined by 𝑓 vap , peak given disc parameters. Initially, within the gas infall time ( ∼ 4 𝑡 soft ), all runs shown in Fig. 12 have a vapour fraction similar to their initial pebble-to-gas ratio ( ∼ 0 . 01) since H 2 O vapour has not yet accumulated significantly. Hence, they share a similar 𝑅 subl / 𝑅 Bondi ∼ 1 . 2.  \nOn the other hand, 𝑅 atm / 𝑅 Bondi shows a clear trend with mass. Initially (at 𝑡 ∼ 4 𝑡 soft ), 𝑅 atm / 𝑅 Bondi is smaller with increasing planet mass (Fig. 12). We can understand this behavior in terms of the massdependence of the horseshoe width. It has been shown that, when ignoring atmospheric feedback, 2 the half-width of the widest horseshoe orbit scales as 𝑥 hs ∝ √ 𝑚 (Paardekooper & Papaloizou 2009; Ormel 2013; Masset & Benítez-Llambay 2016). The corresponding streamline - the critical streamline - passes through the stagnation point at (± 𝑥 sep , 0 ) . For barotropic flows, the Bernoulli constant 𝐵  \n𝐵 = 1 2 𝑣 2 g + 𝑊 + Φ p + 𝜔 ∞ Ω 𝑥 2 (48)  \nis conserved along streamlines (Ormel 2013), where 𝑊 is the enthalpy, Φ p is the planet potential and 𝜔 ∞ = -3 2 Ω is the vorticity in the far field (far from the planet). Under the adiabatic EOS where 𝑃 = 𝜌 𝛾 g (in dimensionless units, see Sect. 2.5), the enthalpy can be calculated as,  \n𝑊 = ∫ d 𝑃 𝜌 g = ∫ 𝑃 𝜌 g d log 𝑃 d 𝑟 d 𝑟 = 𝑚 𝑟 (49)  \nafter substituting the adiabatic solution for a hydrostaticallysupported atmosphere, Eq. (49). Therefore, Φ P + 𝑊 ≈ 0 at the separatix location as well as in the far field. Equating the Bernoulli expression in the far field and at the separatix, respectively, gives:  \n1 2 ( 3 2 𝑥 hs ) 2 -3 2 𝑥 2 hs = -3 2 𝑥 2 sep . (50)  \nHence, 𝑥 sep ∝ 𝑥 hs ∝ √ 𝑚 . And 𝑥 sep is half of the distance of two stagnation points, which denotes the long axis of the atmosphere  \nFigure 13. Characteristic outputs of runs of different 𝑓 p2g , ini evolving with the simulation time. Left panel is like in Fig. 10; Right panel is similar to that in Fig. 10 while the right y-axis (red) represents pebble-p2g0.1 .  \n<!-- image -->  \n( 𝑅 atm denotes the short axis of the atmosphere). If we further assume that 𝑥 sep scales with 𝑅 atm, then 𝑅 atm / 𝑅 Bondi ∝ 𝑚 -1 / 2 accordingly. This means that with increasing planet mass 𝑚 , 𝑅 atm / 𝑅 Bondi will decrease, which explains the trend of initial values of 𝑅 subl / 𝑅 atm as shown in Fig. 12. Furthermore, the gravity at the sublimation front is -𝑚 / 𝑅 subl ∝ 𝑚 -1 . With a larger planet mass, pebbles go through a smaller gravitational focusing near the sublimation front, leading to the behavior closer to the perfectly-coupled runs, which also implies a vapour fraction similar to the perfectly-coupled runs. Therefore, the run pebble-m0.2 has a lower peak vapour fraction compared with the fiducial run (Fig. 12).  \nFor the run pebble-m0.05 , its initial 𝑅 atm / 𝑅 Bondi becomes higher, which reaches the sublimation radius and, 𝑅 subl ≈ 𝑅 atm. Instead of quickly reaching a steady state (recycling-dominated) or entering rapid vapour growth by the positive feedback, it instead reaches a quasi-steady state where all the output quantities evolve slowly due to the comparable rate of sublimation and freeze-out. This 'intermediate' regime has its 𝑅 subl gradually shrink due to the cooling from ice sublimation while 𝑅 atm keeps fluctuating around a value of 1.1 𝑅 Bondi . Compared to higher-mass runs, it has significantly higher vapour amounts in the atmosphere, but it is still not high enough to expand the atmosphere (the vapour fraction reaches ∼ 0 . 2 before 𝑅 atm begins to expand in pebble-5au ). Consequently, the atmosphere does not become significantly heavier, in contrast to the vapour-dominated runs (Sect. 3.4). For its flow pattern, near the Bondi radius (Fig. 9), there is a region where the velocity becomes very low (the arrows almost vanish). Due to this low velocity, the critical streamline is not well defined, causing 𝑅 atm to fluctuate strongly and the horseshoe orbit to become narrower. As a result, the gas near the core seems to be isolated from the horseshoe orbit.  \nIn order to find the reason for the low amount of circulation, we inspect the time evolution of this run. Initially, the flow pattern is similar to the recycling-dominated cases. However, a positive radial gradient of the vapour fraction emerges due to vapour piling up near the sublimation front. In addition, an azimuthal imbalance (both in  \ntemperature and pressure) also emerges caused due to ice preferentially sublimating near the separatix points. This triggers a strong compositional instability (see Sect. 4.2 for more discussions), which transports the piling-up vapour towards the core and excites pressure waves that reshape the flow pattern. In contrast, this instability does not happen in vapour-dominated runs because the atmosphere already expand significantly and the vapour pile-up happens deep inside the atmosphere. It also does not happen in the recycling-dominated runs as, there is no significant vapour pile-up. After this instability, the flows at the horseshoe orbit seem to experience additional pressure support contributed by the vapour, which weakens the horseshoe orbit. The high vapour fraction does not weaken the horseshoe orbit in the vapour-dominated runs, because the sublimation front resides deep inside the atmosphere, far from the horseshoe orbit (Fig. 9, pebble-5au e.g.). The additional pressure support contributed by vapour is tentatively seen in Fig. 11 a, where there is a dip of velocity profile near the sublimation front.  \nThe fact that the vapour fraction drops with increasing planet mass may seem counterintuitive. We also observe this phenomenon for planets at 4.5au (see Table 3) where the pebble-4.5au is classified as intermediate and the higher mass pebble-m0.2-4.5au as recycling-dominated. To further demonstrate this point, we run a simulation starting from 𝑡 = 300 Ω -1 of the vapour-dominated pebble-m0.05 run and let the planetary core grow linearly within 300 Ω -1 to 𝑚 = 0 . 2. In Fig. 12 (grey line), the peak vapour fraction drops significantly and finally relaxes to the same value of pebble-m0.2 . Similarly, the total vapour mass first increases with 𝑅 subl expanding sharply due to gravitational energy release. Then, after the core mass reaches 𝑚 = 0 . 2 ( 𝑡 = 600 Ω -1 ), 𝑀 vap begins to decrease, with recycling flows gradually removing vapour from the atmosphere.  \nIn summary, we found that with increasing mass, planet's atmospheres, while becoming bigger, may hold less (bound) vapour and that (sometimes) our qualitative assessment - recycling, intermediate, or vapour-dominated - may shift towards the former. The reason  \nFigure 14. Characteristic outputs of runs of different St evolving with the simulation time. Similar to Fig. 10 with only one y-axis for the right panels.  \n<!-- image -->  \nis that, although 𝑅 atm increases with 𝑚 , the sublimation radius increases faster.", '3.6 Varying the pebble-to-gas ratio': 'Next, we vary the initial pebble-to-gas ratio while other parameters are kept the same as in the fiducial run (see pebble-p2gxxx in Table 1). The flow patterns of pebble-p2g0.01 and pebble-p2g0.1 are shown in Fig. 9 along the purple line. The run pebble-p2g0.1 is still evolving by the time the snapshots is taken ( 𝑡 = 300 Ω -1 ). Figure 13 presents the time evolution of the characteristic output quantities.  \nVarying the pebble-to-gas ratio alters the pebble flux. With a higher pebble flux, more ice sublimates per unit time and more energy (latent heat) is absorbed from the environment. This mitigates the increase of the temperature caused by the gravitational energy release. Therefore, higher 𝑓 p2g , ini results in smaller sublimation radii initially (Fig. 13, 𝑡 ∼ 4 𝑡 soft ). For the run with 10% pebble-to-gas ( pebble-p2g0.1 ), the positive feedback effect (Sect. 3.4) results in continual atmosphere expansion and simultaneous shrinkage of the sublimation front. On the other hand, for the low pebble-to-gas ratio runs ( pebble-p2g0.01 and pebble-p2g0.02 ), 𝑅 subl > 𝑅 atm initially. Ice contained in pebbles is fully recycled and the vapour content reaches a steady value - the recycling-dominated limit (see Sect. 3.4).  \nFor intermediate values of the pebble-to-gas ratio (0.05), 𝑅 subl ≈ 𝑅 atm initially ( 𝑡 ∼ 4 𝑡 soft ) and the intermediate regimes emerges, which has its vapour content slowly grow and 𝑅 atm fluctuate more widely. Its flow pattern (not shown in Fig. 9) is also characterized by weak amount of circulation near the horseshoe and isolated inner atmosphere, similar to the 𝑚 = 0 . 05 run described above.  \nIn summary, the evolutionary trend of the vapour content is again determined by the respective values of 𝑅 atm and 𝑅 subl . However, instead of shrinking the sublimation front directly by moving the planet to colder regions (Sect. 3.4), increasing the pebble-to-gas  \nratio shrinks the sublimation front as a higher rate of sublimation cools the atmosphere. By varying the pebble-to-gas ratio (from 0.01 to 0.1), we observe a complete transition from recycling-dominated flows to the intermediate and finally the vapour-dominated regimes.', '3.7 Varying the Stokes number': 'The default Stokes number St = 0 . 01 represents cm-sized ( ≈ 0.9 cm) pebbles at 4 au. We vary them by one magnitude (St = 0.003 and 0.03) (see Figs. 9 and 14). The effect of varying the Stokes number is similar to varying the pebble-to-gas ratio, both of which will alter the pebble flux. In figure Fig. 14, the initial 𝑅 subl is smaller with higher Stokes number, which is caused by the higher pebble flux. Pebbles with smaller Stokes number are more easily recycled and the recycling efficacy will be the maximum for the perfectly-coupled pebbles (St = 0 particles), which sets the lower limit of the peak vapour fraction.  \nFor higher Stokes number runs ( pebble-St0.03 ), the intermediate regime emerges with its characteristic weak horseshoe orbits (Fig. 9). Although 𝑅 atm shows significant expansion, the ratio of the 𝑅 subl and 𝑅 atm is not that extreme compared to vapour-dominated runs (Figs. 10 and 13), thus the positive feedback does not emerge to promote further growth of vapour and expand atmosphere. Its sublimation front expands gradually after 𝑡 ≈ 300 Ω -1 , showing a converging trend of 𝑅 atm and 𝑅 subl at late times, which causes the vapour content to decrease. This process is similar to the positive feedback discussed in Sect. 3.4 but in the opposite direction: net ice freezes-out heats the sublimation front region, causing it to expands. However, the vapour fraction is still too small to alter 𝑅 atm significantly, rendering this process much weaker.  \nFigure 16. Predictions of 𝑅 atm and 𝑅 subl for different disc temperatures and materials of interest. Colored symbols denote different materials labelled in the figure while the black star represents water ice under the disc conditions of Johansen et al. (2021) at 𝑡 = 5 Myrs. The black solid line denotes 𝑅 subl = 𝑅 atm.  \n<!-- image -->  \nFigure 15. Summary of the characteristic outpus of all runs on the ( 𝑅 atm / 𝑅 Bondi , 𝑅 subl / 𝑅 Bondi ) plane. The grey lines connect different outputs of the same run with the filled, edge-free symbols denoting the outputs at 𝑡 = 4 𝑡 soft (initial), the filled, solid-edged symbols denoting the outputs after evolving time ( 𝑡 evolv ) and the open, dashed-edged symbols denoting the predicted value for perfectly-coupled grains (Sect. 4.1). Black nodes are also put on the initial outputs for highlighting. Different symbols denote which parameters are varied with respect to the fiducial run (see Table 1) with the opacity indicating different values in the sequence. For example, increasing the opacity of the circles indicates increasing dust-to-gas ratios (0.05, 0.1, and 0.2). The fiducial run is indicated by the cross and run pebble-m0.2-4.5au by a star.  \n<!-- image -->', '4.1 The composition of planets and atmospheres': "In Fig. 15 we summarize all runs listed in Table 1 on the ( 𝑅 atm / 𝑅 Bondi , 𝑅 subl / 𝑅 Bondi ) plane. Symbols represent a certain parameter while the opacity of the color indicates the values that are varied for this parameter. The run pebble-m0.2-4.5au , where two parameters are varied with respect to the default, is indicated by the star. For each run three kind of symbols are plotted: open, dashed edge; filled, no edge; and filled, solid edge. The filled, edge-free symbols (also with black nodes inside) represent the initial state of the runs ( 𝑡 = 4 𝑡 soft ), which are clearly separated by the line 𝑅 subl = 𝑅 atm: the recycling-dominated runs (blue; above the line), the intermediate runs (purple; around the line) and the vapour-dominated runs (red; below the line). The evolved state (filled, solid edge) plots the outputs after time 𝑡 evolv (Table 3, 150 Ω -1 for the recycling-dominated runs and 300 Ω -1 for the intermediate and vapour-dominated runs). For the recycling-dominated runs, both 𝑅 subl and 𝑅 atm hardly change and the symbols tend to overlap. On the other hand, for the vapourdominated runs, the symbols have 𝑅 atm and 𝑅 subl separate divergently ( 𝑅 atm increases, 𝑅 subl decreases) through the positive feedback mechanism explained in Sect. 3.4. The intermediate runs move awayfromthe 𝑅 subl = 𝑅 atm line mostly by the decrease of 𝑅 subl while 𝑅 atm expands less significantly due to the limited accumulation of vapour (see Sect. 3.7 for details).  \nFinally, the open, dashed-edged symbols in Fig. 15 indicate our 'predicted' values for 𝑅 subl and 𝑅 atm. In Sect. 3.5, we already derived 𝑅 atm / 𝑅 Bondi ∝ 𝑚 -1 / 2 and we fit the prefactor as,  \n𝑅 atm 𝑅 Bondi = 0 . 25 𝑚 -1 / 2 . (51)  \nAssuming that sublimation happens in a very narrow radial range and the adiabatic solution still fits well (Eqs. (46) and (47)), 𝑅 subl can be numerically solved given 𝑓 vap , peak , disc parameters ( 𝑇 disc , 𝑃 disc ) and material constants (see Table 2 and App. B). Here we take 𝑓 vap , peak ≈ 𝑓 ice , ini , where 𝑓 ice , ini is the initial ice fraction ( 𝑓 ice , ini = 𝑓 p2g , ini 𝑓 𝑧 , see Sect. 3.2 and Table 1). This approximation is of course only appropriate for the tightly coupled runs. Also, we assume 𝛾 can be taken a constant. The prediction for recycling-dominated runs fits the simulation results well, because the vapour fraction stays low. For intermediate and vapour-dominated runs, the prediction tends to overestimate 𝑅 subl because the vapour fraction at 𝑡 = 4 𝑡 soft is already much higher than that in the perfectly-coupled runs, which also causes the temperature profile to deviate from the adiabatic solution of pure H-He gas (Eq. (47)).  \nWhile we could improve the prediction by correcting it for, e.g., the Stokes number, the trend is never the less clear. Conditions where the predicted value yields 𝑅 subl / 𝑅 atm ≫ 1 will turn out to become recycling-dominated, while for 𝑅 subl / 𝑅 atm ≪ 1 the planet becomes vapour-dominated. We expect this result to apply for different disc conditions and for different materials.  \nAkey application is that of water delivery to Earth. With the waning of the viscous accretion rate as well as the decline of the Sun's luminosity on its descent to the pre-main sequence, it is believed that the disc H 2 O iceline would have crossed Earth's orbit (Garaud & Lin 2007; Oka et al. 2011; Baraffe et al. 2015; Ida et al. 2016). This would have allowed proto-Earth to accrete H 2 Oice-rich pebbles and turn wet, inconsistent with its low water content (2000 -3000ppm Marty 2012). Recently, Johansen et al. (2021) argued that recycling flows prevented accumulation of too much water. In Fig. 16 the point corresponding to the conditions in their simulations 3 is shown by the  \n3 We adopt similar disc parameters ( 𝑇 disc = 130K, 𝜌 disc ≈ 10 -11 g cm -3 , 𝑓 ice , ini ≈ 0 . 0013) at 5Myrs for a proto-earth ( ∼ 0 . 6 𝑀 ⊕ , before the giant impact to produce moon) at 1 au, which is just before the disc dissipation time assumed in Johansen et al. (2021).  \nblack star. As this point lies above the 𝑅 subl = 𝑅 atm line, a recyclingdominated atmosphere is a possibility. If recycling-dominated or intermediate, we expect the vapour contents in the atmosphere to stabilize and even at a level of ≲ 10% (characteristic values for the intermediate regime, see Table 3), the contribution to the bulk composition would still be tiny ( ∼ 10ppm), given that the characteristic envelope mass of terrestrial planets is ∼ 10 -4 𝑀 ⊕ (Johansen et al. 2021). However, a transition to an H 2 O vapour-dominated envelope (and planet) - inconsistent with Earth- can also not be ruled based on our results.  \nThe estimator can also be applied towards materials other than water. Based on the disc model adopted in this work (Sect. 2.5), we investigate retention of carbon and silicates of an accreting Earthmass planet at 1 au. In our solar system, the main carbon carrier that contributes ∼ 80% of the total amount in solids are complex organics which sublimate at ∼ 325 -425 K (Gail & Trieloff 2017). As a proxy, we take the properties of pyrene, one of the common Polycyclic Aromatic Hydrocarbons (PAHs), (see Table 2). We further assume that pebbles roughly contain 5% of their mass in carbon, similar to the CI chondrites (Johansen et al. 2021). For the silicate, we take olivine as the proxy (King et al. 2015) and assume that it accounts for 50% of the pebble mass. The pebble-to-gas ratio is set as in the fiducial run. To account for the changing luminosity of the Sun, we vary the disc temperature from 130 to 300 K (see Fig. 16) and the disc pressure accordingly. With 𝑇 disc increasing, 𝑅 atm / 𝑅 Bondi increases due to the increase of local thermal mass, which lowers 𝑚 with the same planet mass. Except for the lowest disk temperatures, water will be recycled, resulting in a dry planet. To retain the carbon, 𝑇 disc needs to significantly decrease (below ∼ 200 K). Hence, Earth's low-carbon fraction (Marty 2012) may indicate the disc stayed above this temperature. In contrast, silicates always fall in the vapour-dominated regime, implying injection of silicate vapour deep into the atmosphere. This 'vapour-planet' scenario has already been investigated by 1D analytical (Brouwers & Ormel 2020) and numerical (Ormel et al. 2021) modelling highlighting the formation of a ∼ 1 𝑀 ⊕ , vapour-rich atmosphere for sub-Neptune planets. Our 2D work suggest that silicate recycling flows are ineffective to alter these findings.", '4.2 Caveats': "Our pioneering work has adopted several simplifications which may affect our results: simulations were conducted in 2D without accounting for a disc headwind, did neither include self gravity nor radiation transport (RT), and did not consider frictional heating of the gas by the accreted pebbles (dust feedback). The omission of the latter two compensate each other to some extent, such that a (quasi-)steady state is often still reached. In our simulation ommiting the frictional heating, higher pebble flux leads to the vapor-dominated regime because the cooling by ice sublimation pushes the 𝑅 subl inwards (as illustrated in Sect. 3.6 & 3.7). However, in reality, higher pebble flux will increase the temperature of the atmosphere by depositing accretion energy of refractories throughout the atmosphere (Pollack et al. 1996; Rafikov 2006; Zhu et al. 2021). This would probably render the trend in the opposite direction: with higher pebble flux, 𝑅 subl moves outwards and the atmosphere ends up in the recycling-dominated regime. With RT included, the outer regions of pre-planetary atmospheres are likely to be (nearly-)isothermal on top of an optically thick interior - either radiative or convective (Pollack et al. 1996; Inaba &Ikoma 2003; Rafikov 2006; Ormel & Kobayashi 2012; Mordasini et al. 2012). The presence of an isothermal outer layer, therefore, would shift the 𝑅 subl inwards, by an amount that depends on the opacity of the gas or that of the pebbles and recondensed ices (Ormel  \net al. 2021; Brouwers et al. 2021). In contrast, the headwind induced by the disc pressure gradient could possibly shift 𝑅 atm inwards, depending on the disc structure and planet mass (Ormel 2013; Ormel et al. 2015b; Kurokawa & Tanigawa 2018; Moldenhauer et al. 2022).  \nIncluding RT may affect the flow instability seen in the 'intermediate' runs. In these simulations ice preferentially sublimated near the separatix points, which resulted in an azimuthal imbalance (in temperature and pressure). It needs to be verified whether the ensuing instability will also appear in simulation that include RT or whether the simulations presently classified as 'intermediate' will evolve into either the recycling- or vapour-dominated regimes. Interestingly, the interior region in these simulations reveals a weak horseshoe motion similar to that of the 𝛽 -cooling approach adopted by Kurokawa & Tanigawa (2018), suggesting that this low-entropy material is hydrodynamically shielded from the recycling-dominated flows - the buoyancy barrier. In our case, the entropy gradient caused by ice sublimation possibly accounts for the buoyancy barrier. But here, again, RT may act to at least partially erase these features (Moldenhauer et al. 2022). In addition, including RT is likely to ameliorate the entropy violation issue that we occasionally witnessed in our simulations (Sect. A and 3.3).  \nThree-dimensional (3D) flows are known to be qualitatively different from their 2D counterparts. In 3D, material tends to flow in through the polar region and flow out in the midplane (equatorial flows; Ormel et al. 2015b; Fung et al. 2015; Szulágyi et al. 2016; Cimerman et al. 2017; Moldenhauer et al. 2022). In 2D, simulations tend to be characterized by large rotational flows around the core (Ormel et al. 2015a; Béthune & Rafikov 2019 or see e.g. Fig. 11), whereas the amount of rotation in 3D simulations is much weaker due to the polar inflow (Ormel et al. 2015b; Cimerman et al. 2017) at least for low-mass planet and non-isothermal EOS (Fung et al. 2019). This may lead to more efficient mixing, smoothing out the radial gradient of vapour fraction, which again questions the emerge of instabilities seen in the intermediate runs. In 3D, Moldenhauer et al. (2021) showed a fully-recycled atmosphere on a timescale of ∼ 10 3 orbits. But whether the same picture emerges in 3D for sublimated pebbles is still to be investigated. We envision a similar dichotomy as seen in this work with planets accreting modestly volatile pebbles at low rates attaining a steady state and operating in the recycling limit, whereas more refractory materials will be trapped in the atmosphere ( 𝑅 subl / 𝑅 atm ≪ 1).  \nAnother key assumption we adopted is that recondensed ice grains flow back to the disc together with the gas. However, solid recycling flows may be arrested when the condensed particles become too large to settle back to the planet. Whether or not grains will leave the atmosphere depends on how fast ice grains grow into larger particles. The coagulation timescale of grains is (Ormel et al. 2007; Krĳt et al. 2016; Johansen et al. 2021)  \n𝜏 coag ∼ 𝑠 ice 𝜌 · , ice Δ 𝑣𝜌 ice (52)  \nwhere 𝑠 ice is the radius of ice grains, 𝜌 · , ice is its internal density, Δ 𝑣 the particle collision speed and 𝜌 ice the mass density of ice grains. Assuming that the relative velocity is dominated by Brownian motion (small particles ≲ 1 𝜇 m), we obtain,  \n𝜏 coag ∼ 𝑠 ice 𝜌 · , ice √︁ 𝑘 B 𝑇 / 𝑚 p 𝜌 g 𝑓 -1 ice ≈ 47 . 5 yr × GLYPH<18> 𝑠 ice 1 𝜇𝑚 GLYPH<19> 5 2 GLYPH<18> 𝑇 135K GLYPH<19> -1 2 GLYPH<18> 𝜌 · 1 g cm -3 GLYPH<19> 3 2 GLYPH<18> 𝜌 g 10 -10 g cm -3 GLYPH<19> -1 GLYPH<18> 𝑓 ice 0 . 01 GLYPH<19> -1 (53)  \nwhere 𝑓 ice is the ice-to-gas density ratio. In the above expression, we have inserted the values typical for the fiducial run. In addition, we assumedice grains have grown to 1 𝜇 m,whichis approximately the size wheregrowthbyBrownianmotionisreplacedbydrift-inducedmechanisms (Dullemond & Dominik 2005; Ormel & Okuzumi 2013). This time must be compared to the recycling timescale 𝜏 rec. If the sublimation radius lies around or outside the Bondi radius, the outflow motion will have a velocity similar to the shear velocity, 𝑣 rec ∼ Ω 𝑅 subl . Therefore,  \n𝜏 rec ∼ 𝑅 subl 𝑣 rec ∼ Ω -1 ∼ 1 yr (54)  \nAs the coagulation timescale is much longer than the recycling timescale there will be no significant growth of ice grains during their transport away from the sublimation front. Therefore, pebble recycling will not be influenced by dust coagulation unless the coagulated grains have a very fractal structure (low 𝜌 · ; Ormel 2014).  \nIn contrast, for the vapour-dominated regime ice grains accumulate deep in the atmosphere (see Fig. 9) and the coagulation timescale is much shorter. We have measured the coagulation timescale for the pebble-5au run at 𝑡 = 500 Ω -1 , when 𝑓 ice ≈ 0 . 2, 𝑇 ≈ 200 K and 𝜌 g ≈ 10 -10 g cm -3 (see Fig. 11), as 𝜏 coag ∼ 1 Ω -1 . In addition, inside the critical streamline, where the streamlines are closed, the recyling timescale is given by the diffusion timescale 𝜏 rec ∼ 𝑅 2 atm / 𝐷 g ∼ 10 4 Ω -1 . Since this coagulation timescale is much shorter than the recycling timescale, ice grains will grow to larger particles in the atmosphere and experience strong gravitational settling, causing them to fall back to the planet and the amount of grains that escape will be much lower. Instead of an extended ice grain distribution pervading the whole atmosphere seen in Fig. 9, coagulation would limit the particles more to the sublimation front region. This would further exacerbate the growth of the planet's vapour content.", '4.3 Applications to disc chemistry': 'In this work we have conducted local simulations, focusing on the consequences of sublimation for the putative planet atmosphere and its bulk composition. In the recycling-dominated regime the sublimated pebbles simply escape the simulation domain in the form of ice grains by instaneous freezing-out and there is no further interaction between the outflow vapour and the disc. However, the outflow vapour could possibly spread azimuthally through horseshoe orbits to fill the entire co-orbital region given the finite freeze-out rate of vapour (Minissale et al. 2022). Pebble recycling can therefore change the chemical abundance of disks locally.  \nIn particular, the high accretion rates characterizing accreting giant planet would elevate ambient temperature significantly (Machida et al. 2010; Szulágyi 2017; Bergez-Casalou et al. 2020). 𝑅 subl could even expand to ∼ 𝑅 Hill , releasing vapour over a much wider region ( ∼ au) that is observationally accessible. Moreover, giant planets (Lin &Papaloizou 1979; Kley 1999; Kanagawa et al. 2017) tend to open a deep gap, which allows the penetration of high flux UV-photons and promote further photo-chemistry (Bergin et al. 2016; Miotello et al. 2019; Bosman et al. 2021). Cleeves et al. (2015) proposed that the asymmetric features seen in molecular line emission might imply the existence of young giant planets. An example is the observed timeevolving SO asymmetry in HD 100546 (Booth et al. 2023), which is proposed to trace the post shock region where the warm ( ∼ 100 K) environment caused by accreting giants is suitable for photo-chemistry to form SO.  \nQuantitatively, these questions can best be addressed by global  \nsimulations. In a companion paper (Jiang et al.) we have implemented the sublimation module approach accordingly to address the potential of recycling of C-rich ices. This could modify the chemical abundance (e.g., C/O ratio) in the co-orbital region of the planet, which could result in distinct substructures as seen, for example, in the Molecules with ALMA at Planet-forming Scales survey (MAPS, Öberg et al. (2021)). Our phase change module can also be augmented with a chemical network to further understand the connection between young planets and their natal discs.', '5 CONCLUSIONS': "We have developed a phase change module based on the multi-fluid dust module in Athena++. This module treats the mass transfer, energy exchange and momentum conservation processes during sublimation and condensation (freezing-out) self-consistently, enabling us to study the coupled thermodynamic and hydrodyanmical effect of phase change properly in numerical simulations. In this work, we focus on investigating how ice sublimation influence the properties of the planetary atmosphere and the resulting vapour content through pebble accretion.  \nOur main findings are:  \n- (i) The long-term evolutionary trend of the atmospheric vapour content is determined by the relative size of the sublimation radius 𝑅 subl and the atmosphere radius 𝑅 atm. When 𝑅 subl exceeds 𝑅 atm, incoming pebbles can be fully recycled and the vapour content reaches a steady value - the recycling-dominated limit. On the other hand, When 𝑅 subl < 𝑅 atm, vapour tends to be locked deep in the atmosphere and keep accumulating - the vapour-dominated limit.\n- (ii) The accumulation of vapour results in a positive feedback. 𝑅 atm expands due to material accumulation while 𝑅 subl shrinks due to the cooling caused by pebble sublimation. The divergence of 𝑅 subl and 𝑅 atm suppresses recycling and promotes further accumulation of vapour.\n- (iii) Planets locating further out with respective to the star (lower ambient temperature) and higher pebble flux induced by either higher pebble-to-gas ratio or larger Stokes number tend to shrink the sublimation front, resulting in vapour-dominated regime. Conversely, with lower pebble fluxes or higher ambient temperatures, simulations lean more towards the recycling-dominated limit. However, the trend with pebble flux will probably reverse when frictional heating by pebbles is included (Sect. 4.2).\n- (iv) With increasing planet mass, 𝑅 subl increases faster than 𝑅 atm. Consequently, the classification - recycling, intermediate, or vapourdominated - tends to shift towards the former. Planet's atmospheres, while becoming bigger, may hold less (bound) vapour.\n- (v) The intermediate regime is a transition between the recyclingdominated and the vapour-dominated cases. It features an instability triggered by composition gradients near the atmosphere radius and simulations tend to be characterized by weakened horseshoe orbits.\n- (vi) We develop estimators for 𝑅 subl and 𝑅 atm based on disc conditions and material properties. Applying these to an Earth-mass planet at 1 au in an MMSN disc, we find it to be recyling-dominatant in water vapour (i.e., a dry planet) while silicates are always in the vapour-dominated regime. Preventing carbonaceous material to be accumulated, requires the ambient temperature to be ≳ 200 K.  \nThe phase change module developed by us is general. In the future, it can be applied to address questions regarding the chemical compositions influenced by the planet-disc interaction and planetesimal formation near the snowline, among others.", 'ACKNOWLEDGEMENTS': 'The authors thank the referee for their positive assessment of this manuscript. The authors appreciate fruitful discussions with Zhaohuan Zhu, Allona Vazan, Xuening Bai, Zhuo Chen, Kees Dullemond, and Mario Flock. The authors acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper. CWO acknowledges support by the National Natural Science Foundation of China (grant no. 12250610189). RK acknowledges financial support via the Heisenberg Research Grant funded by the German Research Foundation (DFG) under grant No. KU2849/9.', 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'Szulágyi J., Masset F., Lega E., Crida A., Morbidelli A., Guillot T., 2016,  \nSchoonenberg D., Ormel C. W., 2017, A&A, 602, A21 Szulágyi J., 2017, ApJ, 842, 103 MNRAS, 460, 2853  \nVazan A., Ormel C. W., Dominik C., 2018a, A&A, 610, L1  \nVazan A., Ormel C. W., Noack L., Dominik C., 2018b, ApJ, 869, 163  \nVenturini J., Alibert Y., Benz W., Ikoma M., 2015, A&A, 576, A114  \nVenturini J., Alibert Y., Benz W., 2016, A&A, 596, A90  \nWeidenschilling S. J., 1977, MNRAS, 180, 57  \nWinn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409  \nWuchterl G., 1993, Icarus, 106, 323  \nYoudin A. N., Lithwick Y., 2007, Icarus, 192, 588  \nYu C., 2017, ApJ, 850, 198  \nZhang S., et al., 2018, ApJ, 869, L47  \nZhu W., Dong S., 2021, ARA&A, 59, 291  \nZhu W., Petrovich C., Wu Y., Dong S., Xie J., 2018, ApJ, 860, 101  \nZhu Z., Jiang Y.-F., Baehr H., Youdin A. N., Armitage P. J., Martin R. G., 2021, MNRAS, 508, 453  \nvan der Marel N., Dong R., di Francesco J., Williams J. P., Tobin J., 2019, ApJ, 872, 112', 'APPENDIX A: GAS-ONLY RUN': 'In this part, we consider the situation where an inviscid gas flow passes through a planet atmosphere, a problem which has already been studied before (Ormel et al. 2015a; Fung et al. 2015; Kurokawa & Tanigawa 2018; Béthune & Rafikov 2019). The parameters of these gasOnly runs are listed in Table 1. By design, the gas is adiabatic and 𝛾 = 1 . 43 (assume 71% H 2 and 29% He). We also test cases with 𝛾 = 1 . 3 , 1 . 1 (other parameters are the same), see Fig. A1. The first three panels show steady state flow patterns of cases with 𝛾 = 1 . 43 , 1 . 3 and 1 . 1 respectively, where the color denotes gas density (on logarithmic scale).  \nOverall the flow patterns show the typical horseshoe orbit and inner rotational atmosphere shown in previous studies. With the decrease of the adiabatic index, the atmosphere expands and rotational support becomes more and more important, which is clearly shown in the upper panel of Fig. A2. Also, the density increases significantly. As 𝛾 become smaller, pressure support becomes less effective, resulting in more gas falling in, which is also seen when comparing the adiabatic and isothermal atmospheres in Béthune & Rafikov (2019).  \nFor an adiabatic EOS, the entropy should be conserved, which is approximately true for high 𝛾 runs. However, for the run with 𝛾 = 1 . 1, the entropy drops steeply to nearly 0 at the inner boundary, concommitant with the abrupt increase of gas density. We verified that no energy flux was flowing through the inner boundary consistent with the reflective boundary condition. We suspect that the unphysical entropy violation is caused by the large azimuthal velocity that reaches Keplerian values near the inner boundary, which potentially induces significant numerical viscosity. We verified that mitigating velocity gradients by adopting gravitational softening (Fung et al. 2019; Zhu et al. 2021) or a larger value of the inner radius or higher resolution, suppresses the artificial entropy variation, while solving for the entropy equation instead of the energy equation (Mignone et al. 2007, 2012) eliminated the issue. In the simulations runs in the main text we simply monitor the entropy evolution, and re-run simulations at higher resolution when the entropy variation exceeds ∼ 10% and begins to alter the flow pattern.', 'APPENDIX B: ESTIMATION OF SUBLIMATION FRONT': "Combining Eqs. (46) and (47),  \n𝑓 vapr , peak GLYPH<20> 1 + 𝛾 -1 𝛾 𝑅 Bondi GLYPH<18> 1 𝑅 subl -1 𝑟 max GLYPH<19> GLYPH<21> 𝛾 𝛾 -1 𝑃 disc = 𝑃 eq , 0 exp -𝑇 a 𝑇 ' subl ! = ⇒ log GLYPH<18> 𝑓 ice 𝑃 disc 𝑃 eq , 0 GLYPH<19> + 𝛾 𝛾 -1 log 𝜒 = -𝑇 a 𝑇 ' subl (B1)  \nwhere 𝜒 = 1 + 𝛾 -1 𝛾 𝑅 Bondi GLYPH<16> 1 𝑅 subl -1 𝑟 max GLYPH<17> . 𝑇 ' subl is the modified temperature at the sublimation front, which takes into accout the sublimation cooling compared to 𝑇 subl in Eq. (47). Assuming all the ice sublimates at 𝑅 subl (In reality sublimation happens in a finite band, see Fig. 8), latent heat should be deducted from the gas internal energy,  \n𝑃 g 𝜌 g ( 𝛾 -1 ) -𝐿 heat 𝑓 vap , peak = 𝑘 B 𝑇 ' subl 𝜇 g 𝑚 p ( 𝛾 -1 ) (B2)  \nFor 𝑃 g and 𝜌 g we still take the adiabatic solution, which is a very rough zero-order approximation. Then 𝑇 ' subl yields,  \n𝑇 ' subl = 𝑇 subl -𝐿 heat 𝑓 vap , peak ( 𝛾 -1 ) 𝑘 B /( 𝜇 g 𝑚 p ) (B3)  \nCombining Eqs. (B1) and (B3), 𝑅 subl can be numerically solved.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A1. Steady state flow patterns of gasOnly case with different adiabatic index 𝛾 . The color denotes gas density in log scale. Their x and y axies are all plotted in the unit of local scale height. The black dashed line in each plot denotes Bondi radius (Eq. (36)) and the blue solid line is 'critical streamline' (Sect. 3.2).  \n<!-- image -->  \nFigure A2. Comparison of azimuthally averaged radial density and velocity profiles of run gasOnly . The velocity profile is in unit of the planet's Keplerian velocity 𝑣 k = √︁ 𝑚 / 𝑟 and the density profile is normalized by the disc's unperturbed value 𝜌 0.  \n<!-- image -->"}