The Dynamic Proto-atmospheres around Low-mass Planets with Eccentric Orbits

Proto-atmosphere formation on planets with eccentric orbits (e ≥ 0.1) can be modeled as a gravitating spherical object moving supersonically in the inviscid, compressible, and homogeneous gas. In this study, we include radiation transport of energy and assume the disk gas to be an ideal gas. We adopt the assumption of local thermal equilibrium and the linear two-temperature approach⁶ for the flux-limited diffusion (FLD) approximation for radiation transfer (Kuiper et al. 2010; Kuiper & Hosokawa 2018). For convenience, we adopt a frame comoving with the planet in all simulations. The planet object is placed in a wind tunnel instantaneously at the beginning and the subsequent gas evolution is simulated using a radiation-hydrodynamics model until a steady/quasi-steady state is reached.

The problem can be studied with sophisticated hydrodynamics codes (e.g., PLUTO, see Section 2.2) solving the Euler's equations representing the conservation of mass, momentum, and energy:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\rho {\boldsymbol{v}}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}}\otimes {\boldsymbol{v}})+{\boldsymbol{\nabla }}P=\rho {{\boldsymbol{a}}}_{\mathrm{ext}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\boldsymbol{\nabla }}\cdot [(E+P){\boldsymbol{v}}]=\rho {\boldsymbol{v}}\cdot {{\boldsymbol{a}}}_{\mathrm{ext}}\end{eqnarray} \tag{ 3 }$$

where ρ is the mass density of the gas, ${\boldsymbol{v}}$ is the relative velocity between the planet and the gas, P is the gas pressure, ${{\boldsymbol{a}}}_{\mathrm{ext}}$ is the acceleration caused by external forces, E = E_int + E_rad + E_kin is the total energy density of the gas, which is the sum of internal energy density (E_int = c_vρT, since we are using a constant c_v; c_v is the heat capacity at constant volume, T is gas temperature), radiation energy density (E_rad = aT⁴, a is the radiation constant), and kinetic energy density (E_kin = (1/2)ρv²).

With operator splitting we separately solve the transport term ${\boldsymbol{\nabla }}(E{\boldsymbol{v}})$, leaving the time derivatives of the internal and radiation energy density as in Kuiper et al. (2010) and Cimerman et al. (2017):

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({E}_{\mathrm{int}}+{E}_{\mathrm{rad}})}{\partial t}=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{F}}-P{\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{F}}$ is the flux of radiation energy density. Substituting E_int by c_vρT and E_rad by aT⁴, we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{rad}}}{\partial t}=-{\left(\displaystyle \frac{{c}_{{\rm{v}}}\rho }{4{{aT}}^{3}}+1\right)}^{-1}({\boldsymbol{\nabla }}\cdot {\boldsymbol{F}}+P{\boldsymbol{\nabla }}\cdot {\boldsymbol{v}})\end{eqnarray} \tag{ 5 }$$

and from the FLD approximation

$$\begin{eqnarray}&&{\boldsymbol{F}}=-\displaystyle \frac{\lambda c}{\rho {\kappa }_{{\rm{R}}}}{\boldsymbol{\nabla }}{E}_{\mathrm{rad}}\end{eqnarray} \tag{ 6 }$$

where λ is the flux limiter following the choice in Levermore & Pomraning (1981), c is the light speed, and κ_R is the Rosseland mean opacity.

Note that the above calculations do not include irradiation from the central planet because the external sources of luminosity from the planet surface (due to radiogenic heat and/or accretion of planetesimals) are negligible in our setup. The heat generated by accretion of planetesimals with a mass accretion rate of 10⁻⁷ M_pl yr⁻¹ (Stökl et al. 2015) is three to four orders of magnitude lower than that from gas accretion in our models (M_pl is the planet mass). The radiogenic energy with a reference level of 10²¹ erg s⁻¹ ${\text{}}{M}_{\oplus }^{-1}$ (for early Earth, Stacey & Davis 2008) is at least two additional orders of magnitude lower (Stökl et al. 2015). For a more detailed analysis of the FLD approximation in radiation transport please refer to Kuiper et al. (2010).

The external forces considered in the gas system include the gravity from the planet object:

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{\mathrm{ext}}=-\displaystyle \frac{{{GM}}_{\mathrm{pl}}}{{r}^{2}}{{\boldsymbol{e}}}_{{\rm{r}}}\end{eqnarray} \tag{ 7 }$$

where G is the gravitational constant and r is the distance to the center of the planet. Following Cimerman et al. (2017), we do not consider ionization or dissociation of molecules as we focus on solving the flow structure and bulk atmospheric properties. In fact, the models in Desch & Connolly (2002) and Morris & Desch (2010) found that ∼7% of H₂ in nebular shocks can be dissociated, but the total energy of the gas remains unchanged. The radiative forces are neglected because they are small compared to gravity from the planet and the gas pressure.

The set of Euler's equations (Equations (1), (2), and (3)) is closed by the perfect-gas equation of state:

$$\begin{eqnarray}&&P=(\gamma -1){E}_{\mathrm{int}}=\displaystyle \frac{\rho }{\mu {m}_{{\rm{H}}}}{k}_{{\rm{B}}}T\end{eqnarray} \tag{ 8 }$$

where γ is the adiabatic index, which we fix as 7/5 (for molecular hydrogen) throughout the simulations. μ is the mean molecular weight, for which we adopt a constant value of 2.353, corresponding to the gas of solar metallicity. k_B is the Boltzmann constant and m_H is the mass of a hydrogen atom.

The initial values of ρ, P, and ${\boldsymbol{v}}$ of gas in the wind tunnel are set as ${\rho }_{\infty }$, ${P}_{\infty }$, and ${v}_{\infty }$, the same as the properties of unperturbed disk gas far away from the planet. The Mach number is defined as ${ \mathcal M }={v}_{\infty }/{c}_{\infty }$, where the sound speed of the disk is ${c}_{\infty }=\sqrt{\gamma {P}_{\infty }/{\rho }_{\infty }}$.

