The Three Dimensional Flow Field Around Planets on Eccentric Orbits

We investigate the properties of the hydrodynamic flow around eccentric protoplanets and compare them with the often assumed case of a circular orbit. To this end, we perform a set of 3D hydrodynamic simulations of protoplanets with small eccentricities ($e\leq 0.1$). We adopt an isothermal equation of state and concentrate resolution on the protoplanet to investigate flows down to the scale of the protoplanet's circumplanetary disk (CPD). We find enhanced prograde rotation exterior to the CPD for low planet masses undergoing subsonic eccentric motion. If the eccentricity is made large enough to develop a bow shock, this trend reverses and rotation becomes increasingly retrograde. The instantaneous eccentric flow field is dramatically altered compared to circular orbits. Whereas the latter exhibit a generic pattern of polar inflow and midplane outflow, the flow geometry depends on orbital phase in the eccentric case. For even the modest eccentricities tested here, the dominant source of inflow can come from the midplane instead of the poles. We find that the amount of inflow and outflow increases for higher $e$ and lower protoplanet masses, thereby recycling more gas through the planet's Bondi radius. These increased fluxes may increase the pebble accretion rate for eccentric planets up to several times that of the circular orbit rate. In response to eccentric motion, the structure and rotation of the planet's bound CPD remains unchanged. Because the CPD regulates the eventual accretion of gas onto the planet, we predict little change to the gas accretion rates between eccentric and circular planets.


