We present nestedgrid highresolution hydrodynamic simulations of gas and particle dynamics in the vicinity of Mars to Earthmass planetary embryos. The simulations extend from the surface of the embryos to a few vertical disc scale heights with a spatial dynamic range of 1.4  10SUP5SUP. Our results confirm that pebblesized particles are readily accreted with accretion rates continuing to increase up to metresize boulders for a 10 per cent MMSN surface density model. The gas mass flux in and out of the Hill sphere is consistent with the Hill rate   RH24 103 Moplus yr1. While smaller size particles mainly track the gas a net accretion rate of   2 105 Moplus yr1 is reached for 0.31 cm particles even though a significant fraction leaves the Hill sphere again. Effectively all pebblesized particles that cross the Bondi sphere are accreted. The resolution of these simulations is sufficient to resolve accretiondriven convection. Convection driven by a nominal accretion rate of 106 Moplus yr1 does not significantly alter the pebble accretion rate. We find that due to cancellation effects accretion rates of pebblesized particles are nearly independent of disc surface density. As a result we can estimate accurate growth times for specified particle sizes. For 0.31 cm size particles the growth time from a small seed is 0.15 million years for an Earthmass planet at 1 au and 0.1 million years for a Mars mass planet at 1.5 au.
['hydrodynamics (HD)', 'methods: numerical', 'planets and satellites: formation', 'terrestrial planets', 'protoplanetary discs', 'Astrophysics - Earth and Planetary Astrophysics']
Pebble dynamics and accretion on to rocky planets  I. Adiabatic and convective models
{'Pebble dynamics and accretion onto rocky planets. I. Adiabatic and convective models': 'Andrius Popovas, 1 ; 3 ? ˚ Ake Nordlund, 1 y Jon P. Ramsey, 1 z and Chris W. Ormel, 2  \n1 Center for Star and Planet Formation, the Niels Bohr Institute and the Natural History Museum of Denmark,  \nUniversity of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen, Denmark  \n- 2 Anton Pannekoek Institute, University of Amsterdam, Science Park 904, C4.149, 1098 XH Amsterdam, The Netherlands,\n- 3 Rosseland Centre for Solar Physics, Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway  \nAccepted 2018 June 27. Received 2018 June 14; in original form 2018 January 26', 'ABSTRACT': "We present nested-grid, high-resolution hydrodynamic simulations of gas and particle dynamics in the vicinity of Mars- to Earth-mass planetary embryos. The simulations extend from the surface of the embryos to a few vertical disk scale heights, with a spatial dynamic range of GLYPH<24> 1 : 4 GLYPH<2> 10 5 . Our results confirm that 'pebble'-sized particles are readily accreted, with accretion rates continuing to increase up to metre-size 'boulders' for a 10% MMSN surface density model. The gas mass flux in and out of the Hill sphere is consistent with the Hill rate, GLYPH<6> GLYPH<10> R 2 H = 4 10 GLYPH<0> 3 M GLYPH<8> yr GLYPH<0> 1 . While smaller size particles mainly track the gas, a net accretion rate of GLYPH<25> 2 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 is reached for 0.3-1 cm particles, even though a significant fraction leaves the Hill sphere again. Effectively all pebble-sized particles that cross the Bondi sphere are accreted. The resolution of these simulations is sufficient to resolve accretion-driven convection. Convection driven by a nominal accretion rate of 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 does not significantly alter the pebble accretion rate. We find that, due to cancellation effects, accretion rates of pebble-sized particles are nearly independent of disk surface density. As a result, we can estimate accurate growth times for specified particle sizes. For 0.3-1 cm size particles, the growth time from a small seed is GLYPH<24> 0.15 million years for an Earth mass planet at 1 AU and GLYPH<24> 0.1 million years for a Mars mass planet at 1.5 AU.  \nKey words: planets and satellites: formation, terrestrial planets - protoplanetary discs - hydrodynamics (HD) - methods: numerical", '1 INTRODUCTION': "We now know that the majority of stars in our galaxy host at least one planet (Cassan et al. 2012). We also know that planets form inside the gaseous and dusty protoplanetary disks (PPDs) that orbit stars during their youth. As it is now possible to spatially resolve the structure of young PPDs (e.g., with the Atacama Large Millimetre/sub-millimetre Array; ALMA), it is becoming apparent that the majority of PPDs have a rich sub-structure of rings, gaps, vortices and spirals (e.g. van der Marel et al. 2013; Isella et al. 2013; Casassus et al. 2013; ALMA Partnership et al. 2015; Andrews et al. 2016; Meru et al. 2017). Recent, spatially resolved observations of young disks with ALMA (ALMA Partnership et al. 2015; Andrews et al. 2016), as well as meteoritic evidence (Bizzarro et al. 2017), furthermore in-  \n- ? E-mail: popovas@nbi.ku.dk (AP)\n- y E-mail: aake@nbi.ku.dk ( ˚ AN)\n- z E-mail: jramsey@nbi.ku.dk (JPR)  \ndicate that planet formation starts early. What remains to be answered, however, is how planets grow efficiently enough under realistic conditions to agree with the astronomical observations and meteoritic evidence.  \nThere are currently two scenarios for planet growth which are popular in the literature: First, the planetesimal accretion scenario, in which a planetary embryo is impacted by kilometre-size bodies and their mass is added to the embryo (e.g. Pollack et al. 1996; Hubickyj et al. 2005). Eventually, the embryo becomes massive enough that gravitational focusing (Greenberg et al. 1978) becomes important and the embryo enters what is commonly called the 'runaway growth' phase. Planetesimal accretion ends when there is no remaining solid material in an embryo's 'feeding zone', e.g., when the embryo grows massive enough to gravitationally stir up nearby planetesimals, enhancing their eccentricities and effectively kicking them out of the feeding zone (Ida & Makino 1993). The growth of the, now somewhat isolated, protoplanet transitions to the much slower 'oligarchic growth' phase (Kokubo & Ida 1998). In this phase, the  \nfew, largest mass protoplanets grow oligarchically, while the remaining planetesimals mostly remain small. The critical time scale in this context is the lifetime of the PPD, which is of the order of several million years (e.g. Bell et al. 2013). However, numerical simulations of the planetesimal accretion hypothesis have the problem that they predict that it takes much longer than a PPD lifetime for cores to grow to observed planetary masses (e.g. Levison et al. 2010; Bitsch et al. 2015b).  \nThe'pebble accretion'scenario has, meanwhile, demonstrated promise in being able to accelerate the protoplanet growth process significantly (Ormel & Klahr 2010; Johansen & Lacerda 2010; Nordlund 2011; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012; Lambrechts & Johansen 2014; Bitsch et al. 2015b; Chatterjee & Tan 2014; Visser & Ormel 2016; Ormel et al. 2017, etc.). 'Pebbles', in the astrophysical context, are millimetre to centimetre-sized particles with stopping times, t s , comparable to their orbital period, t s GLYPH<24> GLYPH<10> GLYPH<0> 1 K , where GLYPH<10> K = p GM GLYPH<3> GLYPH<157> r 3 is the Keplerian frequency, G is the gravitational constant, M GLYPH<3> and r are the mass of and the distance to the central star. Pebbles are likely to form a significant part of the solid mass budget in PPDs, as indicated by both dust continuum observations (e.g. Testi et al. 2003; Lommen et al. 2009) and by the mass fraction of chondrules in chondritic meteorites (Johansen et al. 2015; Bollard et al. 2017). In this scenario, the combined effects of gas drag due to the difference in speed between the slightly sub-Keplerian dust (and gas) and the Keplerian embryo, and gravity ensure that pebbles are captured and settle to the planet. As the embryo grows, the efficiency of gravitational focusing of pebbles increases, and the effective accretion cross section becomes larger than the embryo itself.  \nThe radius of dominance of the gravitational force of a planet relative to the central star is approximately given by the Hill radius:  \nR H = a 3 s M p 3 M GLYPH<3> ; (1)  \nwhere a is the semi-major axis of the embryo's orbit, M p and M GLYPH<3> are the masses of the embryo and the central star, respectively. Particles of suitable size, passing the embryo even as far away as the Hill radius may be accreted, as has been shown analytically by Ormel & Klahr (2010), using test particle integrations on top of hydrodynamical simulations by Morbidelli & Nesvorny (2012) and using numerical simulations with particles by Lambrechts & Johansen (2012). According to these works, pebble accretion is so efficient that it can reduce planet growth time scales to well within the lifetime of a PPD, even at large orbital distances from the host star.  \nHerein, we present the first-ever high resolution simulations of gas and pebble dynamics in the vicinity of low-mass, 'rocky' planetary embryos (for simplicity, often referred to only as 'embryos' in what follows) and report on the measured accretion efficiency of pebbles. We consider three embryo masses: 0.95 M GLYPH<8> , 0.5 M GLYPH<8> , and 0.096 M GLYPH<8> , and we specifically choose to study conditions expected for pressure traps associated with inside-out-scenarios of planet formation (e.g. Tan et al. 2016). The embryos are thus assumed to be embedded in disks that have a local pressure maximum at the orbital radius of the embryo.  \nThe presence of a pressure bump instead of a flat density distribution provides three advantages:  \n- (i) It eliminates the head wind otherwise associated with the slightly sub-Keplerian motion of the gas;\n- (ii) It helps to trap particles drifting inwards from larger orbital radii (Whipple 1972; Paardekooper & Mellema 2006);\n- (iii) High resolution observations (e.g. ALMA Partnership et al. 2015; Andrews et al. 2016) show disks with rings and gaps, which affect the distribution of solids in the disk. Having the pressure bump allows particles to drift into the pressure maximum, thus mimicking the observed behaviour.  \nThe structure of the paper is as follows. In Section 2 we describe the simulation set up, including initial conditions and their dynamic relaxation to a quasi-stationary state. We then describe, in Section 3, our treatment of particles, their equations of motion and how they are injected into the simulation domain. In Section 4 we describe the gas dynamics that develops in our simulations and compare our results to other studies. In Section 5 we present our main results regarding particle dynamics and accretion efficiency onto planetary embryos. In Section 6 we present the first results of simulations of accretion-driven convection in the primordial atmospheres of rocky planets. Finally, in Section 7, we summarize our main results, discuss their implications and indicate future work.", '2 SIMULATION SET UP': "This study is carried out using the new DISPATCH framework (Nordlund et al. 2018) in a three-dimensional, Cartesian (shearing box) domain, with a set of static, nested patches. We employ a variation of the finite-difference STAGGER solver (Nordlund et al. 1994; Kritsuk et al. 2011) that, apart from the effects of transformation to a rotating coordinate system detailed below, solves the following set of partial differential equations:  \n@GLYPH<26> @ t = GLYPH<0>r GLYPH<1> ' GLYPH<26> u ' ; (2)  \n@GLYPH<26> u @ t = GLYPH<0>r GLYPH<1> ' GLYPH<26> uu + GLYPH<28> = ' GLYPH<0> r' p + p a ' GLYPH<0> GLYPH<26> r GLYPH<8> ; (3)  \n@ s @ t = GLYPH<0>r GLYPH<1> ' s u ' ; (4)  \non a staggered mesh, where GLYPH<26> is the gas density, u is the gas velocity, p is the gas pressure, p a is a stabilising artificial pressure term, GLYPH<28> = is a viscous stress tensor, s = GLYPH<26> log ' p GLYPH<157> GLYPH<26> GLYPH<13> 'GLYPH<157>' GLYPH<13> GLYPH<0> 1 ' is the entropy density per unit volume, GLYPH<13> is the ratio of specific heats and GLYPH<8> is the gravitational potential. The viscous stress tensor is taken to be  \nGLYPH<28> = = GLYPH<28> i ; j = GLYPH<23> 1 GLYPH<1> x ' c s + u ' GLYPH<26> GLYPH<18> @ u i @ x j + @ u j @ x i GLYPH<19> ; (5)  \nwhere GLYPH<1> x is the grid resolution, c s is the sound speed, and the coefficient GLYPH<23> 1 is a small fraction of unity. The artificial pressure is computed as:  \np a = GLYPH<23> 2 GLYPH<26> min ' 0 ; GLYPH<1> x 'r GLYPH<1> u '' 2 ; (6)  \nwhere GLYPH<23> 2 is a coefficient GLYPH<24> 1 (Kritsuk et al. 2011).  \nThis allows a discretisation that explicitly conserves  \nmass, momentum and entropy. The viscous stress tensor stabilizes the code in general, while the artificial pressure term ensures that shocks are marginally resolved.  \nThe equation of state (EOS) considered herein is that of an ideal gas with adiabatic index GLYPH<13> = 1 : 4 and molecular weight GLYPH<22> = 2 . Although not considered here, in a forthcoming paper we include radiative energy transfer and adopt a more realistic, table-based EOS (Popovas et al., in prep.).  \nWe use code units where length is measured in Earth radii, mass is measured in Earth masses, and velocity is measured in units of 1 km s GLYPH<0> 1 . Hence our time unit is R GLYPH<8> GLYPH<157> 1kms GLYPH<0> 1 = 6371 seconds, corresponding to roughly 5,000 time units per year.  \nPartly because the time step is determined locally in each patch, and partly because of the small number of cells per patch ( 33 3 in this case), the hydro solver operates near its optimal speed-updating on the order of 1.4 million cells per core per second-we performed these simulations using only a single, 20-core Intel Ivy Bridge node for each experiment.", '2.1 Grid set up': "We use nested sets of Cartesian patches centred on the embryo. With the embryo as the origin of our coordinate system, the x -axis is parallel to the cylindrical radial direction and points away from the central star, while the y - and z -axes are oriented along the azimuthal and vertical directions, respectively. The + y -axis lies along the direction of orbital rotation. The embryo is assumed to move in a circular orbit. Patches are arranged in a nested, hierarchical 'Rubik's Cube' arrangement,with 216 patches arranged into a 6 GLYPH<2> 6 GLYPH<2> 6 cube. The central 8 patches are then split recursively into 6 GLYPH<2> 6 GLYPH<2> 6 child patches, which have a three times finer resolution than the previous level. This is repeated for up to 7 levels of refinement (8 levels for the smallest embryo mass), while the coarsest set is duplicated 5 times in the y -direction. The size of the innermost, central, patch, exceeds the embryo radius R P by 10-30%. With 7 levels in the hierarchy, the outer scale in the vertical and radial directions is then GLYPH<24> 5000 R P , which is about 4 times a typical disk scale height for an Earth-mass planet, and GLYPH<24> 20x larger than the Hill radius. In the y -direction, the model extends to GLYPH<24> 25000 R P . Each patch contains 33 GLYPH<2> 33 GLYPH<2> 33 cells, leading to a maximum resolution in our simulations, near the embryo's surface, of GLYPH<25> 3% R P . This resolution is higher than has been used in previous works (e.g. D'Angelo & Bodenheimer 2013; Cimerman et al. 2017; Lambrechts & Lega 2017). Since the numerical viscosity is proportional to the cell size (cf. Eq. 5), the numerical viscosity in our simulations is also correspondingly smaller than in these previous studies 1 .", '2.2 Initial conditions': "We consider a local shearing flow in the reference frame that co-moves with the embryo at the local Keplerian frequency,  \nGLYPH<10> K . We ignore the slight curvature of the partial orbit but include the Coriolis force resulting from the orbital motion. The origin of the co-moving frame is, depending on the simulation run, placed at either 1 or 1.5 astronomical units (AU). The initial gas velocity is then given by:  \nu ' t = 0 ' = GLYPH<0> GLYPH<18> 3 2 + GLYPH<16> p GLYPH<19> GLYPH<10> K x ˆ y ; (7)  \nwhere GLYPH<16> p is a pressure trap parameter, defined below. In this study, we do not consider a systematic accretion flow towards the central star, nor the headwind that an embryo experiences when the gas is slightly sub-Keplerian. In analytically constructed accretion disks, with pressure decreasing systematically with radius, the headwind is typically on the order of one percent of the Keplerian velocity, but when an embryo is embedded in a local pressure bump there is no headwind. A systematic accretion flow would allow the smallest size particles to escape the pressure trap, and would thus tend to increase the fraction of large size particles, but since we are mainly concerned with measuring the accretion efficiency of different size particles, the possible evolution of the grain size distribution is not of direct importance for our study.  \nThe background vertical density stratification is set by the balance between the stellar gravity and the gradient of the gas pressure:  \nGLYPH<26> ' z ' = GLYPH<26> 0 exp GLYPH<0> GLYPH<10> 2 K z 2 + 2 GLYPH<16> p x 2 2 p 0 GLYPH<157> GLYPH<26> 0 ! ; (8)  \nwhere GLYPH<26> 0 and p 0 are the disk midplane density and pressure, respectively, and are estimated from the 2D radiative disk models of Bitsch et al. (2015a). The GLYPH<16> p term gives rise to a pressure gradient that balances the extra shear introduced in Eq. 7 above. A pressure bump produces a more self-contained region in the vicinity of the embryo's orbit, which retains the density wake produced by the embryo, preventing it from causing a forceful numerical bounce at the boundaries. This is also mitigated by the wave-killing boundaries that we use (see below), but a less forceful wake means a smaller region at the boundaries needs to be artificially damped.  \nWith typical dust opacities 2 on the order of 0.1 - 1 cm 2 g GLYPH<0> 1 , disks with surface densities similar to the Minimum Mass Solar Nebula (MMSN; Hayashi 1981) are optically thick in the vertical direction, since the MMSN nominally has 1700 g/cm 2 at 1 AU. If early hydrostatic atmospheres that surround planetary embryos are assumed to be nearly adiabatic one can predict their approximate vertical structure from knowledge of only the density, GLYPH<26> 0 , and temperature, T 0 , at the outer boundary. The pressure at the outer boundary of the atmosphere is assumed to match the local pressure of the surrounding disk. Therefore, for reasonable initial conditions, we approximate an adiabatic atmosphere profile outside-in to the embryo. From the EOS we have the  \nFigure 1. Radial temperature profile for a planetary embryo atmosphere using an ideal gas EOS ( GLYPH<13> = 1 : 4 , blue curve) and a realistic, tabular EOS (red curve).  \n<!-- image -->  \npressure at the outer boundary:  \nP 0 = GLYPH<26> 0 k B T 0 GLYPH<22> m u ; (9)  \nwhere k B is the Boltzmann constant and m u is the atomic mass unit. For each layer of the atmosphere, r i , we iterate over the set of equations:  \nGLYPH<26> i = P i GLYPH<22> m u k B T i ; (10)  \nGLYPH<1> ln P i = GLYPH<1> ln r i GM p 2 GLYPH<18> GLYPH<26> i + 1 P i + 1 r i + 1 + GLYPH<26> i P i r i GLYPH<19> ; (11)  \nGLYPH<1> ln T i = GLYPH<1> ln P i GLYPH<13> GLYPH<0> 1 GLYPH<13> ; (12)  \nwhere GLYPH<1> ln r i = 2 ' r i + 1 GLYPH<0> r i 'GLYPH<157>' r i + 1 + r i ' . The temperature and pressure for each layer are then simply T i = T i + 1 e GLYPH<1> ln T i and P i = P i + 1 e GLYPH<1> ln P i , respectively. We use 500 logarithmicallyspaced layers between half the embryo radius and the Hill radius. We then do a polynomial interpolation between these layers to map the density, pressure and temperature to each cell in the simulation box in order to construct the initial conditions in the vicinity of the embryo.  \nFor the purposes of this study, an ideal gas EOS ( GLYPH<13> = 1 : 4 ) is sufficient. As demonstrated by Fig. 1, the thermal profiles of the atmosphere surrounding an M = 0 : 95 M GLYPH<8> embryo given by an ideal gas EOS and a realistic, tabular EOS (Tomida et al. 2013, updated by Tomida & Hori 2016) are only significantly different close to the surface of the embryo. Gas flows in the vicinity of the Hill sphere, as well as embryo particle accretion rates, are insensitive to these differences.  \nWith constant opacities and an internal heating source (due to the accretion of solids) such an atmosphere is convective, and hence has a nearly adiabatic structure (Stevenson 1982). The atmosphere might have a radiative outer shell cf. the works by Piso et al. (2015) on giant planets and by Inaba & Ikoma (2003) for a terrestrial mass planet - but this shell would not substantially change the physical structure of the atmosphere, and ignoring it is a reasonable approximation when constructing the initial conditions. In forthcoming work, however, we include radiative energy transfer in our simulations, and there the initial conditions take this radiative shell into account.  \nThe effective, dynamical boundary between the disk and the embryo atmosphere occurs at distances from the planetary embryo somewhere in the range between the Hill radius  \nand the Bondi radius,  \nR B = GM p c 2 iso ; (13)  \nwhere c 2 iso is the isothermal sound speed. In fluid dynamics, the Bondi radius has no particular dynamical significance, other than being the radius where the pressure scale height of the atmosphere becomes comparable to the distance from the embryo, and the radius where trans-sonic bulk gas speeds would become unbound.  \nNote that the mass in such an adiabatic atmosphere is completely determined (apart from the equation of state) by the pressure and temperature at the outer boundary, and that if the external pressure is reduced the resulting mass is smaller. This implies that mass must be lost from such atmospheres when the disk evolves to lower densities (Nordlund 2011; Schlichting et al. 2015; Ginzburg et al. 2016; Rubanenko et al. 2017).", '2.2.1 Parameter space': "In this study, we consider 8 cases: 3 different embryo masses ( M p = { 0.95, 0.5, 0.096 } M GLYPH<8> ) at 2 distances from the central star ( a = { 1.0, 1.5 } AU), with different disk midplane densities ( GLYPH<26> 0 = { 1 10 GLYPH<0> 10 , 4 10 GLYPH<0> 11 } g cm GLYPH<0> 3 ) and temperatures ( T 0 = { 230, 300 } K), as well as with/without a pressure bump ( GLYPH<16> p = { 0.0, 0.05, 0.1 } ). The basic parameters for all cases are summarised in Table 1.  \nThe radii of the embryos are estimated using the massradius formula from Zeng et al. (2016):  \nGLYPH<18> R P R GLYPH<8> GLYPH<19> = ' 1 : 07 GLYPH<0> 0 : 21 CMF ' GLYPH<18> M M GLYPH<8> GLYPH<19> 1 GLYPH<157> 3 : 7 ; (14)  \nwhere CMF stands for core mass fraction. We take CMF = 0 : 325 , although the exact value is not important; taking CMF = 0 : 0 only changes the radius of the embryo on the order of 5%.", '2.3 External forces (accelerations)': "External forces, arising from the transformation to a rotating coordinate system, and from linearising the force of gravity due to the central star, are added to the hydrodynamics as source terms:  \n@ u @ t = GLYPH<0> GM p ' x ˆ x + y ˆ y + z ˆ z ' r 3 + 3 GLYPH<10> 2 K x ˆ x GLYPH<0> GLYPH<10> 2 K z ˆ z GLYPH<0> 2 GLYPH<10> K ˆ z GLYPH<2> u ; (15)  \nwhere r is the distance from the embryo. The first term in Eq. 15 is the two-body force from the embryo. The second term is the tidal force and has two contributions: the centrifugal force and the differential change of the stellar gravity with radius. Vertically, there is a contribution from the projection of the force of gravity onto the z-axis (the third term). The last term is the classical Coriolis force.  \nTo avoid the singularity of the two-body force, rather than using a 'softening' of the potential (e.g. Ormel et al. 2015a; Lambrechts & Lega 2017; Xu et al. 2017), we truncate it via r = max ' r ; 0 : 1 R P ' . Truncation is better than softening, since the gas outside the radius of the embryo feels the correct force, and the spherical boundary conditions (described below) negates the effect of truncation for the innermost  \nTable 1. Basic parameters for the different simulation runs. See the text for more details.  \ncells, well within the radius of the spherical boundary condition/embryo.", '2.4 Boundary conditions': "Using test simulations, we have ensured that our boundary conditions remain stable over long periods of integration time and under rather extreme conditions (e.g. highly supersonic flows, low densities). The 'azimuthal' y -boundaries are set to be periodic, whereas the vertical z - and 'radial' x -boundaries are set as follows:  \n- GLYPH<15> Zero-gradient conditions for the entropy per unit mass.\n- GLYPH<15> The x - and z -velocity components are set symmetric with respect to the boundaries, while the y -component is extrapolated in the x -direction and z directions.\n- GLYPH<15> The density boundary conditions in the x -direction is constructed from assuming hydrostatic balance:  \nGLYPH<26> x = GLYPH<26> x i exp GLYPH<0> GLYPH<10> 2 K GLYPH<16> p T 0 ' x 2 GLYPH<0> x 2 i ' ! ; (16)  \nwhere GLYPH<26> x is the density in the boundary at location x and GLYPH<26> x i is the density at location x i , which corresponds to the last cell inside a patch's active domain.  \n- GLYPH<15> If a pressure trap is not present we assume a vanishing derivative of the density in the x -direction.\n- GLYPH<15> In the vertical direction, we apply a similar condition:  \nGLYPH<26> z = GLYPH<26> z 0 exp GLYPH<0> GLYPH<10> 2 K 2 T 0 ' z 2 GLYPH<0> z 2 0 ' ! ; (17)  \nwhere GLYPH<26> z is the density in the boundary at location z and GLYPH<26> z 0 is the density at location z 0 , which again corresponds to the last cell inside the active domain.  \nWe additionally apply the 'wave-killing' prescription of de Val-Borro et al. (2006) to reduce reflections at the x -boundaries. It is only applied over a specific range, given by x damp (column 9 in Table 1).", '2.4.1 Spherical boundaries': "In these simulations we use 'jagged', no-slip boundaries to approximate the embryo surface. Fig. 2 shows a 2D representation of these boundary conditions on a staggered mesh. Density and entropy are allowed to evolve freely everywhere; only the mass flux is constrained. We deem a set of cell  \nFigure 2. Treatment of spherical boundary conditions on a staggered mesh. Blue and red squares represent horizontal and vertical mass fluxes outside the boundary, which evolve freely, whereas turquoise and magenta squares represent mass fluxes inside the boundary, which are set to zero. Green squares represent density points that are subject to change, while yellow squares show density points that remain constant throughout the simulation.  \n<!-- image -->  \nfaces that approximate the spherical surface to belong to the boundary. If the cell face is inside the spherical boundary, the mass flux is set to zero, otherwise it is allowed to evolve freely. Hence, density cannot change inside the boundary and cells that have faces belonging to the boundary react exactly as they should: inflows towards the boundary produce a positive @GLYPH<26> GLYPH<157> @ t at the boundary. The corresponding increase in pressure causes an increase in the outward directed pressure gradient at the last 'live' mass flux point, as expected at a solid boundary. This permits the system to find an equilibrium and react to perturbations from it, and since there are never any inflows and outflows through the boundary, mass is conserved.", '2.5 Relaxation of initial conditions': "We allow the system to relax for 30 000 time units (which corresponds to about 6 years), in order to ensure that dynamically consistent conditions are attained. Indeed, as discussed below, tests demonstrate that it takes on the order of only one orbit for the initial conditions to relax.  \nThe density stratification well inside the Hill sphere is almost hydrostatic, but with deviations and fluctuations on short time scales, and it is thus not necessary to include the detailed dynamics of this region while larger structures are relaxing. Therefore, we begin our experiments by fixing the hydrostatic structure there, and set the spherical boundary to a radius of 27 embryo radii. We then 'release' layers (add an extra level of refinement, reduce the radius of the  \nspherical boundary, turn on the dynamics outside of the new spherical boundary) at successively smaller radii (9, 3 and finally, 1 embryo radii). The 'releases' are typically done at times = 10 000, 20 000 and 26 000 in code units. Successive releases are increasingly more expensive due to the smaller cell size and larger sound speeds, but need less time in code units to adjust to the surrounding environment (for essentially the same reasons, i.e., smaller scales and larger sound speeds). This procedure saves significant amounts of computing time during the initial relaxation phase.  \nAfter the last release, the simulations are allowed to relax for an additional 4000 time units before they are deemed ready for the next phase. Figures 3 and 4 show largescale and zoomed-in density slices at the midplane for run m095t00 at t = 0 and after a relaxed state is reached at t = 30 000. Panels on the right show the relaxed state with well established horseshoe orbits, density wakes, and slight density perturbations ahead/behind of the embryo at the same orbital distance (discussed in more detailed below).  \nFig. 4 is a zoom-in of the same run with the Bondi radius more clearly visible (largest red circle); the Rubik's cube nesting of grids is also visible. The dashed lines indicate the successive releases of the radius of the embryo (green - 27, yellow - 9, light red - 3 embryo radii). The tiny black circle corresponds to one embryo radii. The right panel shows the simulation when fully relaxed and, via streamlines, demonstrates that the disk gas penetrates deeply into R B , down to GLYPH<25> 20 R P . At even smaller radii the atmosphere is nearly static, with motions that have a weak rotational tendency (see below).  \nFig. 5 shows the physical structure of the primordial atmosphere at both time = 0 and 30 000. The cuts are done along the x -axis in the midplane and along the z -axis at x = y = 0 . As can be seen, the initial state is very close to the final state of the atmosphere, and there are only negligible variations for different cuts. Although the density and temperature profiles have changed slightly, the changes are minimal and the use of these initial conditions thus saves a lot of relaxation time (compared to, e.g., a sudden insertion of the planet's gravitational potential, which would require waiting tens of orbits with the maximum level of refinement for the dynamics to relax). Although the indnermost parts of the atmosphere still retain some weak random motions, by time = 30 000, the run has nearly reached a steadystate, when injection of particles can commence. Note that, as has been pointed out by Ormel et al. (2015b), a completely steady, static state will likely never be reached in the vicinity of the embryo.", '3 PARTICLES': "After relaxation of the fluid motions and density distribution for 30,000 time units (about 6 years) we add particles for the purpose of tracing the motions of selected sub-samples and reveal, for example, the typical paths of different physical size particles and to be able to measure the rate of accretion onto the planetary embryos. In these simulations, we consider macro-particles (sometimes called super-particles; Rein et al. 2010; Ebisuzaki & Imaeda 2017), each of which represents a swarm of real particles, with a given size and mass. We follow the particle motions, taking into account ex-  \nternal forces and gas drag, until they are accreted onto the planet embryo, or exit the domain through the x -boundaries. We consider particles as 'accreted' when their radial distance to the centre of the planet is smaller than the radius of the embryo.  \nIn order to have statistical coverage everywhere, the initial spatial distribution of macro-particles is chosen to be proportional to the local gas density. The actual dust-to-gas ratio is likely characterized by a strong settling towards the midplane, while the dust-to-gas surface density ratio should be essentially constant in the upstream flow entering the Hill sphere. However, rather than having to make assumptions about the settling, we instead analyse sub-populations of our initial distribution, tagging and following for example only the particles that initially reside within a given distance from the midplane.  \nWe use macro-particles with corresponding real particle sizes normally ranging from 10 GLYPH<22> m to 1 cm, with a constant number of macro-particles per logarithmic size bin. However, note that we never integrate over size. The constant number per bin was chosen to provide good statistical coverage for all sizes when investigating the transport of particles in narrow size bins.  \nWe include 1 cm particles for completeness, even though such large particles are typically found only in CB chondrites (Friedrich et al. 2015), partly because there is no guarantee that the size distribution of chondrules in chondrites spans the entire size range that was available at the time of pebble accretion onto planet embryos. We also made complementary experiments with even larger size particles (up to 1 m), to explore saturation of accretion rates with size.  \nThe shape of the relative size distribution will be modified near the embryo, because large particles accrete more rapidly onto the embryo, while being replenished at the same rate as smaller size particles. The resulting change of the relative size distribution is a result of this study, and does not need to be anticipated in the initial conditions.  \nNote that, since we do not consider particle growth and destruction in this study, nor the back reaction of the particles on the gas, we are free to renormalise the the actual size distribution of particles a posteriori by changing the interpretation of the statistical weight factors of the macroparticles. We can thus match any given initial actual size distribution, and derive what the final distribution would be. We also note that since the span in space and time of our models is very limited (relative to disk extents and lifetimes), ignoring coagulation and shattering is an excellent approximation; such effects may instead be taken into account by modifying a posteriori as per the remarks above - the initial particle distribution.", '3.1 Equations of motion': "The acceleration of macro-particles,  \nd u p dt = a G + a d ; (18)  \nis governed by a G , the differential forces originating from the planet, the star and the Coriolis force, which, in analogy with Eq. 15, are  \na G = GLYPH<0> GM p ' x ˆ x + y ˆ y + z ˆ z ' r 3 + 3 GLYPH<10> 2 K x ˆ x GLYPH<0> GLYPH<10> 2 K z ˆ z GLYPH<0> 2 GLYPH<10> K ˆ z GLYPH<2> u ; (19)  \nFigure 4. Relaxation of initial conditions for the m095t00 run. Intermediate scale density slices at the midplane. The large red circle is the Bondi radius. The solid green circle in panel (a) shows the current 'active' radius, while subsequent dashed circles show where future 'releases' will place a temporary embryo surface. The tiny black circle shows the actual radius of the embryo, which is reached after the final release.  \n<!-- image -->  \n(a) Initial conditions at time = 0.  \n<!-- image -->  \n(b) Relaxed simulation at time = 30 000  \nFigure 3. Relaxation of initial conditions for the m095t00 run. The panels show density slices and velocity streamlines at the midplane, with almost the entire extent of the simulation box in the x -direction and 1/5 of the way in the y -direction. Horseshoe orbits are clearly visible in panel (b). Panel (b) is plotted with a highly enhanced contrast (the colour scale spans only 20% in density) to show the wakes and density increase at the orbital radius in the y -direction. The white and red circles (barely visible in these figures) denote the Hill radius and the isothermal Bondi radius of the embryo, respectively. The embryo itself is smaller than a pixel.  \n<!-- image -->  \n<!-- image -->  \n(a) Initial conditions at time = 0.  \n(b) Relaxed simulation at time = 30 000  \nand the acceleration from the gas drag,  \na d = u gas GLYPH<0> u part t s ; (20)  \nwhich is a product of the relative velocity of the particle, u part , with respect to that of the gas, u gas , and the inverse of the stopping time, t s , which, in the Epstein regime is:  \nt s = GLYPH<26> GLYPH<15> s GLYPH<26> v th ; (21)  \nwhere GLYPH<26> GLYPH<15> is the solid density (taken to be 3 g cm GLYPH<0> 3 ) and v th is the mean thermal velocity of the gas. We further assume that particles are spherical.  \nFor simplicity and speed, we parametrise the stopping  \ntime as:  \nGLYPH<18> t s 1 yr GLYPH<19> = GLYPH<16> s 1cm GLYPH<17> GLYPH<18> c d GLYPH<26> GLYPH<19> ; (22)  \nwhere c d = 10 GLYPH<0> 12 g cm GLYPH<0> 3 is a normalisation constant, defined such that in the sub-sonic Epstein regime, the stopping time is 1 year for a 1 cm particle in a gas density of 10 GLYPH<0> 12 g cm GLYPH<0> 3 . This estimates the stopping time in the disk and outer parts of the Hill sphere sufficiently well, while ignoring the additional decrease of the stopping time due to the higher temperature near the surface of the embryo. Particles that reach such depths of the potential well are in any case doomed to accrete, so nothing is lost by ignoring the  \nFigure 5. Physical structure of the primordial atmosphere for run m095t00 at t = 0 (initial conditions; IC) and t = 30 000 = 6.1 years (fully relaxed). The horizontal and vertical cuts overlay each other nearly perfectly.  \n<!-- image -->  \ntemperature dependence, and doing so reduces the cost to compute the stopping time.  \nAlthough this approximation may overestimate the stopping time slightly, the precise value is not critical, since there is in any case an uncertain factor that depends on the shape and porosity of the particles, and a different value can be compensated for by a corresponding shift of the assumed particle distribution.  \nIn a similar vein, there is no need to expand the parameter space with different assumed disk densities, since that too corresponds to a rescaling of the particle stopping time, and a corresponding rescaling of the disk and adiabatic atmospheric structures with a constant factor in density.  \nBecause of its dependence on gas density and particle size, the stopping time can vary over many orders of magnitude. Given the time step, GLYPH<1> t , that is necessary based on the velocity and acceleration of gravity, a G , one can thus have both GLYPH<1> t GLYPH<28> t s and GLYPH<1> t GLYPH<29> t s . The latter case means that the differential equation governing the particle path becomes stiff .  \nTo handle this, particle positions and velocities are updated with a modified kick-drift-kick leapfrog scheme (e.g. Dehnen & Read 2011), reformulated to give the correct relative speed also in the asymptotic limit when the time step is much larger than the stopping time. Instead of sub-cycling, which is computationally expensive, we use the method that is commonly referred to as'forward time differencing', which ensures that the solution has the correct asymptotic behaviour, both when the stopping time is very long and very short, and that the behaviour is also nearly the correct one in intermediate cases.  \nWith forward time differencing, the velocity update, assuming for now only an acceleration due to the drag force, a d , may be written:  \nv ' t + GLYPH<1> t ' = v ' t ' + GLYPH<1> v = v ' t ' + GLYPH<1> t a d ' t + GLYPH<1> t ' ; (23)  \nwhere the acceleration is time-centred at t + GLYPH<1> t . The resulting update relation is:  \nGLYPH<1> v = GLYPH<1> t a d ' t + GLYPH<1> t ' = GLYPH<0>' GLYPH<1> t GLYPH<157> t s ' ' v + GLYPH<1> v GLYPH<0> v g ' t + GLYPH<1> t '' ; (24)  \nwith solution:  \nGLYPH<1> v = GLYPH<0> GLYPH<1> t GLYPH<157> t s ' 1 + GLYPH<1> t GLYPH<157> ts ' ' v GLYPH<0> v g ' t + GLYPH<1> t '' : (25)  \nIn the limit GLYPH<1> t GLYPH<157> t s GLYPH<29> 1 , the velocity ends up very close to the gas velocity (as it should), but without the need to take very short time steps, resulting in a much more efficient time integration. In the limit GLYPH<1> t GLYPH<157> t s GLYPH<28> 1 , one recovers the normal kick-drift-kick behaviour.  \nWe thus advance the velocity in the first (kick) step with:  \nv 1 = v 0 + GLYPH<1> t ' 1 + GLYPH<1> t GLYPH<157> t s ' ' a G ' r 0 ; t 0 ' + a d ' r 0 + v 0 GLYPH<1> t ; t 0 + GLYPH<1> t '' : (26)  \nThe particle position and time are then updated (drift step) with:  \nr 2 = r 0 + v 1 2 GLYPH<1> t ; (27)  \nand  \nt 2 = t 0 + 2 GLYPH<1> t : (28)  \nFinally, the third step (kick) advances the velocity using the acceleration at the end of the full time step, with  \nv 2 = v 1 + GLYPH<1> t ' 1 + GLYPH<1> t GLYPH<157> t s ' » a G ' r 2 ; t 2 ' + a d ' t 2 '… : (29)  \nWe have tested this formulation on particle drift in disks with known solutions (Weidenschilling 1977), recovering behaviour sufficiently accurate for our purposes for arbitrary values of Stokes numbers St = GLYPH<10> t s .  \nThe particle update procedure is computationally efficient, with the average update cost per cell on an Intel Ivy Bridge core increasing from 0.65 GLYPH<22> s without particles to 0.72 GLYPH<22> s with one particle per cell, indicating an average update cost of about 70 nanoseconds per particle. Most of that time is spend looking up properties needed to calculate the particle-gas drag force.", '3.2 Injection of pebbles': "Once the simulations are considered to be fully relaxed, pebbles are injected into the system. The probability that a particular cell in a particular patch will spawn a macroparticle is proportional to N ppc , the number of particles per cell (which need not be an integer). We set N ppc to 0.05 in most runs, but increase it to 1.0 in a few runs to improve the statistics when studying the distribution of particles in small volumes.  \nThe mass carried by a macro-particle is encoded in a 'weight factor', w = GLYPH<26> GLYPH<1> V GLYPH<157> N ppc , where GLYPH<1> V is the cell volume. The weight factor is thus proportional to the mass of gas in the cell where the macro-particle originates divided by N ppc . An advantage of this approach is that the sum over all particles (or a subset defined by particle size or initial particle location) results in the total mass of those particles (to be renormalized by a preferred initial dust-to-gas ratio).  \nHowever, because the mass contained in any given cell can vary substantially, particles may need to be split to avoid the situation where the particle mass in a patch is represented by too few particles. For purposes of efficiency, the need for splitting is only checked when a particle passes from one patch to another. If an inbound particle carries more weight than, e.g., five times the initial average in the  \n(a) Midplane.  \n<!-- image -->  \n(b)  \nFig. 6 shows the details of the gas flow in the m095t00 run for the M = 0.95 M GLYPH<8> planetary embryo. For clarity, we zoom-in to the Hill radius (234 embryo radii, white circle). The red circle denotes the nominal Bondi radius, R B , which is GLYPH<25> 50 R P . The embryo is marked by a black circle, but is too small to be visible in these figures. The first four panels show the velocity magnitude in the disk midplane ( z =0), and at z = 10, 25, and 75 R P . As expected, the perturbation of the gas due to the embryo decreases with increasing height above the midplane, as indicated by the streamlines 3 . The velocity magnitude in the horizontal plane 75 planet radii above the midplane is nearly uniform along the y -direction, indicating that, at that height, there are no substantial perturbations where the embryo is located, consistent with Ormel et al. (2015b). We also observe that the width of the horseshoe orbits decrease with increasing height above the midplane,  \n<!-- image -->  \n(d)  \nz  \n= 75  \nR  \nP  \n.  \n<!-- image -->  \nz  \n= 10  \nR  \nP  \n.  \n(c)  \nz  \n= 25  \n(e) Vertical component of the vorticity in the midplane.  \n<!-- image -->  \n<!-- image -->  \nR  \nP  \n.  \n(f) xz slice of u z at y = 0.  \n<!-- image -->  \nFigure 6. Gas flows around the M = 0.95 M GLYPH<8> embryo in the m095t00 run. Panels (a-d) show the velocity magnitude at different heights above the midplane. The scale of the colour bar is deliberately held constant. Panel (e) shows the vertical component of the vorticity in the midplane and panel (f) shows u z in a xz slice at y = 0. Velocity streamlines are also shown.  \npatch, it is split into a number of 'child particles'. The child particles continue to carry, also after any additional splits, a memory of their initial identity, which is needed for identifying sub-populations based on initial physical location.  \nWhen a macro-particle is created, it is assigned a random position inside a cell and a randomly chosen size: s = 10 GLYPH<0> 3 R cm ; (30)  \nwhere R GLYPH<18> [0,1] is a random number. The spawned macroparticles are initially given the gas velocity of the cell. We spawn GLYPH<25> 2.6 million macro-particles in most runs, with 53 million in a few runs with N ppc = 1 . We save a complete set of particle data when a snapshot of the gas dynamics is taken, but, for visualization purposes, we also separately save trajectory information for a smaller subset of macro-particles ( GLYPH<24> 15 000) at every particle time step. These 'tracer particles' are used to accurately track the trajectories of accreting particles and correspondingly to show the escape routes of particles that enter the Hill sphere and then leave again.", '4 RESULTS: GAS FLOWS': "We begin by investigating the gas flow patterns around the different mass cores. In the last several years, a number of authors have studied the hydrodynamics in the vicinity of embedded planetary embryos (e.g. D'Angelo & Bodenheimer 2013; Ormel et al. 2015b; Fung et al. 2015; Masset & Ben'ıtezLlambay 2016; Cimerman et al. 2017; Lambrechts & Lega  \n2017; Xu et al. 2017, etc.) across a broad range of different conditions (e.g. different core masses, orbital distances from the parent star, disk structures), physical processes (isothermal/adiabatic/realistic EOS, with/without radiative energy transfer, magnetic fields, etc.). As illustrated below, our simulations are consistent with these previous studies, while offering both a higher resolution and a larger domain around the embryo.  \n<!-- image -->  \n<!-- image -->  \n(a) Midplane.  \n<!-- image -->  \n(b)  \n<!-- image -->  \n(d)  \nz  \n= 75  \nR  \nP  \n.  \nz  \n= 10  \nR  \nP  \n.  \n(c)  \nz  \n= 25  \nR  \n(e) Vertical component of the vorticity in the midplane.  \n<!-- image -->  \nP  \n.  \n(f) xz slice of u z at y = 0.  \n<!-- image -->  \nFigure 7. Same as Fig. 6, but for the M=0.096 M GLYPH<8> embryo in the m01t00 run.  \n<!-- image -->  \n(a)  \nm01t00  \nrun.  \n<!-- image -->  \nFigure 8. Velocity magnitude in the midplane in the vicinity of the embryos as a function of embryo mass.  \n<!-- image -->  \nwhich agrees with the results of Fung et al. (2015), but contradicts the isothermal results of Masset & Ben'ıtez-Llambay (2016). This might be attributed to the fact that we use a non-isothermal EOS, but further analysis of the horseshoe orbits is beyond the scope of this work. Panel (e) of Fig. 6 shows the vertical component of the vorticity measured relative to the rotating coordinate system. A clear vorticity pattern can be seen, with separators splitting the flow into regions that pass the embryo on the left and right sides. The flows passing on the 'wrong' side contribute directly to positive (prograde, counter-clockwise) vorticity, which adds to the rotation of the flow.  \nPanel (f) of Fig. 6 shows the vertical velocity, v z , in  \nan xz plane at y = 0. A vertical gas flow pattern can be distinguished: gas flows towards the embryo from the poles and moves away at the midplane. Zonal flow patterns are visible as closed streamline loops. Vertical velocities inside R B are of the order of a few cm s GLYPH<0> 1 , which agrees well with previous works(e.g. Ormel et al. 2015b; Lambrechts & Lega 2017).  \nQualitatively, the M = 0.5 M GLYPH<8> case is very similar, thus the figure panels for this mass are presented in Appendix A; Fig. A1.  \nThe gas flows around the M = 0.096 M GLYPH<8> embryo ( m01t00 run) behave somewhat differently, as illustrated in Fig. 7. Because of the lower embryo mass, the perturbations  \naway from pure shear flow are smaller, and horseshoe orbits therefore reach further in towards the embryo surface. The strong features in vorticity are correspondingly smaller. The vertical flows of the gas, meanwhile, approach the embryo over a broader extent radially, relative to the size of the Hill sphere. In our simulations, the gas flow patterns in the innermost region change with time, alternating between dominating gas inflows from the leading/outer horseshoe flow and the trailing/inner one. The locations of the stagnation points also vary with time. However, it is possible that part of this behaviour is related to the recursive refinement of the grid. We performed exploratory runs with the same mesh refinement, but without the planet embryo, and found root-meansquare velocities in the innermost region on the order of 10 m s GLYPH<0> 1 . This is similar to the case with the low-mass embryo present, and hence the small scale, low amplitude motions are unreliable. These motions are not caused by noise generated directly by the mesh refinement, but rather by slight differences in the cancellation of the force of gravity and pressure gradients, where small differences due to numerical resolution result in imperfect hydrostatic balance.", '4.1 Gas dynamics close to the planetary embryos': 'One reason why this work stands out relative to previous studies is that we are able to resolve the detailed gas dynamics (and hence accurate particle paths) in the close vicinity of the planetary embryos. Fig. 8 shows the gas velocity magnitude in the vicinity of R B for the M = 0.096 M GLYPH<8> (panel (a)), M = 0.5 M GLYPH<8> (panel (b)) and M = 0.95 M GLYPH<8> (panel (c)) embryo masses. The M = 0.096 M GLYPH<8> core is not massive enough to form a hydrostatic atmosphere - the horseshoe flows constantly flush the atmosphere and only a tiny part of it ( GLYPH<25> 2 R P from the embryo) is nearly static. With increasing embryo mass, more gas becomes bound and the velocity magnitude decreases. Nevertheless, there are always slow gas movements inside R B due to both gas from the disk that reaches deeper interior layers and from small scale turbulence. Therefore, the primordial atmospheres are never fully isolated and always interact with the surrounding disk.', '5 RESULTS: PARTICLES': 'An obvious advantage of our large simulation box is the large particle reservoir. Smaller simulation boxes, especially with non-periodic x -boundaries, often experience a depletion of macro-particles. This makes the accurate evaluation of accretion rates more difficult, and additional injections of particles may be needed (e.g. Xu et al. 2017).  \nIn Fig. 9, we show time-series of particle mass density from run m095t10 across a range of integrated radial layers (shells) from the embryo surface to twice the Hill radius. Black denotes the first snapshot after the injection of the pebbles (after 40 time units) while the other lines are plotted at 4 000 time unit intervals (which corresponds to GLYPH<25> 0.8 of an orbital period). We also plot the integrated density of the gas divided by 100 (the canonical gas-to-dust ratio), in cyan for the same radial bins. The bottom panel shows the number of macro-particles in each radial bin. Here, the initial dips on the black curves stem from changes in the level of refinement. As time goes on, the number of macro-particles  \nFigure 9. Mass distribution of solids in radial shells close to the embryo at different times for the m095t10 run. The vertical dashed and dotted lines show the Bondi and Hill radii, respectively. The cyan line denotes the gas density in each shell, divided by a canonical dust-to-gas ratio of 100. The lower panel shows the number of macro-particles as a function of radius. The bin width dr is 4 : 66 R P . As initially local macro-particles with low weights are replaced by macro-particles with larger weights from larger distances, the statistical noise increases.  \n<!-- image -->  \nFigure 10. Vertical distribution of pebble sizes in the m095t00 run 1.6 orbits after the injection of macro-particles.  \n<!-- image -->  \ndecreases due to accretion and transport, but new particles are continuously coming from further out. These new particles originate in patches with lower numerical resolution (i.e. larger cell sizes), and thus carry larger representative weights. The increasing noise level, meanwhile, is due to the increasing contribution from these higher statistical weight macro-particles.  \nEven after 2.4 orbits, Fig. 9 shows no indication of depletion of solids - the mass influx through a given shell is sufficient to replace any accreted or lost particles. This implies that our initial injection of particles is sufficiently numerous to allow accurate measurements of accretion rates without further injections. Similar trends are seen for the other embryo masses; as the solid density is proportional to the gas density, we choose to not show separate figures for them.  \n(a) Particle distribution over a 2R H vertical extent above/below the midplane.  \n<!-- image -->  \n(b) Same as (a), but with a 2R B vertical extentFigure 11. Gas density and spatial distribution of particles close to the midplane in the m095t00 run. The particles here are split into 3 size bins: 0 : 001 < s GLYPH<20> 0 : 01 cm (yellow), 0 : 01 < s GLYPH<20> 0 : 1 cm (red), and 0 : 1 < s GLYPH<20> 1 cm (green). The black circle denotes the Hill sphere, while the red circle denotes the Bondi radius.  \n<!-- image -->', '5.1 The distribution of particles': "We initialize particles uniformly in the dust-to-gas ratio, and we therefore assume the initial vertical stratification of particles is the same as that of the gas. Particles then settle towards the midplane due to the vertical force of gravity, with settling time scales of the order of 1 GLYPH<157>' t s GLYPH<10> 2 K ' . The settling time scale for small and intermediate size particles particles is long relative to the duration of our experiments, while it is comparable to our run times for the largest particles. Lacking proper turbulence in the disk, the settling of particles is anyway not very realistic-turbulence stirs up the particles and controls their scale height. Additionally, zonal gas flows (e.g. panel (f) of Fig. 6) close to the Hill sphere are capable of stirring up well-coupled particles.  \nTherefore, instead of working with ill-defined and unknown equilibrium populations, we choose instead welldefined, localized sub-populations with initially constant dust-to-gas ratios for our measurements of accretion rates. Tests have shown that it takes only GLYPH<24> 0.5 orbits after the injection of particles for the measured accretion rates to stabilize. The hydrodynamic flow patterns are established quickly as well, taking approximately the same amount of time. This shows that our initial relaxation of the flow patterns, which was allowed to continue for GLYPH<24> 6 orbital periods, was certainly sufficient and allows for accurate measurements of accretion rates.  \nFig. 10 shows the vertical distribution of pebbles in the m095t00 run, 1.6 orbits after the injection of macro-particles. The smallest particles ( s GLYPH<20> 0 : 01 cm; dark colours) continue to show essentially the vertical distribution of the gas, while larger particles ( s > 0 : 1 cm; green to yellow) have begun to settle. This figure also nicely illustrates how small a pertur-  \nbation the Hill sphere is relative to the size of our computational domain.  \nA zoom-in to R B (panel (c)) reveals that there is a concentration of small particles ( s < 0 : 1 cm) inside the Bondi sphere. Larger particles do not linger there for long since they are effectively accreted. The figure shows particle paths based on high-cadence positions of tracer-particles, sub-divided into 3 size bins. The green colour represents 0 : 001 GLYPH<20> s < 0 : 01 cm, blue 0 : 01 GLYPH<20> s < 0 : 1 cm and black s 0 : 1 GLYPH<20> s < 1 cm size particles. The smallest particles spend around 0.1 GLYPH<10> GLYPH<0> 1 inside the Bondi sphere before they reach  \nFig. 11 shows the spatial distribution of particles and their sizes at the midplane, two orbits after the injection of particles. This time, the particles are split into 3 size bins: s GLYPH<20> 0 : 01 cm (yellow), 0 : 01 GLYPH<20> s < 0 : 1 (red) and s > 0 : 1 cm (green). The vertical extent of particles we display is limited to 1 R H in panel (a) and 2 R B in panels (b,c), above and below the midplane. One instantly notices the two leading and trailing arms, stretching away from the embryo. The large concentration of particles there indicates that these particles are moving away from the embryo very slowly. The wakes are located just inside the separatrix sheets of the horseshoe orbit, with the smallest particles mostly concentrating in the outermost regions and the larger ones lying closer to the orbit of the embryo. These macro-particles carry a relatively small statistical weight, because they originate at patches with high resolution and small mass per cell, so the features are somewhat artificially enhanced in the Figure, but it therefore also nicely illustrates the preferred 'leaving' paths of the particles. A large fraction of these particles enter but then later leave the Hill sphere, and the features are therefore concentrated in the midplane.  \n<!-- image -->  \nFigure 12. Trajectories of accreted particles. The axes are given in units of embryo radii. Green, blue and black paths are for particle sizes 0 : 001 GLYPH<20> s < 0 : 01 cm, 0 : 01 GLYPH<20> s < 0 : 1 cm, and 0 : 1 GLYPH<20> s < 1 cm, respectively.  \n<!-- image -->  \nthe embryo. In contrast, particles with s GLYPH<21> 0 : 01 cm, are accreted on average 10 times faster, after they reach the Bondi sphere. Alternatively, some of the smallest particles, if they linger far enough from the embryo, are later carried out by the horseshoe gas flows.  \nIn these runs we do not consider convection driven by the accretion heating (but see below), particle growth and evaporation. Sufficiently strong convection could possibly prevent the smallest particles from reaching the core. Also, as the gas temperature close to the embryo ( GLYPH<25> 4500 K near the M = 0.95 M GLYPH<8> embryo) greatly exceeds the melting point of any solid, evaporation of solids and the resulting enrichment of the atmosphere with heavy elements (e.g. Alibert 2017; Brouwers et al. 2017), is bound to play an important role in models with a realistic treatment of evaporation. As shown by Fig. 1, a more realistic equation of state leads to lower temperatures at the base of the atmospheres, but still large enough to cause evaporation.  \nFig. 13 similarly shows the spatial distribution of pebbles close to the midplane (inside 2 R B ) in the m01t00 run. It shows a very similar picture - escaping particles are trapped in the leading/inner and trailing/outer horseshoe orbit parts. The zoom-in of the Figure (panel b) reveals an interesting, but expected difference - there is no accumulation of particles inside the Bondi sphere. There are two factors that contribute to this: the gas flows may be strong enough to readily carry small and intermediate size particles away, and larger particles are rapidly accreted. Panel (c) in Fig. 12 reveals the answer: the stopping time, even of the smallest size particles we consider in our runs, is short enough that the gravitation of the embryo is not damped by the gas drag enough to prohibit particles from being accreted. The smallest particles in our simulation take merely a 0.04 GLYPH<10> GLYPH<0> 1 to be accreted. And the number of these smallest particles crossing the Bondi sphere is much larger, compared to the M = 0.95 M GLYPH<8> embryo. Here we remind the reader that we do not trace all the particles with this high cadence, therefore only a fraction of all the accreted particles are shown in Fig. 12 (we also, for greater clarity, limit the number of traces shown in these figures to 150).  \nFig. 12 shows evidence for a variation of hydrodynam-  \nical deflection (e.g. Sellentin et al. 2013) - higher mass embryos are accompanied by much denser primordial atmospheres, in which the stopping times become smaller. Whereas the majority of the particles in the distribution we consider here have nearly ballistic orbits when being accreted by the lowest mass embryo (Fig. 12, panel (c)), they become significantly affected by the higher gas density envelopes in panels (b) and (a). In the latter, the highest-mass embryo case, the motion of the gas almost completely controls the motion of particles with size s < 0 : 01 cm.", '5.2 Accretion of pebbles': "When pebbles approach an embryo they are subject to an increasing two-body force and start to veer towards the embryo. Although from Fig. 12 it might seem that they come from all over the place, the figure greatly under-represents the particles (150 traces selected out of 15 000, again randomly selected from GLYPH<25> 2.6 million macro-particles), and the scales are concentrated to the closest 10-20 planetary radii around the embryos. The full ensemble of trace-particles in fact show distinct inflows and outflows through the Hill sphere.  \nFig. 14 plots the mass accretion rates across the Hill sphere, Bondi sphere and the embryo surface over the course of the m095t10 run, displayed as a set of Hammer projections. The macro-particles are sub-divided into 3 representative size-bins ( 0 : 0 GLYPH<20> s GLYPH<20> 0 : 01 cm, 0 : 01 < s GLYPH<20> 0 : 1 cm and 0 : 1 < s GLYPH<20> 1 : 0 cm). The figures show that the bulk of the net accretion comes from near the midplane, and predominantly from the large size particles. As they do not feel much gas drag, some of the largest particles can and do cross the Hill sphere from arbitrary directions. The dominating inflow regions are associated with the inner/trailing and outer/leading horseshoe orbits, and corresponding outflows in the outer/trailing and inner/leading horseshoe orbits. The accretion rates at distances R B are much more isotropic, without clearly preferred trajectories. The modest outflows are mainly close to the midplane region. Not all of the mass accreted through the Bondi sphere reaches the embryo surface during the time scales we consider here. The  \n(a) Spatial distribution of particles over a 2 R B vertical extent above/below the midplane. The black circle denotes the Hill sphere.  \n<!-- image -->  \n(b) A zoom-in of panel (a) to the Bondi sphere. The red circle denotes the Bondi sphere.  \n<!-- image -->  \nFigure 13. Gas density and spatial distribution of particles close to the midplane in the m01t00 run. The pebbles here are split into 3 size bins: yellow s GLYPH<20> 0 : 01 cm, red 0 : 01 < s GLYPH<20> 0 : 1 cm and s > 0 : 1 cm are green.  \nmissing mass is contained in the smallest particles, which trace the gas flows very closely. Having even the slightest turbulence in the primordial atmosphere (which is always the case, as discussed earlier) traps these particles and prohibits them from reaching the surface over sufficiently small time scales, but nonetheless, they eventually are accreted. That said, as the experiments approach a steady-state, the net accretion rate across a sphere of a given radius (e.g. R H or R B ) centred on the embryo approaches parity with the accretion rate onto the embryo. The two top rows of Fig. 14 also show how inefficiently the subsets of smaller particles are accreted. Very similar trends are present in the other two embryo mass cases we consider in this work, therefore,  \nthe figures for M = 0.5 M GLYPH<8> and M = 0.096 M GLYPH<8> , embryos are in Appendix A.  \nTo better illustrate the dependence of accretion on particle size, an animation showing the trajectories of pebbles in the vicinity of an M =0.95 M GLYPH<8> embryo ( m095t10 ) is available on-line as supplementary material (and on YouTube 4 ). In the animation, cyan represents 0 : 0 GLYPH<20> s GLYPH<20> 0 : 01 cm particles, red represents 0 : 01 < s GLYPH<20> 0 : 1 cm particles and yellow represents 0 : 1 < s GLYPH<20> 1 : 0 cm particles. Velocity streamlines are shown in green. The coloured circles meanwhile denote the Hill radius (magenta) and the 'nominal' effective accretion radii for unperturbed flow, r eff (Lambrechts & Johansen 2012, eq. 40; see below) for the three, aforementioned size bins. The animation clearly illustrates that: (i) particles that pass within their respective r eff are accreted; (ii) small, wellcoupled particles are particularly affected by the gas flow and often do not reach r eff .", '5.3 Quantitative measurements of accretion rates': "To allow selection of arbitrary sub-populations, our initial particle population is distributed all over space, with weights proportional to the local cell gas mass. To obtain accretion rates pertaining to a disk where particles have settled towards the midplane, we select the sub-population that initially resides within one Hill radius from the midplane, and subsequently measure accretion rates based on only those particles, renormalising their macro-particle weights w so they correspond to a total mass equal to the assumed dustto-gas ratio (0.01) times the total gas mass in the experiment.  \nBy sub-dividing the particle distribution into 6 size bins, we can study the dependency of the accretion rates on size. Fig. 15 shows that the net accretion rate onto the embryo scales linearly with the size s , which is different than the s 2 GLYPH<157> 3 scaling predicted by e.g. Ormel & Klahr (2010); Lambrechts & Johansen (2012), assuming unperturbed shear flow conditions. As acknowledged in Lambrechts & Johansen (2012), these conditions only hold when the mass of the planet embryo is very small. The key reason this scaling does not hold in our simulations is because the flow pattern (i.e. horseshoe orbits with stagnation points/separatrix surfaces) and density (with contrasts of the order of 10 3 ) in the vicinity of the planet are very different from unperturbed shear flow. Small, well-coupled particles are most strongly affected by these differences; they tend to follow gas streamlines, and therefore follow horseshoe orbits or are deflected around the planet. As a result, the dependence of the accretion rate on particle size becomes steeper than 2 GLYPH<157> 3 . The measured linear scaling is valid even for larger particles, all the way to dmsized particles, as illustrated in Fig. 16. The accretion rate saturates (no longer increases with particle size) for particles larger than about 1 meter.  \nA crude measure for the importance of flow effects for pebble accretion can be obtained by equating the pebble accretion impact radius in the shear limit ( r eff GLYPH<24> GLYPH<28> 1 GLYPH<157> 3 R H , where GLYPH<28> is the Stokes number) with the Bondi radius R B -the scale on which flow effects are expected to become important. For r eff < R B flow effects may be expected to change accretion  \nFigure 14. From left to right, the mass accretion rate across the Hill sphere, Bondi sphere and embryo surface, respectively, over the course of the m095t10 run, displayed as a Hammer projection. From top to bottom are the particle accretion rates for three size bins 0 : 0 GLYPH<20> s GLYPH<20> 0 : 01 cm, 0 : 01 < s GLYPH<20> 0 : 1 cm and 0 : 1 < s GLYPH<20> 1 : 0 cm. A longitude of 270 GLYPH<14> corresponds to the y -direction and a latitude of 0 GLYPH<14> corresponds to the disc midplane.  \n<!-- image -->  \nFigure 15. Pebble accretion rates onto the embryo, for particle sizes within each size bin, drawn from an initial distribution consisting of particles within GLYPH<6> R H from the midplane, for the three embryo masses, 0.95 M GLYPH<8> (black), 0.5 M GLYPH<8> (red), and 0.096 M GLYPH<8> (blue). The green line shows a mass accretion rate proportional to particle size, s .  \n<!-- image -->  \nrates that have been derived under the assumption of an unperturbed gas flow. Equating the radii, we find a threshold mass of M P GLYPH<3> GLYPH<24> h 3 GLYPH<28> 1 GLYPH<157> 2 M star , where h is the disk aspect ratio. Hence above M P GLYPH<3> flow effects may modify pebble accretion estimates based on assumptions of pure shear flow. This expression also states that smaller particles (smaller  \n2  \nFigure 16. Pebble and boulder accretion rates onto the embryo, for particle sizes between 0.1 - 100 cm, subdivided into 6 size bins, drawn from an initial distribution consisting of particles within GLYPH<6> R H of the midplane for the 0.95 M GLYPH<8> embryo mass case.  \n<!-- image -->  \nGLYPH<28> ) are more severely affected, qualitatively consistent with what we find. However, to explain the reduction quantitatively, i.e., to explain the linear dependence of GLYPH<219> M on GLYPH<28> that we find, requires a more thorough analysis, which will be considered in a future work.  \nThe runs with 56 million macro-particles give similar results as the corresponding runs with 2.6 million particles, with less than 10% difference in the accretion rates through  \nFigure 17. Local pebble accretion efficiency (the net accretion rate divided by the radially inward mass flux across a sphere of a particular radius) as a function of particle size, for spheres of several different radii, drawn from an initial distribution of particles within GLYPH<6> R H of the midplane for the m095t10 run.  \n<!-- image -->  \nthe Hill sphere and less than 1% difference at the Bondi sphere. Thus the results are rather robust, even when using only 2.6 million macro-particles. Naturally, however, the greater the number of particles, the smaller a fraction of the population one can sample from and still retain statistical significance; i.e., one has the ability to perform more focused studies of particle behaviour, as a function of size and initial position. The M = 0.95 M GLYPH<8> case, for example, initially has only about 20% of the macro-particles within one Hill radius of the midplane. Thus, to study midplane accretion rates it is useful to have a large particle population available.  \nAn estimate of the maximum accretion rate may be obtained by taking the dust-to-gas ratio times the average gas mass flux through the Hill sphere. For our M = 0.95 M GLYPH<8> embryo immersed in our 1/10 MMSN disk, this flux is GLYPH<25> 4 10 GLYPH<0> 3 M GLYPH<8> yr GLYPH<0> 1 . With a dust-to-gas ratio of 0.01, we thus estimate a maximum accretion rate of about 4 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 . This rate is exceeded by about a factor of five in the 0.3 - 1 m size bin in Fig. 16. By comparison, the 'Hill rate', GLYPH<6> p GLYPH<10> K R 2 H , is 3 : 9 10 GLYPH<0> 3 M GLYPH<8> yr GLYPH<0> 1 and nearly identical to the actual, 3-D gas mass flux.  \nThe accretion rates for the largest size bin ( GLYPH<24> 0.3-1 cm) are approximately 3 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 , 7 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 and 2 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 , for the three embryo masses, respectively (cf. Fig. 15). The accretion rates scale, as expected, approximately as the Hill rate time the Stokes number (proportional to M 2 GLYPH<157> 3 a GLYPH<0> 3 GLYPH<157> 4 with MMSN scalings).  \nFig. 17 illustrates the local efficiency of accretion of particles across particular spheres in the m095t10 simulation. The 'local accretion efficiency' is defined as the net accretion rate through a sphere of radius R , normalised by the accretion rate obtained by counting only the solids transport inwards through the same sphere (note that this definition differs from other uses of 'accretion efficiency' in the literature; cf. Lambrechts & Johansen 2012; Guillot et al. 2014; Ormel 2017). Evidently, a large fraction of all the particles entering the Hill sphere leave again without being accreted. As illustrated in Fig. 17, the local efficiency is particularly low (near 0%) for the lowest mass particles, which are well coupled to the gas and, due to temporal fluctuations in the accretion rate, can in fact temporarily show negative values. The local efficiency, meanwhile, rapidly grows for smaller radial distances, and at GLYPH<24> 30 R P the local efficiency is nearly  \n100%, except for particle sizes below 0.003 cm, which are still very well coupled to the gas there. Over size bins where the efficiencies are statistically well determined, they also scale approximately linearly with size. This is a consequence of the denominator (the unsigned, radially inward transport rate) being approximately independent of particle size (at a given radius) in our regime.  \nWe do not see significant systematic differences when different pressure bump parameters GLYPH<16> p are used for the same mass embryo. This is to be expected, since the results will only change when the radial settling has become significant. One would expect on the one hand a tendency for reduction of the accretion rates, because of the smaller y -velocities closer to the orbital radius of the embryo, and on the other hand a tendency for increase, due to the accumulation of large size particles. Ultimately, the accretion rates would be constrained by the radial influx of particles from larger orbital radii into the pressure trap.", '6 ACCRETION DRIVEN CONVECTION': "The high resolution close to the embryo surfaces enables us to resolve convection driven by accretion heating. Accretion of solids onto the embryo releases a lot of potential energy, which is converted to heat via the friction force. To simulate the effects of the heating due to the accretion, we added an extra heating term Q accr to one of our simulations, m095t10conv . Considering a shell of extent dr at some distance r from the embryo, the volume dV = Adr , where A = 4 GLYPH<25> r 2 , then has a potential energy difference (which is per unit mass) given by GLYPH<8> ' r + dr ' GLYPH<0> GLYPH<8> ' r ' . To get an energy per unit volume we need to multiply with the mass accreted per unit time:  \nQ accr = GLYPH<219> MGM p 4 GLYPH<25> r 4 ; (31)  \nwhere Q accr is the heating rate per unit volume and GLYPH<219> M is the mass accretion rate (which we set to 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 ). The heating rate, Q accr , is then added to the entropy equation (Eq. 4) via:  \n@ s @ t = ::: + GLYPH<26> Q accr P gas : (32)  \nFig. 18 illustrates the consequences of the addition of the accretion heating term in the m095t10-conv run. It shows extensive convective patterns, which considerably modify the gas flows around the embryo (left panel). The horseshoe orbits have somewhat receded, and the weak random flows that appear in non-convective runs are effectively overwhelmed by the much stronger convective flow patterns. The middle panel shows a vertical projection of the convective motions, which extend out to nearly 40 embryo radii. The vicinity of the embryo is no longer in a quasi steady state the gas flow patterns constantly change, and the interconnection between the primordial atmosphere and the disk is thus strengthened; i.e., no parts of the atmosphere are effectively isolated. The convection is very dynamic, with local velocities reaching GLYPH<24> 1 km s GLYPH<0> 1 , inducing sound waves when the convective cells reach the outskirts of the atmosphere. The direction of the outflows are constantly changing and now lack any significant constant patterns, but do have an  \n<!-- image -->  \n<!-- image -->  \nFigure 18. From left to right: xy slice in the midplane, yz slice at x = 0 of the entropy per unit mass and particle paths in the vicinity of the core in a single snapshot of accretion driven convection in the m095t10-conv run. For clarity, in the middle panel, the gas streamlines are omitted. In the right panel, black curves show traces for particles of size 0 : 1 < s GLYPH<20> 1 cm, blue 0 : 01 < s GLYPH<20> 0 : 1 and green 0 : 001 < s GLYPH<20> 0 : 01 cm, respectively.  \n<!-- image -->  \napproximately isotropic distribution of the velocity dispersion.  \nWhen we see or model convection at a stellar surface (e.g. in the form of 'granulation' on the Sun), the cellular patterns that develop are very much constrained by conservation of mass. In a stratification where mass density drops very rapidly with height, mass conservation forces the flows to 'turn over', from up to down, which again forces the horizontal size to be limited (lest the ratio of horizontal to vertical speed become very large). Here, we have a situation that is quite different, with a 'small' body at the centre of a large spherical volume. This situation allows the flows much greater freedom, with a possibility of expansion in all directions.  \nThe heating driving the convection in this study, computed assuming GLYPH<219> M = 10 GLYPH<0> 6 M GLYPH<8> yr GLYPH<0> 1 , is not consistent with the actual, about 20 times larger measured accretion rates (e.g. Fig. 15). Moreover, radiative energy transfer and a more realistic EOS would tend to change the properties of convection, so to reach true consistency, improvements are necessary. Nevertheless, it is possible to learn already from the current experiment, which can illustrate to what extent the convective motions are able to affect accretion rates relative to an analogous non-convective case.  \nThe right hand panel of Fig. 18 shows a small subset of high cadence paths for accreted particles, colour coded by size. Although their paths are affected, the largest particles are not significantly stirred up by the convective motions. The smaller the particles are, however, the more they are affected by the convective motions. Quantitative measurements show that the average of the unsigned mass fluxes through spherical surfaces where convection is active is much larger than in the non-convective case.  \nHowever, as illustrated in Fig. 19, we find that convection in the current case does not significantly reduce the net accretion rate onto the embryo. Formally, the local accretion efficiency (the net mass accretion rate as a fraction of the unsigned, radially inward solid mass flux; Fig. 20) is reduced significantly for smaller particles. This is because convective motions gives rise to larger unsigned rates of radial mass transport, and hence increase the denominator in the defi-  \nFigure 19. Pebble accretion rates in the accretion driven convection run ( m095t10-conv ). Particle sizes are split into 6 representative logarithmic size bins. The green line shows a mass accretion rate proportional to particle size, s .  \n<!-- image -->  \nFigure 20. Local pebble accretion efficiency in the accretion driven convection run ( m095t10-conv ). Particle sizes are split into 6 representative logarithmic size bins.  \n<!-- image -->  \nnition of local accretion efficiency. The convective motions also cause larger temporal fluctuations, and as a result we measure negative net accretion rates and a local accretion efficiency of less than zero at R B for some of the smallest size bins. The convective velocities in this case are generally on the order of the free fall drift velocities of the particles, as illustrated in Fig. 21, but the systematic drift of particles rel-  \nFigure 21. Free fall drag velocities of particles of size 1 cm (blue), 0.1 cm (green), 0.01 cm (red) and 0.001 cm (cyan). The velocity dispersion of the m095t10-conv run is shown in magenta.  \n<!-- image -->  \native to the gas continues, and since the convective motions do not result in a net transport of mass, their effect on the particle drift is a second-order effect. It is possible, however, that this could change if the convective motions were sufficiently resolved and included more turbulent, small-scale motions.  \nThe effect on the particle accretion rate is also likely to increase, and possibly become quite significant, when accretion heating is included self-consistently by including radiative energy transfer and a realistic EOS. The resulting, self-consistent accretion rates are likely to be lower than the current values, but the convective motions are likely to be stronger than in our current, nominal case.", '7 SUMMARY AND CONCLUSIONS': "In this paper, we have reported the results from high resolution, nested grid, 3D hydrodynamic and particle simulations which focus on solid accretion onto Mars- to Earth-mass planetary embryos. The primary goals of this work were to establish realistic three-dimensional flow patterns inside and outside the Hill sphere, to use these to determine accurate accretion rates, and to investigate how the accretion rates depend on embryo mass, pebble size, and other factors.  \nOur excellent spatial resolution and efficient time integration, using locally determined time steps, made it possible-using only a single, 20-core compute node per run-to perform the first-ever simulations of pebble accretion covering scales from several disk scale heights to a few percent of a planet radius, including atmospheres that selfconsistently connect to a surrounding accretion disk. It is also the first time that accretion heating driven convection has been modelled in a domain with such a large dynamic range.  \nWeconducted a series of simulations with different planetary embryo masses and strengths of a pressure trap. To be able to accurately measure accretion rates, and to study how they are affected by convection, we injected millions of macro-particles with representative sizes ranging from 10 micron to 1 m. The motion of the particles were followed, taking into account drag forces relative to the gas,  \nwhile feedback from the particle ensemble on the gas was neglected-consistent with the modest dust-to-gas ratios occurring in the models. The motions of a subset of particles were recorded with high cadence, allowing for their trajectories and histories to be studied in great detail.  \nOur main results may be summarized as follows:  \n- (i) Starting from approximately hydrostatic initial conditions, relaxation periods on the order of one orbit are sufficient to establish near quasi-stationary conditions in the vicinity of the Hill sphere.\n- (ii) The relaxed atmospheric structures do not differ dramatically from spherically symmetric hydrostatic atmospheres matched to surrounding disk conditions at the Hill radius.\n- (iii) Adiabatic and hydrostatic atmospheres with an ideal gas and GLYPH<13> = 1.4 do not differ substantially from atmospheres computed with a more realistic, tabular EOS (Tomida & Hori 2016).\n- (iv) The presence of convection changes conditions well inside the Hill sphere from near-steady flows to dynamically evolving flows, albeit stationary statistics (e.g. velocity dispersion as a function of radius) develop in less than an orbital period.\n- (v) The convective motions result in an enhanced mass exchange between different layers inside the Hill sphere, but do not significantly affect the mass exchange through the Hill sphere with the surroundings, nor do they affect the net accretion rates significantly.\n- (vi) In our fiducial 1/10 MMSN model, the (unsigned) gas mass flux through the Hill sphere for our 3D horseshoe flows is about 4 10 GLYPH<0> 3 M GLYPH<8> yr GLYPH<0> 1 , and thus agrees closely with the canonical'Hill rate', GLYPH<6> GLYPH<10> R 2 H = 3 : 9 10 GLYPH<0> 3 M GLYPH<8> yr GLYPH<0> 1 (obtained assuming GLYPH<6> = 170 g cm GLYPH<0> 2 and R H GLYPH<157> R GLYPH<8> = 231 ).  \n(vii) The gas mass flux, scaled with the dust-to-gas ratio, is reflected in similar rates of particles passing both in and out of the Hill sphere, with inflow increasingly dominate for increasing particle size.  \n(viii) As illustrated by Figs. 13 and 14, the particles that are not accreted and leave the Hill sphere have preferred paths along the separators of the horseshoe flow and are concentrated towards the midplane.  \n- (ix) In contrast, at the canonical Bondi sphere, the particle fluxes are mainly inward (except for the smallest particle sizes), and hence changes in the detailed structure of the layers close to the embryo (e.g. due to the EOS, including from evaporation of pebbles) are not important for the final accretion rates.\n- (x) The accretion rates for particles is observed to scale linearly with particle size, different than for unperturbed shear flow conditions, where it scales as s 2 GLYPH<157> 3 (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014). We attribute the steepening of the power-law to aerodynamic deflection of well-coupled particles under the higher densities and very different flow patterns in our simulations relative to unperturbed shear flow conditions.  \n(xi) We find that the accretion rate of 0.3 - 1 cm particles for an 0.95 M GLYPH<8> embryo immersed in a 1/10 MMSN disk, assuming a dust-to-gas ratio of 0.01 and assuming that the dust has settled to a midplane layer with height H < R H , is about 2 10 GLYPH<0> 5 M GLYPH<8> yr GLYPH<0> 1 . Since the accretion rates scale essentially linearly with particle size and inversely with gas  \ndensity, the accretion rate for the same size particles would be essentially the same in higher density disks; the larger supply of solids in a denser disk will tend to be cancelled by a correspondingly slower accretion speed.  \n(xii) For our 1/10 MMSN solar mass disk, the accretion rate saturates (no longer increases with particles size) for particles larger than about 1 m, while for denser disks the particles can be even larger. Thus, only for such large particle sizes will the accretion rates onto an embryo increase with the disk surface density. Conversely, if mm-size pebbles dominate the mass budget, the accretion rates depend only weakly on the disk surface density.  \nHaving such a high spatial resolution means that we could only cover a limited parameter space and run the simulations for a relatively small number of orbital periods. We have also relied on an ideal gas EOS, while in reality the high temperatures near the embryo surface would lead to the dissociation of molecular hydrogen (e.g. D'Angelo & Bodenheimer 2013), and therefore result in a much denser and somewhat cooler inner atmosphere. This would not significantly alter the accretion efficiency of the larger particles, however, as these effects occur at small radii from which these particles are not able to escape.  \nGiven the accretion rates determined in Section 5.3, we can estimate the growth time scales from low mass seeds to full planets, under the specified conditions (i.e. a gas surface density of 170 g cm GLYPH<0> 3 , a dust-to-gas ratio of 0.01, and solids from within H < R H . As per the results discussed above, to lowest order, the effects of disk density on mass supply and accretion speed cancel. The remaining dominant scaling is the dependence on embryo mass, which is approximately M 2 GLYPH<157> 3 , corresponding to a mass that grows approximately as t 3 . This implies that, as long as this approximate scaling persists, growth from a small seed takes on the order of three times the instantaneous mass divided by the instantaneous accretion rate.  \nWe conclude that if growth from low mass embryos to full size planets is dominated by accretion of 0.3 - 1 cm size particles, the accretion from a low mass seed would have taken of the order of 0.15 million years in the case of Earth, and of the order of 0.1 million years in the case of Mars. With chondrule size particles-if we take them to be 0.3-1 mm-the growth times are 10 times longer. These estimates assume a local ratio of dust-to-gas surface densities of 0.01. Increased dust-to-gas ratios, e.g., due to radial accumulation of solids in a pressure trap, would lead to correspondingly reduced accretion times.  \nDuring the growth phase, embryos would inevitably be surrounded by hot, primordial atmospheres of much larger mass than the current day atmospheres of the rocky Solar System planets. The decrease in disk densities required to stop pebble accretion would also result in removal of much of the primordial atmosphere as it expanded to match the decreasing pressure at its outer boundary (Nordlund 2011; Schlichting et al. 2015; Ginzburg et al. 2016). The isotopic signatures present in noble gases in the Earth's atmosphere may provide indirect evidence related to this (Pepin 1991, 1992a,b).  \nA decrease in gas and dust densities will lead to increased rates of radiative cooling, at first dominated by cooling in the vertical direction. The consequential, gradual flat-  \ntening of the hydrostatic structure inside the Hill sphere could lead to the formation of a circumplanetary disk with prograde rotation (cf. Johansen & Lacerda 2010). In order to model such an evolution, it is crucial to realistically include the major sources of heating, energy transport, and cooling: Accretion heating, the resulting convection, and radiative cooling while also using a realistic EOS (e.g. Tomida et al. 2013; cf. Fig. 1) and opacities.  \nTo this end, we are currently conducting simulations that include these effects. In the future, we will also model the accretion heating in a manner consistent with the solid accretion rate (cf. Sect. 6) and include pebble destruction via ablation (e.g. Brouwers et al. 2017; Alibert 2017).", 'ACKNOWLEDGEMENTS': 'We would like to thank the anonymous referee for in-depth and helpful comments that have helped to significantly improve the quality of this manuscript. The work of AP and ˚ AN was supported by grant 1323-00199B from the Danish Council for Independent Research (DFF). CWO is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). The Centre for Star and Planet Formation is funded by the Danish National Research Foundation (DNRF97). Storage and computing resources at the University of Copenhagen HPC centre, funded in part by Villum Fonden (VKR023406), were used to carry out the simulations presented here. The authors are grateful to Kengo Tomida and Yasunori Hori, who kindly provided the tabular EOS. AP is grateful to the Anton Pannekoek Institute at the University of Amsterdam for a three month stay during which part of this work was carried out. This research made use of the YT project for analysis and visualisation (Turk et al. 2011).', 'REFERENCES': 'Nordlund ˚ A., Ramsey J. P., Popovas A., Kuffmeier M., 2018,', 'APPENDIX A: ADDITIONAL FIGURES': 'Fig. A1 shows the details of the gas flow in the m05t00 run for the M = 0.5 M GLYPH<8> planetary embryo. The similarity with the higher mass case is apparently because embryos of such masses have well established primordial atmospheres, which are very close to being hydrostatic, leading to very similar flows in the vicinity of the Hill sphere.  \nFig. A2 shows the mass accretion rates over the course of the m05t10 run. The accretion rates are about a factor of two lower than in the m095t10 case, and the preferred accretion paths through the Hill sphere are once again closer to the midplane and the stagnation points (separators).  \nFig. A3 shows the mass accretion rates over the course of the m01t10 run. Note that the contour levels are modified, in order to enhance the features. The simulation shows a similar scenario to the cases with higher mass embryos, except the mass accretion rates are once again smaller. Also, there is effectively no loss of mass from R B . This is because  \nFigure A1. Same as in Fig. 6, but for the M =0.5 M GLYPH<8> embryo from the m05t00 run.  \n<!-- image -->  \n(a) Midplane.  \n<!-- image -->  \n(b) z = 10 R P .  \n(c) z = 25 R P .  \n<!-- image -->  \n<!-- image -->  \n(e) Vertical component of the vorticity in the midplane. (f) xz slice of the vertical velocity component at y = 0.  \n<!-- image -->  \nthe R B for such a small mass embryo is just GLYPH<24> 12.3 embryo radii, or just 4% of R H , and therefore actually deeper in the potential well of the embryo than in the higher mass cases.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \n<!-- image -->', '22 A. Popovas, ˚ A. Nordlund, J. P. Ramsey, C. W. Ormel': 'Figure A2. Same as Fig. 14, but for the m05t10 run.  \n<!-- image -->', 'Pebble dynamics and accretion onto rocky planets.': 'Figure A3. Same as Fig. 14, but for the m01t10 run.  \n<!-- image -->'}