2.2.1. The Grid

2.2.2. Boundary and Initial Conditions

The boundary conditions of the computational domain are similar to those in Thun et al. (2016). The inner radial boundary of the domain is set as reflective so that no mass or radiation is transported through the planet surface. The outer radial boundary condition is configured according to where the gas flows in and out of the domain. When 0 ≤ θ < π/2, we use the inflow boundary condition to set the ghost cells with the unperturbed gas values. When π/2 ≤ θ ≤ π, a zero-gradient condition is implemented instead to ensure the values of the ghost cell stay the same as the boundary cells so no reflections would occur. We adopt axisymmetric boundary conditions for both boundaries in the polar dimension.

We use the well-known "minimum mass solar nebula" (MMSN) model (Weidenschilling 1977; Hayashi 1981; Thommes & Duncan 2006) to provide the approximate initial values for the gas properties in the unperturbed disk. In our canonical simulation, we set up a planetary object (M_pl = 1 M_⊕, R_pl = 1 R_⊕) embedded in the homogeneous disk gas (${\text{}}{T}_{\infty }=300$ K, ${\rho }_{\infty }={10}^{-9}$ g cm⁻³) located approximately at a = 1 au from the Sun. The velocity of the planet relative to the gas ${v}_{\infty }$ is set as 4.8 km s⁻¹ (or ${ \mathcal M }=3.93$). Planets with orbital eccentricities between 0.17 and 0.35 are able to reach such a velocity at certain positions in their orbits. Adopting values provided from different disk models would slightly alter the background density for the location of the planet. For example, the updated MMSN model (Desch 2007) would increase the nebular density at 1 au by 1–1.5 orders of magnitude, while the "minimum mass extrasolar nebula" (Chiang & Laughlin 2013; Lee et al. 2014) would increase it by 0.5–1 orders of magnitude. These variations are within the range of background gas density considered in our parameter study (Section 3.2.3).

More accurately, the density profile of the background gas medium follows an approximated vertical hydrostatic stratification structure for protoplanetary disks (e.g., Armitage 2010):

$$\begin{eqnarray}&&{\rho }_{\infty }(z)={\rho }_{\infty }\ \exp \left(-\displaystyle \frac{{z}^{2}}{2{H}^{2}}\right)\end{eqnarray} \tag{ 9 }$$

where z is the vertical distance from the midplane, $H\,=\sqrt{{{ \mathcal R }}_{{{\rm{H}}}_{2}}{{Ta}}^{3}/({{GM}}_{\odot })}$ is the pressure scale height of the disk gas (${{ \mathcal R }}_{{{\rm{H}}}_{2}}$ is the specific gas constant for H₂, M_⊙ is the solar mass). For simplicity, we do not take into account the vertical density profile but adopt the midplane gas density ${\rho }_{\infty }$ as the background value. This is because we concentrate on relatively small planets or planetary embryos (e.g., objects of sub-Earth or Earth size). The radii of these objects (including the proto-atmospheres) are at least four orders of magnitude smaller than H in our simulations.

The initial setup is slightly different in our orbit-dependent simulation (see Section 2.3) because we start the simulation when the planet is at its perihelion. The initial values of gas properties (${\rho }_{\infty },{T}_{\infty },{P}_{\infty },{v}_{\infty }$) are adjusted to match the corresponding disk environment for an orbit with semimajor axis a = 1 au and eccentricity e = 0.2.

2.2.3. Opacities

The main goal of this work is to explore how an eccentric orbit affects the dynamics and structure of the planetary proto-atmosphere. To achieve this goal and answer the four science questions in Section 1 we carry out a list of simulations with different initial conditions, which can be divided into the "snapshot" simulations and the orbit-dependent simulation (see Table 1).

Table 1.  List of Hydrodynamic Simulations

"Snapshot" Simulations
Name	${v}_{\infty }$ (km s−1)/${ \mathcal M }$	Mpl (M⊕)	${\rho }_{\infty }$ (g cm−3)	Comments
Canonical	4.8/3.93	1.0	10−9	The canonical case
VEL3.2	3.2/2.64	1.0	10−9	Sensitivity study on Mach number/relative velocity
VEL3.7	3.7/3.02	1.0	10−9
VEL4.2	4.2/3.47	1.0	10−9
VEL5.3	5.3/4.38	1.0	10−9
VEL5.9	5.9/4.84	1.0	10−9
VEL6.4	6.4/5.29	1.0	10−9
VEL7.0	7.0/5.75	1.0	10−9
VEL7.5	7.5/6.2	1.0	10−9
PM0.1	4.8/3.93	0.1	10−9	Sensitivity study on planet mass
PM0.3	4.8/3.93	0.3	10−9
PM0.5	4.8/3.93	0.5	10−9
PM0.7	4.8/3.93	0.7	10−9
PM2.0	4.8/3.93	2.0	10−9
PM3.0	4.8/3.93	3.0	10−9
DEN-7	4.8/3.93	1.0	10−7	Sensitivity study on gas density
DEN-8	4.8/3.93	1.0	10−8
DEN-10	4.8/3.93	1.0	10−10
Orbit-dependent Simulation
Name	Mpl (M⊕ )	a (au)	e	Comments
ORB	1.0	1.0	0.2	Time-dependent simulation to track orbital evolution

Download table as:  ASCIITypeset image