INTRODUCTION
While classic 1D models (Mizuno 1980;Pollack et al. 1996) model core accretion as the quasi-static contraction of a surrounding envelope, recent simulations of protoplanets are painting an increasingly hydrodynamic and three-dimensional picture. Early 3D simulations (Machida et al. 2008;Ayliffe & Bate 2012;Tanigawa et al. 2012) identified a robust pattern of polar inflow coupled with midplane outflow. This meridional circulation is sourced by high latitude horseshoe flows that descend and penetrate deep into the atmosphere before outflowing in the midplane (Fung et al. 2015). This circulation occurs for a wide range of planet massesfrom super-Earths (Ormel et al. 2015) and giant planet cores (Lambrechts & Lega 2017) to gas giants (Schulik et al. 2019). The pattern continues to persist even with more sophisticated treatment of the thermodynamics (D'Angelo & Bodenheimer 2013). It has been argued that for super-Earth sized protoplanets, gas recycling associated with meridional circulation is dominant enough to stave off the transition to runaway accretion (Ormel et al. 2015;Cimerman et al. 2017). Even in simulations where a bound atmosphere emerges (Lambrechts & Lega 2017;Béthune & Rafikov 2019), recycling can alter the region's size, making it smaller than 1D models would assume. This suggests that the 3D hydrodynamic nature of the flow field should be taken into account in building a comprehensive theory of planet formation.
Though less robust than meridional circulation, a number of studies have also seen the emergence of a circumplanetary disk (CPD). This inherently multidimensional phenomenon tends to occur preferentially towards higher mass planets (Machida 2009;Ayliffe & Bate 2009;Tanigawa et al. 2012). More recent studies have shown that CPDs are also sensitive to the thermodynamics of the gas (Szulágyi et al. 2016;Szulágyi 2017;Fung et al. 2019), occurring most readily in isothermal atmospheres (Tanigawa et al. 2012;Wang et al. 2014;Fung et al. 2019). When and if a CPD forms, in turn changes the nature of accretion onto the protoplanet (Schulik et al. 2020) and relates to questions of moon formation.
All of these simulations however, adopt the assumption of a circular orbit. In this study, we investigate how robust this and other results of the 3D simulations are when the circular orbit assumption is relaxed and replaced with a low to moderate eccentricity (e ≤ 0.1). Whether we are justified in allowing for eccentric protoplanets requires examining the balance of torques which either damp or excite eccentricity. Low mass embedded planets and planetesimals are expected to have their eccentricity damped through interactions with the surrounding disk. In linear theory, gravitational torques exerted at exterior Lindblad resonances serve to increase a planet's eccentricity, opposite signed torques exerted at co-orbital Lindblad resonances are large enough to overwhelm the exterior torques and make the net evolution one of eccentricity damping (Artymowicz 1993). The damping time for planet with mass M p embedded in a disk with aspect ratio h p is (Artymowicz 1993): with additional scalings dependent on disk properties omitted and assuming e h p . For e h p , Papaloizou & Larwood (2000) suggest an additional scaling (e/h p ) 3 making the circularization times longer, but still shorter than the formation timescale for planets. Torques from numerical simulations of embedded planets further suggest that eccentricity is damped (Bitsch & Kley 2010;Fendyke & Nelson 2014). For planets sufficiently massive to open a gap, the co-orbital Lindblad contributions vanish, removing the main source of eccentricity damping (Goldreich & Sari 2003). At this point in the planet's evolution, the response is highly non-linear, requiring full hydrodynamic simulations to adequately compute the torques. Simulations suggest that planets can become eccentric in gaps (D'Angelo et al. 2006;Dunhill et al. 2013;Duffell & Chiang 2015;Muley et al. 2019). Though these re-sults are generally for higher mass planets than those investigated here, similar behavior could apply to lower masses given a sufficiently low viscosity.
Furthermore, the damping mechanisms described above assume a fairly quiescent and passive formation process. In all likelihood, the formation pathways are far more complex, allowing for modes of exciting eccentricity. Towards low masses for example, a protoplanet swarm in equlibrium with its protoplanetary disk is expected to produce characteristic eccentricities on the order of h p (Papaloizou & Larwood 2000). A stellar companion or partner giant planet interacting through secular or mean-motion resonances could also excite these moderate eccentricities. Even heating effects can result in eccentricities close to h p (Eklund & Masset 2017;Chrenko et al. 2017). Though the relative importance of these effects is uncertain, it seems unlikely that the circular orbit assumption applies to all protoplanets. In Section 5.3, with insight gained from our simulations, we will further call into question the assumption of damped eccentricity.
Few works have explored the modifications that eccentric protoplanets would make to modern theories of planet formation. First there is the work of ;  who perform 2D and 3D pebble integrations to understand how pebble accretion efficiencies might change given an accreting eccentric protoplanet instead of the often assumed circular one. Their work suggests an enhancement to the pebble accretion efficiency for low eccentricities e ∼ 0.05 with an eventual suppression as one goes to even higher eccentricity e ∼ 0.1. This does not however, take into account perturbations to the flow field arising from the planetary potential-perturbations which, at least in the circular case, have been shown to dramatically suppress the accretion of small pebbles. More recently, Mai et al. (2020) who perform 2D radiation hydrodynamic simulations of planets on highly eccentric e > 0.1 orbits. At such high e, the dynamics becomes dramatically different from circular models due to the emergence of a bow shock. What happens in the e < 0.1 regime however? This regime is physically interesting because it spans the transition from subsonic to supersonic orbital motion. Is the flow field altered by the planet as in the circular case or is it determined by eccentric motion?
To answer these questions we perform global 3D isothermal hydrodynamical simulations of planets on circular and eccentric (e ≤ 0.1) orbits and study how the flow field is altered. Our simulations end up being entirely complementary to those by Mai et al. (2020). We sacrifice their more thorough treatment of thermodynamics for the ability to study the inherently 3D flows. We work in a global framework whereas they adopt the local approximation. We end our simulations with e = 0.1, precisely where their models begin. In this way, we are able to begin to build a unified picture of eccentric protoplanets and test the relative importance of relevant physics.
Because the work on eccentric protoplanets has been fairly sparse, this study is largely exploratory. We explore the most information dense regions of parameter space where the planet is both transitioning from subsonic to supersonic (e ≤ 0.1) and transitioning from the linear to non-linear mass regime. This allows us to explore a maximal amount of interesting regimes with the fewest number of computationally expensive models. The drawback of course being that we are unable to explore any of these regimes in particular detail. We instead extrapolate based on the few models run in this transition regime. This strategy necessitates further confirmation of the extrapolated results by more focused studies but provides an exploratory glimpse of the how eccentricity can modify the 3D structure of the protoplanetary envelope.
We begin this paper with a description of our methods (Section 2). In Section 3 we discuss the circular orbit models and demonstrate their consistency with previous works. Section 4 gives the results of our eccentric models and compares them with the established patterns of circular orbit models. Finally in Section 5 we discuss potential implications of these eccentric models and questions of their applicability.

Equations and Code
To determine how protoplanet eccentricity might alter current theories of planet formation, we perform a series of 3D global hydrodynamic simulations of a protoplanet orbiting its host star through a protoplanetary disk. In a frame rotating with uniform angular velocity Ω 0 , the motion of nebular gas at position r is governed by the following fluid equations: where ρ, v, P , Φ are the standard fluid density, velocity, diagonal pressure tensor, and gravitational potential respectively. For the sake of simplicity and because protoplanets generically form circumplanetary disks (CPDs) under isothermal conditions (Fung et al. 2015), we adopt an isothermal equation of state, with sound speed c, to close the above equations. Because of this choice, these models are able to test whether CPDs are disrupted or altered by eccentric motion.
To solve equations (2)-(4), we employ the Athena++ code 1 (Stone et al. 2020). Because Athena++ does not natively solve the angular momentum conserving form of the equations in a rotating frame, we implement all the terms on the right side of equation (3) as source terms. Though implementing coriolis forces in this manner has been shown to produce inaccurate results for long-term integrations (Kley 1998), for simplicity we adopt the source term method and restrict our simulations to tens of orbits which is sufficient to study flows near the planet operating on fast dynamical times. This method remains untenable for studying processes such as gap-opening.
We solve the equations in a spherical coordinate system (r * , θ * , φ * ) centered on the star and placing the orbital plane at θ * = π/2. The angular frequency of the rotating frame is set equal to the constant Keplerian frequency of the orbit. For planet with semi-major axis a and star with mass M * , this aligns Ω 0 with the polar axis of the (r * , θ * , φ * ) system and sets Ω 0 = GM * /a 3 . Under this construction, the gravitational potential Φ takes the form: where M p is the planet mass and (r p , φ p ) denote the radial and azimuthal coordinates of the planet. is a gravitational softening length, added to prevent the gravitational acceleration from becoming disruptively large near the planet. The third term in equation (5), known as the indirect potential, is included for the sake of completeness but results in unnoticeable differences for these models.
To avoid the unnecessary computational complexity from integrating Kepler's equations, we approximate the eccentric orbital motion to first order in eccentricity e and restrict the simulations performed to e ≤ 0.1. Under this approximation the planet coordinates in the rotating frame take the form, With the aforementioned choice of Ω 0 , this eccentric motion reduces to epicycles around the fixed point at r * = a, φ * = 0 in the rotating frame.

Length Scales and Model Parameters
In the case of a planet orbiting a star and interacting hydrodynamically and gravitationally with a gaseous protoplanetary disk, there exist several established length scales characterizing the relevant physics. As we refer to these quantities heavily in this work, it is prudent to define and discuss them here. First, the planet's Bondi radius: defined here in terms of Newton's gravitational constant G, planet mass M p and sound speed c, arises as the natural length scale comparing the planet's gravity to the thermal state of the nebular gas. Similarly, the planet's Hill radius: defined here in terms of the planet semi-major axis a and the planet-to-star mass ratio q, arises as the natural length scale comparing the strength of planet's gravity to that of its host star. With r b as a length scale set by the planet-disk interaction and r h as a scale set by the planet-star interaction, it is natural to also consider some scale set by the star-disk interaction -namely the disk scale height: defined here as the ratio of sound speed to orbital frequency Ω. Because only two of these three scales are independent it is sensible to uniquely characterize a model by two parameters. For this work we use the planet's thermal mass q t defined by as one of these dimensionless parameters. Choosing q t as the dimensionless parameter offers unique advantages over other potential options. In particular, it categorizes the dimensional length scales of the problem into the following well-defined hierarchy: For the ordering in equation (12), where the thermal mass is small (i.e. subthermal), q t reduces to the only parameter necessary to characterize the problem (Machida et al. 2008). We ensure that our models cover this subthermal regime by simulating q t = 0.25 and incrementing q t until computational challenges prevent us from obtaining interpretable results (q t ≈ 2). In order to accurately capture as wide a range of q t as possible, these simulations were designed to be global in nature. In the global framework, uniquely specifying a model requires a second dimensionless parameter which we choose to be the disk aspect ratio evaluated at the planet's semimajor axis: Because these types of simulations are predominantly dependent on the choice of q t and not h p (Fung et al. 2015) so we fix this secondary parameter to h p = 0.05 for all models. The crux of this work lies in introducing a third dimensionless parameter -the eccentricity e -and determining how this new parameter alters circumplanetary flows. The eccentricity is chosen to investigate a wide range of behaviors as the planet goes from circular (e = 0) to subsonic eccentric motion (e 0.05) to supersonic eccentric motion (e 0.05).
An additional parameter, the gravitational softening length though not physically motivated, is necessary to avoid computational difficulties. The softening length can have a dramatic effect on planet scale flows and can prevent CPDs from forming entirely by suppressing rotational velocities . In a physical sense, one desires to make the softening length as small as the planet's physical radius while still large enough to avoid any computational troubles. To this end, we run a set of simulations with small = 0.015R b which is small enough to form a flattened CPD around the planet. To verify the robustness of any results we run an additional set of simulations with large = 0.1R b which is too large to form a flattened disk. We summarize the range of parameters explored in Table 1 and provide labels to uniquely designate each model. The full suite of models covers all permutations of the parameters listed and future references to particular models are indexed by their parameter labels. For example, q025e1 refers to model with q t = 0.25, h p = 0.05, e = 0.1, and = 0.015R b . All models are run for 20 orbits as this tends to be enough time for flow fields to reach steady state (or quasi-steady state for eccentric models).

Initial Conditions
We initialize the disk in axisymmetric hydrostatic equilibrium at time t 0 with the same density and velocity profiles as Fung et al. (2019). In the rotating frame they are: We set ρ 0 = 1 because the self-gravity of the gas is ignored. This initial profile corresponds to a vertically integrated gas surface density that scales as r −3/2 * . We also introduce the component of the potential due to the planet Φ p gradually with the following form: This form increases the potential from 0 at t = 0 to its full value after 2 orbits.

Domain and Boundaries
For all models the same domain is used. The grid extends from 3a/5 to 5a/3 in radius with a logarithmic spacing. For radial boundary conditions, the fluid variables in the ghost cells are fixed to the equilibrium values obtained from equations (17)-(18). Our refinement prescription for the mesh (see 2.3.3) introduces sufficient diffusion near the boundaries to render wavekilling boundaries unnecessary. Because the equations have vertical symmetry across the midplane, we save a factor of 2 in computation by modeling only half the domain in θ p . The polar grid is then set only from θ p = π/2 − 4H p to π/2 with a linear mesh spacing and reflecting boundary conditions on both boundaries. The azimuthal grid covers the full 2π with linear mesh spacing. With this domain, the root grid is set to have 64 cells in r * , 16 cells in θ * , and 512 cells in φ * to keep cells roughly square.

Resolution
To study flows near the planet within a global simulation of the full protoplanetary disk necessarily requires a refinement scheme to concentrate resolution on the scales of interest. Athena++ comes equipped with a block-based adaptive refinement scheme that can be leveraged for these types of simulations. In the blockbased approach, meshblocks made up of grid cells are refined in a self-similar manner based on some userdefined criteria. For our criteria, we set a minimum cell size desired over a control volume centered on the planet and allow the mesh to adjust itself to meet these requirements. Though this means the resolution cannot be controlled in regions outside the specified control volume, the number of cells per meshblock can be tuned to either relax or steepen how quickly the resolution changes as one moves away from the planet. For our fiducial resolution, the minimum cell width is set to W min = R b /256 within a control volume of N vol = 32 cells in each direction from the planet. This corresponds to a resolution of 1/W min = 256 cells per R b in the innermost r < R b /8. The meshblock size is fixed at N i = 16 cells for each direction. With this refinement prescription, the cell width W i in the the i direction at coordinate where x p is one of [r p , θ p , φ p ] depending on the choice of i. To verify the accuracy of the fiducial models, higher resolution simulations are also performed. The high resolution models simply double the resolution everywhere in the domain. The fiducial models with = 0.015R b are not fully converged near the planet for some quantities like the integrated mass (Section 4.1). In these cases, the converged hydrostatic models with large are used to inform our understanding. Nevertheless, we find that the character of flow and conclusions presented are robust with respect to all models as the resolution is increased.

Frames and Coordinates
In Section 2.1, we introduced a spherical coordinate system for the non-inertial rotating frame centered on the star in which we perform our simulations. In this computational frame, e = 0 planets lie fixed at the point (r * , θ * , φ * ) = (a, π/2, 0) while e > 0 planets follow a periodic epicyclic path about the same point. Because we are predominantly interested in the character of flows near the planet, in the following analysis it becomes advantageous to work in a coordinate system centered on the planet rather than the star. We construct a Cartesian system on the planet indexed by coordinates (x, y, z). We align this new frame to be curvilinear with respect to the inertial motion of the planet. In particular, in an inertial frame,ŷ is chosen to always be aligned with the instantaneous direction of orbital motion.ẑ is chosen to lie perpendicular to the orbital plane (by this design it points in the same direction as the starcentered polar axis). Finally,x is chosen to ensure the usual right-handedness of Cartesian (x, y, z). We will also frequently refer to this frame using spherical coordinates to discuss the rotational and radial aspects of the flow. The spherical coordinates are defined with the usual Cartesian-spherical transformation: and are unscripted (r, θ, φ) to distinguish them from the star-centered coordinates (r * , θ * , φ * ). With the same coordinate transformation, equations (22)-(24), we also define a set of Cartesian coordinates for the frame centered on the star (x * , y * , z * )

THE CIRCULAR ORBIT
There already exist a number of studies investigating the interaction of the circular orbit planet with its isothermal environment in 3D (Machida et al. 2008;Ormel et al. 2015;Fung et al. 2015;Kurokawa & Tanigawa 2018;Kuwahara et al. 2019;Béthune & Rafikov 2019). Within this body of literature and across various numerical implementations (e.g. PEnGUIn (Fung et al. 2014) vs. PLUTO (Mignone et al. 2007) vs. Athena++, spherical vs. Cartesian coordinates, global vs. local approximation), several key aspects of the flow remain consistently observed. We demonstrate here that, in the limit of a circular orbit, our results are consistent with theirs.
3.1. Density Figure 1 highlights the midplane density perturbation caused by the introduction of a planet, using model q1e0 as an example. On the global scale, we find the expected spiral waves (Ogilvie & Lubow 2002) and a shallow gap forming at the planet's orbit. As one zooms into the region near the planet, the density begins to increase rapidly, owing to the isothermal equation of state. Whereas the hs models take on the character of a hydrostatic symmetric sphere in the region r R b /2, the small models have non-negligible rotational support and a vertically flattened density structure near the planet. The steady state vertical density profile for q1e0 is displayed in Figure 2. Zooming into the inner R b /2 near the planet in the third panel, one can see the existence of low density polar region and a flattened high density CPD near the planet. This agrees with the work of Fung et al. (2019) where isothermal CPD are shown to be ubiquitous given a small enough .
A feature of these models, that we have not found sufficient reference to in most existing works is the presence of two low density vortices in the polar region above the planet. There is some evidence for these vortices in Fung et al. (2015) and Fung et al. (2017). They find vortices occurring at approximately the same location but weaker depending on the choice of code. Here they are more dramatic, perhaps owing to the smaller choice of . Because vortices create pressure minima at their core, the isothermal EOS produces corresponding density minima. As a result, these vortices can be seen in Figure 2 as the sliced low density lobes appearing near (x * , z * ) = (0.99a, 0.015a) and (x * , z * ) = (1.01a, 0.015a). Though these vortices exist for all models with small , they become more tenuous as q t decreases.

Flow Field
The most salient point to emerge from circular orbit models is a pattern of meridional circulation near the planet. Both the models here and in other works (Machida et al. 2008  find that the flow near the planet is one in which gas flows in through the poles and out through the midplane. This finding is, in turn, used to inform full evolutionary models of planet and satellite formation (Batygin & Morbidelli 2020). An example of the meridional circulation in these models is shown in Figure 2 as arrows in the lowermost panel. The same pattern of meridional circulation is present for all q t . The circulation is fed from high latitude horseshoe orbits which then descend toward the planet. The polar inflow then reaches supersonic speeds, once the gas descends to within ∼ R b /2 of the planet. Near the midplane the flow turns and spirals outwards exiting near the planet's Lagrange points. This can be seen in Figure 3 where we plot representative midplane streamlines for ciruclar orbit model q1e0. The blue curves correspond to the aforementioned outflows sourced from high latitude horseshoe orbits. The other families of midplane flows are also pictured in Figure  3 with horseshoe flows in red and green and the background disk flows in purple and orange. Together, these three families of streamlines define the well established pattern of midplane flow for 3D isothermal simulations.
The same meridional circulation prevents solids in the midplane from being readily accreted due to the outflowing nature of the gas there. Because solids tend to settle to the midplane, this so-called outflow barrier (Kuwahara et al. 2019), could insulate the planet from a significant reservoir of solids. This behavior has been confirmed in 3D simulations of pebble accretion (Kuwahara & Kurokawa 2020) where midplane particles with dimensionless Stokes number τ s q 2 t are partially prevented from accreting.
The flow field of our circular models is more explicitly provided in the next section where we make direct comparison with the eccentric models.

COMPARISON WITH THE ECCENTRIC ORBIT
Within several orbits, eccentric models reach a quasisteady state where the state is approximately periodic with period 2π/Ω 0 . Because the states of the eccentric models are the same modulo one orbit, it is simplest to refer to a particular time in a quasi-steady model by the phase of its orbit ψ. In what follows, our convention uses ψ = π/2 as the phase of periapse and ψ = 3π/2 as the phase of apoapse.

Density
While the eccentric motion is still subsonic (e 0.05), the density field in the eccentric models is not appreciably altered from the circular profile. When the eccentric motion becomes supersonic (e 0.05), a bow shock appears, making the density field a function of orbital phase. In the subthermal models, this bow shock tends to increase the overall mass and density by increasing the pressure downstream of the shock.
This phase-dependent density field is apparent in the midplane and vertically sliced density profiles provided in Figure 4. The bow shock can be seen in the q t = 1 models particularly at ψ = 3π/2 ( Figure 4) but because the flow is already sonic near R b in those circular models, the bow shock is not nearly as robust or influential as in lower q t models. The bow shock forms between R b /2 and R b and is nearly perpendicular to the direction of eccentric motion, meaning its orientation rotates a full 2π around the planet over the course of one orbit.
In Figure 5, we plot the total mass interior to R b as a function of time for the hs models. The masses for the e1hs planets tend to track the circular cases with some additional periodic variation, but in the q t = 0.25 and 0.5 cases, the masses of these low q t planet atmospheres show an average uniform systematic enhancement of ∼ 40%. We use the hs models here in particular because the mass measurements are not entirely converged for smaller . The mass in the unconverged cases is dominated by a contribution close to the planet (r < ), which we cannot properly resolve even with several additional levels of mesh refinement.

Rotation
In the circular case, gas well inside R b rotates in a prograde sense, and gas on the scale of R b merges with the background shear (which has a negative, or retrograde, sense of "rotation", even though it does not circulate the planet.) In the eccentric case, we find that the rotation of the CPD (r 0.05R b ) is not significantly affected by eccentricity. In both the circular and the eccentric cases, the CPD appears to remain bound. Though determining what gas is 'bound' in these types of simulations is difficult, the CPD gas is bound in a kinematic sense. By this we mean: orbit-averaged midplane radial velocities are directed towards the planet within the CPD and away from the planet exterior to the CPD. This result is remarkably consistent with those of Fung et al. (2019).
Exterior to the CPD, the rotation rate appears sensitive to eccentricity when the planet is subthermal. In this parameter space, eccentricity, while still subsonic, may increase prograde rotation, but it rapidly flips the rotation to retrograde once it crosses into supersonic.
We say subsonic eccentric motion may increase prograde rotation because we only detect significant spin-up in one of our models, q025e05, which also has the lowest planet mass in our parameter space. To understand this behavior, consider the acceleration the gas experiences due solely to the eccentric motion. In the frame of the planet, and to first order in e, this acceleration takes the form a e = 2aΩ 2 0 e (sin Ω 0 tx − cos Ω 0 tŷ) of forced harmonic motion. Gas in the circular case that orbits the planet at angular velocity Ω greater than Ω 0 , experiences the orbit averaged acceleration. Because the gas on the scale of the CPD does rotate with Ω Ω 0 and the orbit averaged a e vanishes, the rotation curve of the CPD remains unchanged between circular and eccentric models. Further out from the planet, Ω drops and the harmonic forcing of Equation (25) can operate on the gas in nontrivial ways. This case of spin-up can be seen in Figure 6 where we present the orbit averaged specific angular momentum as a function of radius for the q025 models.
Though the spin-up occurs for only one fiducial model, we suspect this behavior is more generic at lower q t and that its limited observance here is due only to our choice of parameter space. This hypothesis is primarily motivated by the hs models where a pattern of spin-up is more prevalent (Figure 6). Although these models with a large softening length are more akin to planets with large radii relative to R b (i.e. super-Earths), they can nonetheless be used to inform trends about other regions of parameter space. Of the hs models, q025e025, q025e05, and q05e05 exhibit increased prograde rotation with the increase being greater as one goes to higher e or lower q t . This suggests the existence of some criterion as a function of q t , e that determines whether an eccentric model reproduces the circular rotation curve or whether it shows an enhanced level of rotation. We suspect that a similar criterion exists for the fiducial, smaller models, but is more stringent due to the increased gravity. If the fiducial models are following the same pattern as the hs models, then q025e05 just happens to be the only spin-up model that falls into our searched parameter space and any models with q t < 0.25, e = 0.05 should also have spin-up. Testing of the fiducial circular models with a smaller value of returns the same rotation curve on the scale of R b . Therefore, we conclude that the flow is converged on the scale where spin-up is occurring. Explicit confirmation of the spin-up ubiquity in lower q t < 0.25 models however remains inaccessible with our current numerical setup.
In the supersonic case, the rotational behavior of the gas changes and the above discussion of spin-up does not apply. For superthermal planets, the orbit averaged rotational velocities show slight enhancement on the scale of R b but interior to 0.1R b , the circular rotation curves are reproduced, suggesting again that the CPD is unaltered. For the subthermal models, gas within R b has a lower magnitude of rotation than the circular orbits with the q025 models even rotating in the retrograde sense. This state of suppressed rotation is possibly due to gas losing angular momentum as it passes through the bow shock. It is uncertain whether the retrograde behavior continues for lower q t but the absolute value of the rotation curve in all models tracks the circular CPD rotation precisely.

Meridional Circulation
As the eccentricity increases, models get further from the established paradigm of meridional circulation. No longer is the planet's atmosphere a steady polar inflow that turns and exits through the midplane as in Figure  2. The changes to the flow in the eccentric case can be summarized as follows: • The flow is highly variable in both time and space.
• The eccentric flow is periodic with period 2π/Ω 0 .
• The inflow and outflow can be larger in magnitude than the equivalent circular orbit flows.
• Though there is never truly significant outflow through the pole, the dominant source of inflow can actually be through the midplane (instead of the pole for circular orbits).
• The earlier points only become more apparent as one goes to model with lower q t or higher e.     radial mass flux ρv r exhibits the standard pattern of inflow/outflow with outflow occurring in the midplane primarily along φ = 3π/4 and φ = 7π/4. Adding a small e simply moves the projected mass flux pattern laterally in epicycles. As e increases, the flow becomes highly variable. Orbit-averaging would lead to large cancellations between in-and out-fluxes. In this sense, the orbit averaged mass flux is misleading-the amount of gas processed through the planet's Bondi radius can be considerably larger than that implied by the average mass flux.
To truly describe the flow in these eccentric models it becomes necessary to consider the inflow and outflow independently instead of just the average. In Figure 8, this is done explicitly by averaging the inflowing and outflowing mass fluxes separately over azimuth. The result is a time series of the mass flux as a function of polar angle for a particular radius (R b /2). This Figure in particular highlights many of the essential differences found between eccentric and circular orbits. We see that the magnitude of the flow becomes highly variable and exhibiting periods of significant increase. The geometry is also changed, with the dominant inflow contribution actually coming from the midplane.
For models with higher q t or e, the inflow/outflow magnitudes increase, becoming several times larger than the circular model in some cases. This is shown in Figure  9 where we plot the amount of inflowing and outflowing mass through the Bondi radius averaged independently over the course of an orbit. These mass flow ratesṀ are also normalized by their circular equivalentsṀ 0 to give the relative enhancement to the mass of processed gas as a function of our dimensionless model parameters q t and e. We refer to this eccentric flux relative to the circular flux as the enhancement factor and find that it scales inversely proportional with q t and linearly with e. This scaling can be interpreted as follows.
If the velocity in eccentric models is set by the eccentric motion, it scales as v ∼ aΩ 0 e. Meanwhile, the velocities in the circular case at R b are set by the velocity of the background shear which is on the order of of v ∼ Ω 0 R b ∼ aΩ 0 q t h p at r = R b . Because the background densities are similar between eccentric and circular cases, dividing the two velocities gives the scaling for the measured increase in mass. The result is: where the prefactor of 2 is chosen empirically as an approximate value necessary to fit the data in Figure 9. This enhancement to the processed gas mass is then simply a consequence of the planet sweeping up more gas as it moves upon epicycles than it can on circular orbit.

Midplane Streamlines
In section 4.2.1 and particularly in Figure 6, it is evident that eccentric motion can enhance the rotational state of gas in the midplane. The eccentric motion also acts to enhance the radial motion of gas in the midplane.
In the eccentric cases, the inflow is no longer restricted to φ = π/4, φ = 5π/4 but is augmented by inflow from the planet plowing through the surrounding disk. This inflow changes orientation based on the phase of the orbit, traversing the full 2π in azimuth. These general characteristics of the midplane inflow/outflow can be seen in the top row of Figure 10 where we plot the radial velocity in the midplane for several phases of an eccentric model along with the corresponding circular model. To decent approximation, the flow field can be considered as just the sum of the background shear and the epicyclic motion of the planet. This approximation is shown as the middle row in Figure 10 and can be compared to the top row for an estimate of its accuracy. We also plot the residual in the bottom row, which represents the change to the flow field arising from the interaction with the planetary potential, which is typical for all eccentric models. In all cases, the planet induces additional outflow in the direction of epicyclic motion and additional inflow in the opposite direction. These flows act to oppose gas motion supplied by epicyclic motion, as one might expect from the fact that the planet is creating a steady atmosphere around itself. The magnitude of this induced inflow/outflow tends to increase with larger q t and larger e.
To capture the full midplane velocity field, we integrate tracer particle trajectories. We do this by postprocessing models q1e0, q1e05, q1e1. 10, 000 particles are seeded randomly in the region 0.65 < r * < 1.65 and −π/4 < φ * < π/4 but excluding a 1.5R b around the planet. The trajectories are then integrated for 2 orbits using the full time-dependent flow field. We only use a first order method for this integration but with a timestep equal to the Courant condition, updating the flow field every timestep. In Figure 11 we plot a subset of these randomly seeded particle trajectories in the frame of the planet. For a more faithful representation of pathlines, the number of pathlines in each panel is kept proportional to number of particles that entered the frame. Due to our short integration time, some par-ticle trajectories are "unfinished" -the integration ends shortly after they enter the frame. We replace these few pathlines with other randomly selected pathlines that exhibit a more complete trajectory in Figure 11.
Upon integrating particle trajectories, it becomes clear that some fluid elements are able to make it through the time dependent inflow/outflow and travel all the way from R b to circulate near planet. Two such particle paths are highlighted in Figure 11 in red. This behavior is not possible if the planet were on a circular orbit due to the robust outflowing nature of midplane gas near the planet. In orange, we highlight a couple of particles that were able to make it fairly close to the planet but were ultimately swept up by outflow. Blue and green pathlines show particles that entered R b but are promptly diverted away either by outflow induced by the planet or outflow from the changing epicylic motion.
These particle trajectories illuminate several other key points of the eccentric midplane flow field. Most notable is just how different the geometry of the paths are between circular and eccentric models. Even for e = 0.05, it becomes impossible to associate an eccentric pathline with one of the families of streamlines in the circular case ( Figure 3). It is possible that the dearth of pathlines for e = 0.05 just below the planet could be associated with the circular outflowing horseshoe and recycling streamlines however. We also draw attention to the fact that the e = 0.1 contains a higher density of pathlines than e = 0.05, which in turn has a higher pathline density than the circular, e = 0 case. Because the same number of pathlines were randomly seeded in each model and we keep the number of displayed pathlines tied to the number of particles entering the frame, this higher density is partially a consequence of the mass flux enhancement discussed in 4.2.2. The model with higher e sweeps up more material into R b as it undergoes its epicyclic motion.

Solid Accretion
Planets immersed in a Keplerian shear are expected to readily accrete solid pebble sized particles due to the combined effects of gas drag and the planetary potential. If the pebbles have settled into a disk thinner than ∼ R b , we expect eccentricity to increase pebble accretion rate.  show that when the relative velocity between the planet and a pebble with stopping time t s is set by eccentric motion (v ∝ e), it leads to a smaller accretion cross section b = GM p t s /v ∝ e −1/2 with increasing e. The net pebble mass flux however, is increased:Ṁ = ρbv ∝ e 1/2 . On the other hand, if the pebble layer is thick, the mass flux scales asṀ = ρb 2 v so that eccentric motion causes no change to the accretion rate. For sufficiently high eccentricity (e/h p (q t /τ s ) 1/3 ), pebbles are not sufficiently coupled to the gas and pebbles with dimensionless stopping time τ s , enter the ballistic regime . In the ballistic regime, the cross section is set by the planet's physical radius, leading to a severe suppression of the accretion rate. The result is an increase in the 2D pebble accretion rate by a factor of ∼ 4 for modest e ∼ 0.05 but not too large eccentricities.
This picture of eccentric pebble accretion does not include changes to flow field induced by the planetary potential however. Using our high resolution 3D simulations, we have been able to test whether the eccentric flow field is sufficiently determined by the eccentric motion or whether deviations arise from interactions with the planet.
A first and obvious amendment to the estimates put forth by  is motivated by the dramatic density profile surrounding the planet. In both circular and eccentric cases, the atmosphere density is orders of magnitude higher than the background disk. For a pebble in the Epstein drag regime, this density increase leads to an equivalent decrease in the stopping time t s ∝ ρ −1 and thereby a smaller accretion cross section b ∝ t 1/2 s ∝ ρ −1/2 . This skews the accreted pebble distribution toward larger pebble sizes. On the other hand, a smaller t s means a pebble is less likely to fall into ballistic regime, leading to a substantial increase in the accretion cross section for pebbles near the ballistic regime transition. With eccentric motion, we find evi-dence of mass increase interior to a bow shock which only adds to these effects. For pebbles experiencing Stokes drag, t s is independent of ρ, and these effects vanish.
For our steeply increasing density profile, pebbles able to make it close to the planet can have their stopping times reduced by 2-3 orders of magnitude. The fate of pebbles with sufficiently low τ s to make them wellcoupled to the gas becomes determined by the topology of the velocity field. In the circular case there is no hope to be accreted -they are simply diverted away by the horseshoe flows (Ormel 2013;Kuwahara et al. 2019) or recycling streamlines seen in our Figure 3. In the eccentric case, the flow topology changes drastically and judging from our red curves in Figure 11, there exist trajectories that advect particles significantly closer to the planet. If a pebble is advected by one of these trajectories, the possibility of accretion is increased, owing both to the increased gravity and the opportunity to circulate the planet for an extended period. We therefore expect the accretion rate of sufficiently low τ s pebbles to be greater than the circular case. On the other hand, we find the planet induces outflow always oriented to oppose the epicyclic motion ( Figure 10). Small inflowing particles can be slowed by this induced flow and then swept up by outflow as the planet changes the direction of its motion. At least with tracer particles, we see a significant fraction penetrate into the Bondi radius, become stalled near the planet induced outflow, and then exit as the epicyclic motion changes large scale inflow to outflow. This would suggest some suppression for the accretion efficiency predicted by  toward the small t s end. Because we still see some tracer particles making it close to the planet, the net result is still one of increased accretion relative to the circular case.
For larger particles and those in the Stokes regime, the small scale flow topology and density profile become irrelevant. The accretion of these particles is instead determined by the large scale flows set by the planet's eccentric motion. As we demonstrate in Section 4.2.2, these large scale fluxes are enhanced by a factor of up to ∼ 4, simply due to the eccentric orbit sweeping up more gas as it moves through the disk. A similar enhancement may be applied to the pebble accretion rate. Though they use different methods,  also find a pebble accretion enhancement of up to ∼ 4, likely because the underlying physics is similar. The details depend upon both pebble size and eccentricity, but our estimates show that the pebble accretion rate for eccentric planets should be a factor of a few greater than their circular counterparts. Figure 10. Midplane radial velocity vr in the planet frame for models q25e0 and q25e05. The top row shows vr in our simulations. The middle row shows vr one expects arising purely from the disk's background shear (Equation (18)) and the epicyclic motion of the planet (Eqns. (6)- (7)). The bottom row plots the residual of the first two rows, where positive (red) is planet-induced outflow and negative (blue) is planet-induced inflow.

Gas Accretion
Traditional 1D models of core accretion treat the evolution as a hydrostatic growth sequence whose contraction is regulated by the envelope's ability to cool. Our hydrodynamic simulations exhibit two departures from this classic picture: 1) A continuous recycling flux of unbound gas exterior to r 0.1R b and 2) A robust bound CPD occurring within r 0.1R b . Eccentric motion tends to increase the former while leaving the latter unchanged.
To the first point -circular orbit atmospheres both here and in previous studies (Ormel et al. 2015;Cimerman et al. 2017;Lambrechts & Lega 2017) exhibit significant recycling of unbound gas through at least some portion of their atmosphere. To quantify the rate of this recycling, Cimerman et al. (2017) define a recycling time as the mass of the envelope interior to some radius divided by the inflowing mass flux: t recycling = M (< r)/Ṁ in (r). When this recycling time is less than the cooling time for the gas, gas cannot cool efficiently. In our eccentric models, the recycling time can be made significantly smaller by the mass flux enhancement discussed in Section 4.2.2. Because the mass of an eccentric model is roughly the same as its circular equivalent, the eccentric t recycling can be obtained simply by dividing the circular recycling time by our Equation (26). If the mass flux enhancement varies as a function of radius then this should also be taken into account. However we find that eccentric mass flux enhancement also tends to be fairly constant over a large range of radius. The Figure 11. Midplane pathlines for models q1e0, q1e05, and q1e1. Color is chosen to highlight particular pathlines (see text for details). extent of the recycling region and its importance varies depending on the equation of state and how one defines gas as 'bound'. In our case, recycling is best defined as occurring exterior to the CPD, r 0.1R b . In any case, eccentric motion should only serve to hasten the recycling of gas irrespective of one's precise definition of the recycling region.
To the second point -interior to 0.1R b , the flow is dictated by the CPD. Similar to Fung et al. (2019), the gas in our eccentric CPDs is bound in a kinematic sense -orbit-averaged midplane radial velocities are directed towards the planet within the CPD and away from the planet exterior to the CPD. In this case, the eventual accretion of gas onto the planet would be regulated by details of the CPD with the recycled outer envelope acting as an outer boundary condition. In this sense, eccentricity only matters for setting the boundary condition on the CPD.
If eccentricity were able to modify the rotation curve of the CPD, it would also directly modify the accretion rate, but as we show in Section 4.2.1, the CPD rotation curve remains insensitive to changes in e. This is at least partially a consequence of our choice of an isothermal equation of state. Under a different equation of state, rotational velocities tend to be lower than in our isothermal models ) and often fail to form CPDs until well above our mass range (Schulik et al. 2020). In those cases, where gas is rotating more slowly, the rotation curve could be more dramatically modified by eccentric motion (akin to our hs models). If the rotational velocities were increased/decreased (as in our subsonic/supersonic case), the accretion rate could be suppressed/enhanced. But because the 3D kinematics and thermodynamics are inextricably linked, simulations including the proper prescriptions for cooling are required to address these questions.

Eccentricity Damping
With some idea of the flow field around eccentric planets, we call into question the often made assumption that sub-Jovian planets should have their eccentricity damped. As we have seen, our model planets are able to significantly perturb the flow field in non-linear and nontrivial ways. Though one can attempt to estimate nonlinear torques in the method of Ward (1991) by including the angular momentum exchanged between the planet and horseshoe streamlines, our simulations demonstrate their inadequacy. Figure 11 shows that the notion of horseshoes, separatrices, families of streamlines, all fail to apply to even the modest eccentricity of e = 0.05. This suggests a need for full 3D simulations studying the eccentricity evolution of protoplanets. To our knowledge however, only one such study exists - Bitsch & Kley (2010). Though their work finds a general pattern of eccentricity damping, their models suffer from a resolution of ∼ 5 cells per R b and the excision of a significant fraction of R b . Though this excision was a sensible effort to disregard contributions from gravitationally bound material, our models show that most of this material is actually unbound. Judging from the complicated flow pattern of our Figure 11, it seems that ∼ 5 cells per R b is insufficient to resolve the torque contribution near the planet. With no other such simulations performed in the past decade, we find insufficient evidence to conclude that sub-Jovian protoplanets should have their eccentricities damped. We are working on calculating the torques for our own simulations in a self-consistent way, in order to determine the resulting orbital evolution.

SUMMARY
We have carried out a set of global 3D isothermal simulations of eccentric planets with e ≤ 0.1 and thermal mass 0.25 < q t < 2.0 corresponding to 10M ⊕ < M p < 80M ⊕ around a 1M star. Our main findings are: 1. Supersonic eccentric motion leads to the formation of a bow shock. This increases the overall mass interior to the shock (Section 4.1).
2. The structure and rotation of the CPD is resilient, remaining kinematically bound and similar to the circular orbit models even with a moderate eccentricity (Section 4.2.1).
3. Exterior to the CPD, subsonic eccentric motion leads to increased prograde rotation. This occurs for subthermal models with sufficiently large yet still subsonic e. As the motion becomes supersonic, the rotation in subthermal models becomes increasingly retrograde (Section 4.2.1).
4. Flow geometry, now dependent on orbital phase, is altered from the established circular picture. Instead of a pattern of polar inflow/midplane outflow, it can vary wildly. In some cases, the amount of midplane inflow becomes greater than the polar inflow (Section 4.2.2).
5. Flow magnitudes and fluxes are amplified with velocity set by the eccentric motion. The flux of material processed through the Bondi radius is enhanced by a factor ∼ 2e/(q t h p ). Recycling times exterior to the CPD are decreased by this same factor (Section 4.2.3).
6. The increased flux can increase the pebble accretion rate up to several times the circular rate. Planetary perturbations to flow field can suppress this accretion rate for the smallest solids, but the net result is still one of increased pebble accretion (Section 5.1).
7. Because the CPD is unaltered by eccentric motion, the rate of gas accretion onto eccentric planets should be comparable to the circular case (Section 5.2). Though we make some extrapolations to lower mass protoplanets, a global framework limits the ability to explore the subthermal regime more thoroughly. Future models, employing the local approximation, will be able to test whether the conclusions made here still hold when extrapolated to lower q t . Future models will also implement a more sophisticated treatment of the thermodynamics. This will undoubtedly change the nature of the bow shock interaction among other things. Additionally, though we make some qualitative predictions about the accretion of solids, quantitative verification requires integrating solid trajectories with the included effects of gas drag. With such a time-dependent flow, it becomes necessary to integrate the particles in conjunction with the hydrodynamics. With the development of a dust module in the Athena++ code, we will be able to more adequately address this question of solid accretion onto eccentric protoplanets. Regardless, this paper serves a departure from the oft assumed circular orbit and demonstrates that small to moderate eccentricities can alter established characteristics of planet formation.