For the "snapshot" simulations, we set the relative velocity between the planet and disk gas (${v}_{\infty }$, or represented by the Mach number ${ \mathcal M }$), the planet mass (M_pl), and gas density (${\rho }_{\infty }$) as free input parameters to specify the disk–planet environment. The gas pressure is adjusted to match the gas temperature at 300 K. The simulations are run until they reach steady/quasi-steady states so that the solutions represent "snapshots" of the proto-atmospheres in the corresponding environments. For the orbit-dependent simulation, we use the planet mass (M_pl), the orbital semimajor axis (a), and orbital eccentricity (e) as inputs to determine the specific eccentric orbit for the planet. Other parameters such as gas density, pressure, temperature, and relative velocity are determined upon the time-dependent star–planet distance, which in turn is dependent on these input orbital parameters. It is therefore a time-dependent simulation in which we read outputs regularly to track the evolution of the proto-atmosphere throughout the eccentric orbit. Given that these simulations are significantly more expensive in terms of computation hours, we do not perform a parameter study involving full-orbit simulations in this work.

In reality, the orientation of the planetary bow shock is also changing as the planet moves to different locations on the orbit. It is directly determined by the vector difference between the planet orbital velocity (${{\boldsymbol{\upsilon }}}_{{\rm{pl}}}$) and the disk gas velocity (${{\boldsymbol{\upsilon }}}_{{\rm{gas}}}$) at each location (see Figure 2). Both quantities decrease in magnitude with a larger distance from the star. At perihelion, the planet travels faster than the gas and generates a shock front aligned with ${{\boldsymbol{\upsilon }}}_{{\rm{pl}}}$. At aphelion, the gas in sub-Keplerian motion is faster and the shock front is in the opposite direction as ${{\boldsymbol{\upsilon }}}_{{\rm{pl}}}$. As shown in Figure 2, the change of bow shock orientation throughout the orbit is limited (max. ∼20° for an orbit e = 0.2). Because the timescale on which the orientation changes is much longer than the dynamical timescale of gas accretion, the variation of bow shock direction is negligible in the local-framed hydrodynamic simulations.

Zoom In Zoom Out Reset image size

3.1.1. The Flow Structure

When given enough time, all the "snapshot" simulations are able to reach steady states where the atmospheres reach their final mass and the accretion rates approach zero. Note that we do not take into account the contraction of the atmosphere occurring on the Kelvin–Helmholtz timescale, which is too long to be modeled in radiation-hydrodynamic simulations. Depending on the specific setup, the time required to reach the steady state varies. To avoid extensive computational time for some simulations to reach steady states, they are run only until they reach quasi-steady states.⁷ The time required for the system to reach a quasi-steady state is also the dynamical timescale. For the canonical case, the timescale is ∼6 × 10⁶ s. We then fit the temporal growth of the atmosphere mass with an asymptotic function to estimate the approximate final mass (see Appendix for a detailed description).

In Figure 3 we present the hydrodynamic solutions (gas density and temperature) of the canonical simulation in a quasi-steady state (M_pl = 1.0 M_⊕, ${\rho }_{\infty }={10}^{-9}$ g cm⁻³, ${v}_{\infty }=4.8$ km s⁻¹). Panels (a) and (b) show the 2D solutions while (c) and (d) show the 1D profiles along different angles.

Zoom In Zoom Out Reset image size

A strong bow shock is formed in front of the planet with the shock front located at 3.6 R_pl (planet radii). The ram pressure from the nebular gas and the bow shock prevent a large amount of gas from being accreted to the planet on the side directly facing the shock front. But the planet still retains a proto-atmosphere from the post-shock gas that is largely hydrostatic in a steady/quasi-steady state. The proto-atmosphere is largely pressure-supported. The nebular gas is heated at the shock front due to strong compression and preheated beyond the shock front by the radiation emitted by the hot post-shock gas. This is a type of radiation shock called a supercritical shock (Sincell et al. 1999) and features a temperature spike at the density shock front (Figure 3(d)). The stronger the shock is, the higher the temperature spike can be and the larger the preheated region beyond the shock front is. Figure 3(d) also shows that the inner proto-atmosphere is mainly heated by the accretional energy, and the temperature of the outer layer (two to three planet radii) is mildly raised by the compression of the bow shock. There is also a small high-temperature region behind the planet caused by the converging flows from vortices. We illustrate the gas flow with the linear integral convolution (LIC) (Figure 3(a)) and streamlines (Figure 3(b)). Both illustrations show a clear boundary between two types of flow—gas entering the bow shock with larger impact parameter and flowing directly past the planet, and gas with relatively small impact parameter flowing back toward the planet in the downstream region, forming multiple stable vortices ("egg-shape" region, enclosed by the yellow dashed line in Figure 3(a)).

The "egg-shape" flow region resides in the post-shock subsonic area, with its boundary in front of the planet coinciding with the stagnation point of the bow shock. Interestingly, the proto-atmosphere is accreted not directly through the gas encountering the front of the planet, but through the gas flowing back toward the planet in the "egg-shape" pattern. Such patterns are clearly distinct from those of gas accretion on planets with circular orbits, where gas is accumulated from all directions (e.g., Ormel et al. 2015a, 2015b; Cimerman et al. 2017). The flow structure becomes stable soon after the bow shock is established.

As the simulation reaches a quasi-steady state, we dye the proto-atmosphere region with a tracer to track the gas movement over ∼5 × 10⁶ s, a time long enough to see the gas replacement inside the "egg-shape" region (Figure 4). Although the tracer diffuses over time, it is clear that the majority of the dyed gas flows back to the planet inside the egg-shape region. A portion stays bound close to the planet, while the rest eventually flows out of vortices. These vortices then get replenished by the next batch of incoming gas. The gas inside the "egg-shape" region is constantly recycled and exchanges with the disk gas, with an estimated recycling timescale between 5 × 10⁵ and 5 × 10⁶ s,⁸ about half of the dynamical timescale.

Zoom In Zoom Out Reset image size

3.1.2. Definitions of the Atmosphere

Defining the boundary of the proto-atmosphere of an eccentric planet is far less straightforward than in the cases where planets are on circular orbits. The commonly used definitions for symmetric atmospheres such as the Hill radius and Bondi radius are no longer applicable in the supersonic environment. Like those around planets on circular orbits, the proto-atmosphere is always connected to and actively interacting with the disk gas (Section 3.1.1).

For the purpose of calculating the atmosphere mass, we provide three different ways to define a boundary for the proto-atmosphere based on the flow structures resolved by the hydrodynamic simulation:

• 1.  
  A gas density contour around the planet that encloses the hydrostatic atmosphere and is not too distorted in shape because of the bow shock. We found that the 10^−8.5 g cm⁻³ contour is most suitable when the background gas density is ${\rho }_{\infty }={10}^{-9}$ g cm⁻³. The mass enclosed by this boundary is 1.40 × 10²¹ g.
• 2.  
  The boundary of the "egg-shape" region where the gas recycles in and out. There is 1.43 × 10²¹ g of gas in this region when the system reaches a steady/quasi-steady state.
• 3.  
  The boundary of the gas that stays gravitationally bound to the planet, according to Figure 4. The mass of this part of the gas is about 4.31 × 10²⁰ g.

Figure 5 demonstrates the different boundaries of the atmosphere given by these definitions with dashed lines in different colors. As we will see in Section 3.2, the gas enclosed by the 10^−8.5 g cm⁻³ contour line and the gas in the "egg-shape" region are of similar mass, both about half an order of magnitude higher than the bound mass. By subtracting the bound mass from the mass in the "egg-shape" region, we find the amount of recycling gas behind the planet to be 9.98 × 10²⁰ g.

Zoom In Zoom Out Reset image size

We would like to point out that the atmospheric pressure on the planet surface is probably a better metric of how much gas is retained, because in the current context the calculated atmosphere mass is relatively subjective—it depends on which boundary is used. In comparison, the surface pressure does not rely on a specific definition. The larger the surface pressure, the thicker the atmosphere. In cases where the planet surface is partially or completely molten, the surface pressure directly determines how much volatiles can be dissolved in the magma ocean (e.g., Wu et al. 2018; Sharp & Olson 2019). The average surface pressure for the canonical case is 4 × 10⁻² bar. We use both surface pressure and the atmosphere mass to measure the amount of proto-atmosphere throughout the study.

Does the proto-atmosphere of an eccentric planet always have a similar flow structure to that seen in the canonical simulation? How does it change in response to different planet–disk environments? To answer these questions, we vary the three input parameters—the relative velocity between the planet and the gas (${v}_{\infty }$), the planet mass (M_pl), and the background gas density (${\rho }_{\infty }$)—separately and carry out the corresponding "snapshot" simulations (see Table 1).

3.2.1. The Effect of Relative Velocity

The relative velocity between the planet and the disk gas directly determines the strength of the bow shock around the planet when it exceeds the sound speed (∼1.2 km s⁻¹ for the assumed background disk gas). In the parameter study, we vary the relative velocity between 3.2 and 7.5 km s⁻¹, corresponding to planets with orbital eccentricities between 0.1 and 0.5.

The hydrodynamic solutions of these simulations reveal that the flow structures around the planet are similar to those in the canonical case. The "egg-shape" region forms behind the planet with part of the gas bound to the planet and the rest exchanging with the disk, regardless of how fast the planet is traveling through the gas (Figure 6). As the relative velocity gets larger, the opening angle of the bow shock becomes narrower. The larger ram pressure coming from the nebular gas suppresses the accretion of the atmosphere further and reduces the stand-off distance of the shock front. As a result, the planet keeps a smaller atmosphere with an even smaller portion of the bound part. Figure 7 shows how the atmosphere mass and planet surface pressure change with the variation of the relative velocity. The atmospheric properties show most dramatic changes when the planet is traveling relatively slowly in the disk gas (but is still supersonic). In the cases of high relative velocity, the "egg-shape" region is dominated by the recycling movements of gas.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. The Effect of Planet Mass

The planet mass directly affects the amount of atmosphere that can be gravitationally accreted onto the planet. We vary the planet mass from 0.1 to 3 M_⊕, covering a range of terrestrial planets from the sub-Earth to super-Earth regime (Table 1). We assume that all simulated planets have Earth-like composition and interior structures (e.g., similar element ratios and temperature structures). Therefore, we calculate the planet radius as a function of planet mass using the ExoPlex mass–radius software package (Unterborn et al. 2018) for planets ≤0.7 M_⊕ and the mass–radius relationship described in Dressing et al. (2015) for planets >0.7 M_⊕.

Figure 8 shows the hydrodynamic solutions for planets with 0.1, 0.7, and 2.0 M_⊕. We find that except for the case where M_pl = 0.1 M_⊕ (e.g., the case of Mars), all other simulated systems maintain similar flow structures—a stable atmosphere with an "egg-shape" region. In the given context, a Mars-sized planet or planetary embryo is not able to hold a regular proto-atmosphere with the little gravity it provides. As can be seen from the leftmost panel in Figure 8, the supersonic gas forms a "secondary shock front" on the planet surface. The planet's gravity is not high enough to focus the incoming gas, i.e., no bow shock is formed, and hence the incoming gas is able to directly strike onto the planet's surface where it is shocked. The yellow dashed lines in the figure mark where the gas velocity reduces below sound speed. A planetary embryo with such a small size cannot hold a bound atmosphere. Recycling vortices are still present behind the planet, though this region can no longer enclose the whole planet. Note that the situation can change with different background setups, e.g., reducing the relative velocity between the planet and the disk gas could help a small planet object such as Mars to retain a normal atmosphere. If Mars was once scattered on an eccentric orbit (Hansen 2009), it is likely to experience dramatic atmosphere stripping and re-accretion repeatedly throughout the orbit. This can result in the loss of any outgassed atmosphere, which is replaced by the reducing H₂/He atmosphere.

Zoom In Zoom Out Reset image size

On the other hand, more massive planets are able to greatly slow down the supersonic gas and hold atmospheres of a larger mass. The "egg-shape" regions around these planets also reveal more complicated vortices and other flow structures.

Figure 9 summaries the trends of atmosphere mass and averaged planet surface pressure with increasing planet mass, similar to Figure 7. Both properties increase as the planet gets more massive, most dramatically in the sub-Earth mass range. In our setups, planets with sufficient orbital eccentricity do not possess a bound atmosphere when their mass is lower than 0.5 M_⊕. The hydrodynamics around planets with mass <0.7 M_⊕ are dominated by the recycling gas flows.

Zoom In Zoom Out Reset image size

3.2.3. The Effect of Background Gas Density

The background gas density can have values different from 10⁻⁹ g cm⁻³ for a number of reasons. As was mentioned in Section 2.2.2, different models for the solar/stellar nebula would provide augmented nebular density profiles (0.5–1.5 orders of magnitude higher). Besides the uncertainties from theoretical models, the age of the protoplanetary disk can also affect the nebular density in the background. Within 1–10 Myr, the nebular gas density can be lowered by one to two orders of magnitude as the disk evolves through accretion, wind, and photoevaporation (e.g., Mordasini et al. 2012; Bai 2016; Nakatani et al. 2018a, 2018b).

Unlike the relative velocity and the planet mass, the impact of the nebular gas density on our simulation solutions seems less straightforward. In our study, we perform the simulations with the background gas density varying from 10⁻⁷ to 10⁻¹⁰ g cm⁻³ while keeping other parameters the same. The resulting flow patterns and atmospheric properties do not change monotonically with the nebular density. When we lower the background gas density, the hydrodynamic solution presents the same flow pattern—the "egg-shape" recycling region (the right panel of Figure 10). However, the region expands and the planet is able to hold a more massive atmosphere. The decrease of the nebular density leads to a lower gas pressure in the background. With less ram pressure pushing against the moving planet, the bow shock becomes weaker and facilitates gas accretion onto the planet. In practice, we decrease the background gas density to as low as 10⁻¹¹ g cm⁻³, in which case the planet retains an even larger atmosphere. However, we exclude the results from this simulation because it converges too slowly and it becomes too expensive computationally to reach a quasi-steady state.

Zoom In Zoom Out Reset image size

When we elevate the background density, the hydrodynamics of the system seem to switch to a different regime. The solutions present a different flow structure (the left and middle panels in Figure 10). We no longer observe the stable recycling of gas through the "egg-shape" region. Instead, small-scale vortices constantly emerge and disintegrate behind the planet. In the case of 10⁻⁸ g cm⁻³ (DEN-8), we observe the presence of Kelvin–Helmholtz instability along the direction parallel to the shearing flow. The instability invokes propagating density waves, which are most prominent at two to eight planet radii behind the planet and weakened further downstream. In the case of 10⁻⁷ g cm⁻³ (DEN-7), some portion of gas appears to be flowing back to the planet near the symmetry axis. The formation of large vortices to initiate the recycling is unstable, however—different scales of vortices are frequently shredded and re-emerge. In both cases, no nebular gas is being stably recycled but a portion of atmosphere stays bound to the planet. Nonetheless, we can still define an atmosphere using a suitable density contour for each case⁹ . As the simulation converges, the atmosphere mass fluctuates around some constant value. By averaging these fluctuating values we have an estimate of the amount of atmosphere around the planet.

The fact that the hydrodynamics change with the background density is an effect of radiative transfer. Because the equations of hydrodynamics do not set any absolute scale of density, the solutions of simulations without radiative transfer should be self-similar, with the background density being a scaling factor. A typical optical depth for the proto-atmosphere within the bow shock is ∼200 when the background density is 10⁻⁹ g cm⁻³. This quantity reaches ∼600 when the background density is 10⁻¹⁰ g cm⁻³ and drops to ∼20–30 when the background density is 10⁻⁸ or 10⁻⁷ g cm⁻³. The drastic drops in opacity and optical depths of atmospheres in the cases of high background density are caused by the evaporation of dust—the gas temperatures near the planet surface exceed 1500 K. This leads to a different hydrodynamics regime to the recycling regime in the cases of lower background density.

Figure 11 summarizes how the atmosphere mass and surface density change with the background gas density. We observe two hydrodynamics regimes operating in cases of different background density, which are likely invoked by the opacity jump in the atmospheres. In either regime, there are two competing mechanisms that affect how much atmospheres the planet can hold. In the regime where instabilities are present and the large-scale gas recycling is absent (DEN-7 and DEN-8), the atmosphere is more massive because of the sufficient gas supply from a dense nebular environment. This factor wins over the stripping effects from a larger ram pressure against the bow shock. The atmosphere is thin and compact. In the recycling regime, however, the drop in background gas density does not significantly limit the availability of nebular gas for accretion, but rather facilitates gas accretion by weakening the ram pressure. The atmosphere is thick and puffy. We conclude that the recycling regime is favorable for planets to keep a stable atmosphere in the planetary bow shock.

Zoom In Zoom Out Reset image size

Lastly, we would like to point out that the parameter study in background gas density also provides qualitative estimates for the cases where planets are located at different distances from the star. If a planet is on an orbit further away from the star (e.g., 10 au), the decrease in the surrounding gas temperature and density can result in a reduction in the ram pressure against the planetary bow shock. The flow structure of the proto-atmosphere likely resembles that of the low-density case study, where the planet retains a puffy and massive atmosphere in the recycling regime. On the other hand, a planet located closer to the star (e.g., 0.1 au) is exposed to a denser and hotter nebular environment. With the surrounding temperature approaching 1000 K and the nebular density raised by 2–3 orders of magnitude, the proto-atmosphere likely presents instabilities as in those high-density case studies.

In the second part of our hydrodynamic simulations, we focus on tracking the evolution of the proto-atmosphere as the planet moves on an eccentric orbit. All parameters relevant to the planet–disk environment are time-dependent, including the background gas density, the temperature, and the relative velocity between the planet and the gas (only the speed is time-dependent, the direction of velocity is assumed to be fixed given that the change of bow shock orientation is negligible in local-framed simulations—Figure 2). Therefore, the outcomes of the simulation are solely controlled by the orbital parameters (the semimajor axis and the eccentricity) as well as the planet mass.

We model a planet of 1 M_⊕ traveling on an orbit with semimajor axis a = 1 au and eccentricity e = 0.2. The simulation starts with the planet at its perihelion. Figure 12 presents the periodic features reflected in the atmospheric properties as the planet completes 1.9 orbits. At small phase angles (<π/3), the planet undergoes a fast-accretion stage, trying to reach an equilibrium with the surrounding nebula. To ensure that the fluctuations of all calculated quantities are caused merely by orbital effects, we exclude the first half of an orbital period. From the first aphelion to the second aphelion (phase angle = π–3π, shaded by light gray in Figure 12), the planet completes a full orbit.

Zoom In Zoom Out Reset image size

We can see that the relative velocity between the planet and the gas rises from a minimum as the planet leaves its aphelion. However, the atmospheric properties—the surface pressure, temperature, the mass, and the mass accretion/loss rate—do not respond to the change immediately. The proto-atmosphere should have been undergoing stripping with the increasing relative velocity/Mach number, but the surface pressure, temperature, and mass continue to grow for a while before decreasing. On the other hand, the atmosphere is still experiencing mass loss even when the relative velocity starts to decline after reaching a peak (phase angle = 3π/4). To sum up, the system shows a delay in response to the time-dependent boundary conditions. The oscillations in atmospheric properties show a phase shift of ∼π/4, corresponding to an average reaction time¹⁰ ∼7 × 10⁶ s, which matches the dynamical timescale of the system.

Note that in the orbit-dependent simulation, the relative velocity is not the only parameter affecting how much of the atmosphere can be retained. The time-dependent background gas temperature, density, and pressure can also make a difference, although their variations within one orbit are much smaller than those of the relative velocity. For example, the relative velocity/Mach numbers are similarly small for the planet at aphelion and perihelion. But the planet seems to retain more atmosphere after it passes aphelion, because the denser and hotter disk environment at perihelion leads to a larger ram pressure toward the planet, making the gas accretion less effective. We also notice that the changes in the surface temperature with phase are not entirely synchronized with the other properties—it drops much faster than it rises. Spontaneously, one could expect the surface temperature to change in phase with the change of surface pressure for a pressure-supported atmosphere. But when the surface pressure is high, the surface temperature can reach 1500–1600 K and drastically lower the local opacity (dust is evaporated). So as the atmosphere gets stripped (at high relative velocity/Mach number) it undergoes fast radiative cooling. However, the high relative velocity also produces a strong bow shock that dramatically heats up the compressed gas at the shock front. The heat can be transferred to the atmosphere from outside to inside, creating a temporary temperature inversion. Eventually, the heat from shock gas prevents the surface temperature from further decrease and causes it to rise again (Figure 13).

Zoom In Zoom Out Reset image size

The orbit-averaged values of different quantities—drawn as horizontal solid lines in Figure 12—are consistent with the results acquired from the "snapshot" simulation VEL4.2 (horizontal dotted–dashed lines). With the appropriate selection of input parameters, the "snapshot" simulation is also able to represent the "averaged" picture of the dynamic proto-atmosphere over the corresponding orbit.

The proto-atmosphere on low-mass planets with circular orbits has been studied in detail to explore the formation of super-Earths and Earth-like planets (see literature review in Section 1). By comparing our results with the outcomes from hydrodynamic simulations for planets on circular orbits, we identify a few major differences in atmosphere accretion on planets with the two types of orbits.

(1) The flow structure. The relative motion between the planet and gas is subsonic in cases where planets are on circular orbits. The gas accretion does not generate a bow shock. In 1D hydrodynamic simulations (e.g., Stökl et al. 2015, 2016), proto-atmosphere accretion on the planet is assumed to be in azimuthal symmetry—the infall of the nebular gas lies in the radial direction. In 2D and 3D simulations (e.g., Fung et al. 2015; Ormel et al. 2015a, 2015b; Cimerman et al. 2017; Kuwahara et al. 2019), the planet is considered to be embedded in a Keplerian disk. A local shearing flow with its reference frame is comoving with the planet as the gas travels faster on the side facing the star and slower on the other side. Both 2D and 3D simulations find a pressure-supported rotating atmosphere around the planet. In 2D, the atmosphere becomes isolated from the disk gas when the system reaches a steady state. In 3D, the atmosphere stays connected to the disk as the atmosphere is constantly exchanging gas with the nebula. The incoming gas enters the Bondi sphere around the planet at high latitudes and leaves at the equatorial regions. In all simulations, two horseshoe flow regions are seen on the front and rear sides of the planet. The situation of a planet experiencing a headwind from the disk gas in sub-Keplerian motion has also been considered. In such cases, the flow pattern becomes asymmetric around the planet as the streamlines are compressed toward the headwind. The nebular gas enters from the equatorial region and leaves at high latitudes. For figures of the detailed flow patterns of the above simulations, please refer to Fung et al. (2015), Ormel et al. (2015a, 2015b), Cimerman et al. (2017), and Kuwahara et al. (2019). In summary, proto-atmospheres on planets with both circular (subsonic) and eccentric (supersonic) orbits are actively interacting with the nebular gas through large-scale recycling. The recycling/replenishment timescale in the case of a circular orbit is generally longer (>10⁷ s) than that of a typical eccentric planet simulated in this study (5 × 10⁵–5 × 10⁶ s).

(2) The steady state. In 2D simulations for planets on circular orbits, the gas accretion seems to be able to reach a steady state where the solutions no longer fluctuate. The flow pattern and other physical quantities (e.g., density, velocity), however, could be time-variable in some 3D simulations (Cimerman et al. 2017). But some other simulations suggest that they can reach steady states (T. Moldenhauer et al. 2020, in preparation). In comparison, gas accretion onto planets with eccentric orbits can also reach steady/quasi-steady states as shown in this study.

(3) Atmospheric properties. Planets on circular orbits are able to hold much more massive atmospheres. The 1D models in Stökl et al. (2015, 2016) point out that a 1 M_⊕ planet can accrete up to 300 bar (or 10²⁵ g) of atmosphere, which is ∼0.2% of the planet mass. The 2D simulations in Ormel et al. (2015a) and 3D models (T. Moldenhauer et al. 2020, in preparation) give similar or slightly lower values. In terms of surface temperature, a 1 M_⊕ planet on a circular orbit can be ≥3000 K (T. Moldenhauer et al. 2020, in preparation; Stökl et al. 2016), whereas for an eccentric orbit we find 1000–2000 K. This is due to the accretion of a much smaller atmosphere and implies that dust grains dominate the opacity in such an atmosphere.

A summary of the above comparison is presented in Table 2.

The accretion on a gravitating object traveling supersonically in a homogeneous medium is a classical problem called the BHL accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944). For decades, numerous analytical and numerical studies have been exploring the flow structure and the stability of BHL accretion (see Edgar 2004; Foglizzo et al. 2005 for reviews). Typical applications for these BHL accretion models include wind accretion in binary systems (e.g., Taam et al. 1991; Xu & Stone 2019), accretion on other compact objects such as black holes (e.g., Pogorelov et al. 2000), and the protostellar clusters and galaxy clusters (e.g., Sakelliou 2000; Bonnell et al. 2001).

Gas accretion on an eccentric planet is also a BHL accretion problem in principle, although the complications from radiative properties and boundary conditions can deviate the solutions from the classical derivatives. In fact, the flow structure and stability of BHL accretion are proven to depend strongly on the chosen grid method (2D planar, 2D axisymmetric, or 3D), the given parameters, and the accretor's properties (Foglizzo et al. 2005). Our results in this study are aligned with the previous finding that simulations in a 2D axisymmetric case are in general stable (e.g., Pogorelov et al. 2000; Foglizzo et al. 2005). The "egg-shape" flow pattern—where a portion of the gas flows back to the planet in the post-shock region—is also commonly observed in other axisymmetric simulations (e.g., Koide et al. 1991; Pogorelov et al. 2000; MacLeod & Ramirez-Ruiz 2015).

Though planets and planetary embryos can be scattered and/or migrate in the young protoplanetary disk, their phases on eccentric orbits are believed to be temporary as long as the disk gas is present. Planets traveling in a gaseous disk can experience both aerodynamic and gravitational gas drag. The aerodynamic drag is the pressure resistance that any gas or fluid exerts on the object moving inside it, and it is more effective for smaller bodies such as planetesimals (e.g., Kominami & Ida 2002). The gravitational drag arises from torques as the planet interacts with the gaseous disk (e.g., Papaloizou & Larwood 2000). It is more prominent in massive objects such as planetary embryos and planets.

The aerodynamic damping timescale can be calculated from the formulae in Tanaka & Ida (1999) and Kominami & Ida (2002) as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\tau }_{{\rm{a}}{\rm{e}}{\rm{r}}{\rm{o}}} & \approx & \displaystyle \frac{{M}_{{\rm{p}}{\rm{l}}}}{\pi {{\boldsymbol{ \mathcal R }}}_{{\rm{p}}{\rm{l}}}^{{\bf{2}}}{\rho }_{\infty }{v}_{\infty }}\approx 3\times {10}^{5}{\left(\displaystyle \frac{e}{0.2}\right)}^{-1}{\left(\displaystyle \frac{{M}_{{\rm{p}}{\rm{l}}}}{{M}_{\oplus }}\right)}^{1/3}\\ & & \times \,{\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-9}{\rm{g}}{{\rm{c}}{\rm{m}}}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{r}{1{\rm{a}}{\rm{u}}}\right)}^{1/2}\,{\rm{y}}{\rm{r}},\end{array}\end{array}\end{eqnarray} \tag{ 10 }$$

where we approximate the relative velocity between the planet and the gas ${v}_{\infty }$ as (1/2)ev_K, and v_K = (GM_⊙/r)^1/2 is the Keplerian speed.

For low eccentricities (e.g., e ≤ H/r or ${v}_{\infty }\leqslant {c}_{\infty }$), the damping timescale for the gravitational drag is given as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\tau }_{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{v}},{\rm{l}}{\rm{e}}} & \approx & \left(\displaystyle \frac{{M}_{\odot }}{{M}_{{\rm{p}}{\rm{l}}}}\right)\left(\displaystyle \frac{{M}_{\odot }}{{{\rm{\Sigma }}}_{\infty }{r}^{2}}\right){\left(\displaystyle \frac{{c}_{\infty }}{{v}_{{\rm{K}}}}\right)}^{4}{{\rm{\Omega }}}_{{\rm{K}}}^{-1}\\ & \approx & 5\times {10}^{2}{\left(\displaystyle \frac{{M}_{{\rm{p}}{\rm{l}}}}{{M}_{\oplus }}\right)}^{-1}{\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-9}{\rm{g}}{{\rm{c}}{\rm{m}}}^{-3}}\right)}^{-1}\\ & & \times \,{\left(\displaystyle \frac{r}{1{\rm{a}}{\rm{u}}}\right)}^{-3/4}\,{\rm{y}}{\rm{r}}.\end{array}\end{array}\end{eqnarray} \tag{ 11 }$$

${{\rm{\Sigma }}}_{\infty }={\rho }_{\infty }H\sqrt{2\pi }$ is the surface density of the disk gas. Ω_K = (GM_⊙/r³)^1/2 is the Keplerian frequency (Kominami & Ida 2002).

However, Papaloizou & Larwood (2000) found that when e > 1.1H/r, the sign of the torque on the planet reverses and leads to an eccentricity damping timescale strongly dependent on e. The ratio of the disk scale height to its radius, H/r, ranges from 0.02 to 0.07 throughout the protoplanetary disk (0.1–10 au). For planets in the supersonic regime as we consider in this study (e ≥ 0.1), we need to calculate a damping timescale different from Equation (11). For such cases, the authors suggest to replace τ_grav,le by [e/(H/r)]³τ_grav,le. So we have

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{grav},\mathrm{he}} & \approx & 1\times {10}^{5}{\left(\displaystyle \frac{e}{0.2}\right)}^{3}{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{{M}_{\oplus }}\right)}^{-1}{\left(\displaystyle \frac{{\rho }_{\infty }}{{10}^{-9}{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1}\\ & & \times \,{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-3/2}\mathrm{yr}.\end{array}\end{eqnarray} \tag{ 12 }$$

The ratios of the these timescales are

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{aero}}}{{\tau }_{\mathrm{grav},\mathrm{le}}}\approx 600\times {\left(\displaystyle \frac{e}{0.2}\right)}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{{M}_{\oplus }}\right)}^{4/3}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{5/4}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{aero}}}{{\tau }_{\mathrm{grav},\mathrm{he}}}\approx 3\times {\left(\displaystyle \frac{e}{0.2}\right)}^{-4}{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{{M}_{\oplus }}\right)}^{4/3}{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{2}.\end{eqnarray} \tag{ 14 }$$

Therefore, for planets more massive than 10²⁵ g on orbits with low eccentricities, the orbits are mainly circularized by the gravitational drag (Kominami & Ida 2002). For planets on orbits with higher eccentricities, aerodynamic gas drag dominates if M_pl < 0.1 M_⊕, while gravitational drag dominates if M_pl > 1 M_⊕. These calculations are consistent with calculations of eccentricity damping by Morris et al. (2012).

The orbital eccentricity can also be damped by the dynamical friction coming from the swarms of planetesimals scattered throughout the disk. The damping timescale has a similar expression to Equation (11), except that ${{\rm{\Sigma }}}_{\infty }$ is substituted by the surface density of solids in the disk Σ_solid and ${c}_{\infty }$ is substituted by the velocity dispersion of planetesimals v_dis. Therefore, the mechanism can either be incorporated into the gravitational drag regime (Kominami & Ida 2002) or neglected in cases where e > 0.03 because the gravitational drag dominates over dynamical friction (Morris et al. 2012).

The damping of orbital eccentricity follows

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{{dt}}=-e{\tau }_{\mathrm{damp}}^{-1}\end{eqnarray} \tag{ 15 }$$

where τ_damp is the damping timescale. If it is substituted by τ_grav,he, de/dt becomes linearly proportional to e⁻² and M_pl. The lower the eccentricity (still larger than 0.1) and the larger the planet mass, the faster its orbit gets damped.

For a typical planet simulated in our models, its eccentric orbit is circularized in about 0.1 Myr, which is much longer than the relevant timescales in this study. The damping of eccentricity is therefore negligible in both "snapshot" and orbit-dependent simulations. The planets may be able to hold more atmospheres as their orbits get circularized if the ambient disk gas has not yet dissipated. But this procedure can be complicated, because the damping also vanishes with the disappearing disk gas.

The study does not consider the cases where the relative velocity between the planet and the gas is below sound speed (e.g., orbital eccentricity is very small, e < 0.08). In such cases, the planet does not generate a bow shock that compresses and heats up the gas in front of the planet, and the stripping of the proto-atmosphere will be much less effective. In fact, it is no longer suitable to model planets in this regime with 2D axisymmetric simulations because the relative velocity (${v}_{\infty }$) could be comparable to the shear or headwind velocity of the ambient gas. To model gas accretion on planets traveling at subsonic speed, we suggest using the grid methods adopted for planets experiencing a headwind on circular orbits (Ormel et al. 2015a, 2015b, T. Moldenhauer et al. 2020, in preparation). The proto-atmosphere is predicted to be asymmetric and less massive than in the circular cases, but much more massive than in the supersonic cases. As our paper solely focuses on eccentric planets in the supersonic regime, the modeling of gas accretion on subsonic planets is left as future work. It will be interesting to compare the results from these models with the circular cases and the supersonic cases.

Another limitation of this study comes from the assumption of axisymmetry. A number of numerical simulations on BHL accretion have adopted the full 3D grid method to represent the more realistic geometry. Unfortunately, most of them are not stable (Foglizzo et al. 2005). The 3D simulations usually break the axial symmetry of the flow structure and develop instabilities downstream of the bow shock (e.g., Ruffert 1994; Ruffert & Arnett 1994). These instabilities are not as prominent and strong as the "flip-flop" behaviors observed in 2D planar simulations and could arise from numerical artifacts. Nonetheless, the gas is still observed to be flowing back toward the accretor in the post-shock region, similar to the 2D axisymmetric simulations. We expect that the hydrodynamic solutions for the atmospheres (e.g., density, pressure, temperature, mass) will be within the same orders of magnitude as obtained in this study, regardless of the adopted grid methods.

Lastly, we would like to point out that we neglect solid accretion that could occur at the same time as gas accretion. We assume that the planet never grows in mass and the planet luminosity is barely affected by pebble/planetesimal accretion. Solid accretion is common and important for planet formation in young disks (≤5 Myr), and luminosity feedback could be important (Mordasini 2013). There are several studies that have considered the impacts of orbital eccentricity on solid accretion (e.g., Lissauer et al. 1997; Liu & Ormel 2018). Given that the timescales relevant to gas accretion in this study are days to weeks, the neglect of long-term effects from solid accretion can be justified.