missing

# Enhanced Self-Organizing Map Solution
for the Traveling Salesman Problem

Joao P. A. Dantas1, Andre N. Costa1, Marcos R. O. A. Maximo2, Takashi
Yoneyama3

###### Abstract

Using an enhanced Self-Organizing Map method, we provided suboptimal solutions
to the Traveling Salesman Problem. Besides, we employed hyperparameter tuning
to identify the most critical features in the algorithm. All improvements in
the benchmark work brought consistent results and may inspire future efforts
to improve this algorithm and apply it to different problems.

##  1 Introduction

The Traveling Salesman Problem (TSP) is one of the most studied routing
problems within the combinatorial optimization field [Laporte 1992], which is
defined as a salesman that must visit a series of cities once and only once,
returning to the starting city [Brocki and Koržinek 2007]. Although coming
from a relatively simple idea, finding the solution for the problem, i.e., the
shortest path, is NP-Hard [Huang et al. 2017], thence the great interest to
find efficient ways to solve it.

Due to the problem’s complexity and to the fact that many of its applications
require fast ways to solve it, it is common to employ heuristics to generate
approximate (suboptimal) solutions to the TSP. Examples of these methods are
Genetic Algorithm (GA), Simulated Annealing (SA), tabu search, and ant colony
optimization [Xu et al. 2008, Calado and Ladeira 2011]. Another form to find
approximate solutions to complex optimization problems is to use neural
networks [Peterson 1990], which can be efficient due to their adaptability.

According to the literature, the first neural network solution for the TSP was
the one presented in [Hopfield and Tank 1985], based on the minimization of an
energy function. From that, many methods have been proposed to improve the
Hopfield neural network [Xu et al. 2008]. Another use of neural networks for
solving the TSP is the one proposed in [Kohonen 1998] through competitive
unsupervised learning based on winner-take-all and winner-take-most
algorithms, also called Self-Organizing Map (SOM), due to how the adaptation
of the neurons works [Brocki and Koržinek 2007]. The SOM adjusts its neurons
to fit the input cooperatively by inspecting the list of cities. This
localized response to the input list generates a neighborhood preserving map,
resulting in a near-optimal path [Xu et al. 2008].

In summary, most of the TSP heuristics solutions in the literature have used
one of the following methods: Hopfield network, Kohonen’s SOM, GA, and SA
[Markovic et al. 2012]. In comparison with other heuristic methods, the SOM
has presented low computation complexity and promising performance. Therefore,
these neural networks have attracted researchers to explore and enhance their
performance when applied to the TSP [Favata and Walker 1991, Budinich 1996].

In our work, the main contribution is to modify the method of SOM as proposed
and implemented in [Martin 2018], adding some features and tuning
hyperparameters. We offer these modifications to solve a particular problem
involving TSPs varying in size from $50$ to $200$. The remainder of this paper
is structured as follows: Section 2 provides details with respect to the
experiment setting; Section 3 states the improvements on the method for this
context; Section 4 presents the results obtained; Section 5 brings the
conclusions and the possible future work.

##  2 Experiment

The proposed experiment aimed at predicting edges belonging to optimal
solutions of $2{,}000$ two-dimensional Euclidean TSP ranging from $50$ to
$200$ city nodes. We evaluate the solution considering the F1 score of the
predicted adjacency matrix compared to the optimal solution. The F1 score is
the harmonic mean of the precision and recall, where an F1 score reaches its
best value at $1$ and worst score at $0$ [Zhang et al. 2015]. The relative
contribution of precision and recall to the F1 score is equal, which is
preferred over accuracy as this problem is often very imbalanced [Mele et al.
2021].

##  3 Improvements

In this section, we discuss the main improvements performed in the method
adopted. We employ hyperparameter tuning to adapt the algorithm. Besides, we
created two additional features to find solutions with better results.

###  3.1 Hyperparameter Tuning

When exploring the algorithm used in this work, we first identified its
primary hyperparameters (population size, number of iterations, learning rate,
and discount rates for the latter two) and their effects on the final scores.
To understand the influence of each of these factors and find their best
configuration, we employed a single-factor design to tune the proposed
technique, i.e., varying approximately one feature at a time to examine its
isolated influence. The search for the hyperparameters that better suit the
dataset led us to the baseline algorithm, which received all modifications
discussed next.

###  3.2 Additional Modifications

The first improvement to the baseline algorithm was to change the way to
choose the first node considered. The standard technique randomly selects the
starting point, which may lead the solution into a local minimum, depending on
the location of this first node. In addition to the randomly chosen initial
city, we forced the algorithm to start from the city in the centermost
position and the furthest position from the centroid of all cities. After
running the algorithm in these three different initial positions, we chose the
one that presents the shortest path.

As the second modification, after employing the hyperparameters tuning and
identifying the most significant feature as the population size, we improved
the algorithm based on the variation of this hyperparameter in each SOM
iteration. Figure 1 exemplifies the SOM iterations until the result with the
shortest path is obtained. We used 20 different population sizes in each TSP,
from 1 to 20 times the number of cities, and we calculated the route length of
all of these TSP solutions, choosing the one with the shortest path. Indeed,
that additional feature turned the process computationally more costly.
However, it brought better results since the population size was constantly
changing, and we chose the one that produced the shortest distance in all
solutions.

Figure 1: SOM algorithm iteration samples.

##  4 Results

Table 1 shows the F1 Scores results of the hyperparameters tuning and
clarifies the importance of the population size feature. We chose as our
baseline the following configuration: $100{,}000$ iterations, $0.9997$ for the
discount rate of initial neighborhood, $0.8$ for the learning rate, $0.99997$
for the discount rate of the learning rate, and $6$ for the population size
multiplier factor.

Table 1: Hyperparameters tuning with F1 scores

Number of | Discount Rate of | Learning | Discount Rate of | Population Size | F1
---|---|---|---|---|---
Iterations | Initial Neighborhood | Rate | Learning Rate | Multiplier Factor | Score
100 | 0.9997 | 0.8 | 0.99997 | 8 | 0.06878
100,000 | 0.9997 | 0.8 | 0.99997 | 8 | 0.07800
100,000 | 0.9997 | 0.8 | 0.99997 | 16 | 0.07784
100,000 | 0.9997 | 0.8 | 0.99997 | 4 | 0.07855
100,000 | 0.9997 | 0.8 | 0.99997 | 2 | 0.07853
100,000 | 0.9997 | 0.8 | 0.99997 | 6 | 0.07885
100,000 | 0.9997 | 0.8 | 0.99997 | 5 | 0.07869
100,000 | 0.99997 | 0.8 | 0.99997 | 6 | 0.07854
100,000 | 0.997 | 0.8 | 0.99997 | 6 | 0.07834
100,000 | 0.997 | 0.01 | 1.00000 | 6 | 0.07769
100,000 | 0.9997 | 0.8 | 0.99997 | 7 | 0.07878

After choosing the baseline hyperparameter configuration, we employed the
additional modifications described in subsection 3.2, achieving slightly
better results with $0.07891$ of the evaluation metric proposed. Therefore,
these results show that, through subtle and well-planned changes in algorithms
already consolidated by the scientific community, it is possible to obtain
better results in classical optimization problems.

##  5 Conclusions

This work brings an enhanced way to use SOM to find suboptimal solutions for
the TSP. Firstly, we employ hyperparameter tuning to identify the baseline
algorithm configuration. Secondly, we combined two different modifications to
the benchmark technique, leading to consistent route length results. For
future work, we suggest additional improvements, mainly considering the
stochastic aspects of the method, and advanced techniques to choose the best
hyperparameters, which may bring better results to the SOM algorithm.

## References

  * [Brocki and Koržinek 2007] Brocki, Ł. and Koržinek, D. (2007).  Kohonen self-organizing map for the traveling salesperson problem.  In Recent Advances in Mechatronics, pages 116–119. Springer.
  * [Budinich 1996] Budinich, M. (1996).  A self-organizing neural network for the traveling salesman problem that is competitive with simulated annealing.  Neural Computation, 8(2):416–424.
  * [Calado and Ladeira 2011] Calado, F. M. and Ladeira, A. P. (2011).  Problema do caixeiro viajante: Um estudo comparativo de técnicas de inteligência artificial.  e-xacta, 4(1).
  * [Favata and Walker 1991] Favata, F. and Walker, R. (1991).  A study of the application of kohonen-type neural networks to the travelling salesman problem.  Biological Cybernetics, 64(6):463–468.
  * [Hopfield and Tank 1985] Hopfield, J. and Tank, D. (1985).  Neural computation solutions for tsp applications.  Biological Cybernetics, 52:141–152.
  * [Huang et al. 2017] Huang, L., Wang, G.-c., Bai, T., and Wang, Z. (2017).  An improved fruit fly optimization algorithm for solving traveling salesman problem.  Frontiers of Information Technology & Electronic Engineering, 18(10):1525–1533.
  * [Kohonen 1998] Kohonen, T. (1998).  The self-organizing map.  Neurocomputing, 21(1-3):1–6.
  * [Laporte 1992] Laporte, G. (1992).  The traveling salesman problem: An overview of exact and approximate algorithms.  European Journal of Operational Research, 59(2):231–247.
  * [Markovic et al. 2012] Markovic, D., Madic, M., Tomic, V., and Stojkovic, S. (2012).  Solving travelling salesman problem by use of kohonen self-organizing maps.  Acta Technica Corviniensis-Bulletin of Engineering, 5(4):21.
  * [Martin 2018] Martin, D. (2018).  Solving the traveling salesman problem using self-organizing maps.  https://github.com/DiegoVicen/som-tsp.
  * [Mele et al. 2021] Mele, U. J., Gambardella, L. M., and Montemanni, R. (2021).  A new constructive heuristic driven by machine learning for the traveling salesman problem.  Algorithms, 14(9):267.
  * [Peterson 1990] Peterson, C. (1990).  Parallel distributed approaches to combinatorial optimization: benchmark studies on traveling salesman problem.  Neural computation, 2(3):261–269.
  * [Xu et al. 2008] Xu, X., Jia, Z., Ma, J., and Wang, J. (2008).  A self-organizing map algorithm for the traveling salesman problem.  In 2008 Fourth International Conference on Natural Computation, volume 3, pages 431–435. IEEE.
  * [Zhang et al. 2015] Zhang, D., Wang, J., and Zhao, X. (2015).  Estimating the uncertainty of average f1 scores.  In Proceedings of the 2015 International Conference on The Theory of Information Retrieval, pages 317–320.


# Growing Planet Envelopes in Spite of Recycling Flows

Avery P. Bailey Department of Physics and Astronomy, University of Nevada, Las
Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154-4002, USA  Nevada
Center for Astrophysics (NCfA), University of Nevada, Las Vegas, NV, USA
Zhaohuan Zhu Department of Physics and Astronomy, University of Nevada, Las
Vegas, 4505 South Maryland Parkway, Las Vegas, NV 89154-4002, USA  Nevada
Center for Astrophysics (NCfA), University of Nevada, Las Vegas, NV, USA

###### Abstract

The hydrodynamic exchange of a protoplanet’s envelope material with the
background protoplanetary disk has been proposed as one mechanism to account
for the diversity of observed planet envelopes which range in mass fractions
of $\sim 1\%$ for super-Earths to $\sim 90\%$ for giants. Here we present and
analyze 3D radiation-hydrodynamics models of protoplanet envelopes to
understand how the exchange of mass and energy between the protoplanet and
background disk influences the formation process. Our protoplanet envelope
simulations show an exchange of material bringing the outer $\gtrsim 0.4R_{b}$
envelope to steady state. This exchange provides a continuous source of
energy, which acts to increase the observed luminosity beyond that inferred
from the binding energy liberated from Kelvin-Helmholtz contraction alone – a
finding important for potential protoplanet observations. The inner $\lesssim
0.4R_{b}$, on the other hand, appears insulated – growing in accordance with
1D quasi-static theory. We incorporate these 3D hydrodynamic effects into an
extensible 1D framework with a physically motivated three-layer recycling
parameterization. Specializing to the case of Jupiter, recycling produces
minimal changes to the growth rate with the planet still entering runaway
accretion and becoming a gas giant in $\sim 1$ Myr. Even in the inner disk
($0.1$ AU), our 1D models suggest that recycling is not so robust and
ubiquitous as to stop all cores from becoming giants. At the same time
however, this recycling can delay a runaway phase by an order-of-magnitude
depending on the inner disk conditions and core mass.

planets and satellites: formation, planets and satellites: gaseous planets,
planet–disk interactions

##  1 Introduction

The core accretion theory of planet formation exists as a favorable mechanism
to form gas giants akin to Jupiter. In the core accretion model, gas accretion
and the consequent setting of planet formation timescale is mediated by the
Kelvin-Helmholtz contraction of the protoplanet envelope (Mizuno, 1980;
Pollack et al., 1996). Traditionally, this slow contraction has been modeled
as a 1D quasi-static process where the envelope structure is considered to be
nearly hydrostatic with small deviations arising from gravitational
contraction which drive the long term envelope evolution (Bodenheimer &
Pollack, 1986; Ikoma et al., 2000). Historically core accretion and its
subsequent extensions (e.g. pebble accretion) have been the de-facto standard
for solar system formation and planet formation in general. To its credit,
core accretion naturally explains the metal-rich interiors of Jupiter and
Saturn (Guillot et al., 1997; Wahl et al., 2017; Militzer et al., 2022) and
the correlation between giant planet occurrence and host star metallicity
(Fischer & Valenti, 2005). More recently, observational challenges and
theoretical considerations (driven by increased computing power and
advancement in the methods for numerical solution of the fluid equations) have
prompted a re-examination of the quasi-static core accretion framework.

From an observational point of view there appear at least three pieces of
evidence incompatible with the traditional picture of core accretion. First is
the existence of planets with core masses far exceeding what is expected from
the core accretion scenario. HD 149026b, for example, is a transiting Saturn
mass planet with unusually high density that internal models would estimate as
having a core of $\sim 70M_{\odot}$. While the existence of such a core argues
for a core accretion origin over something like gravitational instability, the
inferred core mass far exceeds the relatively robust 5-20 $M_{\oplus}$
critical core mass predicted by quasi-static core accretion models (Piso &
Youdin, 2014; Piso et al., 2015). Second is the existence of numerous super-
Earths and mini-Neptunes in exoplanet systems. These planets have atmosphere
to core mass ratios at the order of $1\%$, massive enough to suggest they
formed in the gaseous protoplanetary disk. At they same time, they obviously
did not runaway to become gas giants even though many have core masses larger
than the theoretical critical core mass. Finally is an observed abundance of
sub-Saturn mass planets compared to the predictions of core accretion. In
theory, there ought to be a bimodal distribution of planets consisting of a
population of gas giants that were massive enough to enter runaway accretion
and a low-mass population that failed to runaway before the disk dispersed
with the minimum occurrence rate landing in the sub-Saturn/Neptune regime.
While population synthesis models exhibit some indication of this bimodality
(Ida & Lin, 2004; Mordasini, 2018), sub-Saturns appear to be nearly as common
as their Jovian counterparts (Dong & Zhu, 2013; Petigura et al., 2018).

Several solutions have been put forth to address these observational tensions.
Late formation or high opacity would allow super-Earths to avoid the
transition to runaway accretion (Lee et al., 2014; Lee & Chiang, 2015). A
hierarchical merging formation channel where cores grow through giant impacts
(Ikoma & Hori, 2012; Ginzburg et al., 2016) or gap-opening (Fung & Lee, 2018;
Ginzburg & Chiang, 2019) would similarly slow gas accretion and ease these
tensions. The advective transport of energy and mass between the protoplanet
envelope and the protoplanetary disk, termed atmospheric recycling, has also
been advanced as a mechanism to explain the abundance of sub-Saturns (Lee,
2019) and super-Earths (Ormel et al., 2015; Cimerman et al., 2017; Ali-Dib et
al., 2020; Moldenhauer et al., 2021, 2022) and is the primary subject of this
paper.

From a theoretical perspective, some recycling must occur as all protoplanet
envelopes are embedded in the Keplerian shear of a background protoplanetary
disk at some point in their evolution. As an inherently multi-dimensional
phenomenon however, only recently, with advances in computing power and
techniques has it been possible to simulate this recycling interaction in 3D
with the appropriate physics. While much impressive work has been done on the
subject, the full scope of recycling and its application to the full diversity
of planet forming environments is not entirely understood. If super-Earths and
gas giants are thought to form in the same core accretion framework with
recycling stalling the super-Earth evolution but not the giants for example,
there must be some differential recycling efficiency or preferential behavior
surrounding the circumstances of formation. The magnitude of this efficiency
can be best understood from 3D radiation-hydrodynamical simulations which
themselves suffer from short simulation times ($\sim 10-10^{3}$ orbits)
compared to the evolutionary timescale. This necessitates careful methods and
analysis to properly interpret 3D models in the context of long-term evolution
and eventual formation outcome. In this work we examine and refine some of our
existing 3D radiation-hydrodynamics models of protoplanetary envelopes, taking
care to study the mass and energy budget of the envelope as these diagnostics
inform the long term evolution of these protoplanets.

In section 2 we present numerical methods and an overview of the dynamics
present in our radiation-hydrodynamics envelope models. Section 3 quantifies
the mass exchange between planet and background disk by examining the envelope
mass budget and tracking the circumplanetary flow topology with passive
scalars. Because hydrodynamic flows also transfer energy between envelope and
disk, we also examine the envelope energetics and present them in section 4.
In section 5 we summarize our findings and interpret these models in the
context of core accretion and the formation of planets like Jupiter along with
close-in super-Earths/mini-Neptunes.

##  2 Simulation Details

The simulations and setup employed here are similar to those of Bailey et al.
(2023) but generally at higher resolution, extended simulation time, and with
the inclusion of passive scalars. We summarize the most important information
here but additional details may be found in the preceding work. The radiation-
hydrodynamics simulations are performed with Athena++ in a local frame
centered on the planet with boundary conditions fixed to those of the
background disk. The mass and momentum equations for this system are the
standard inviscid hydrodynamics ones with rotational source terms arising from
Coriolis and centrifugal forces along with gravitational source terms from the
planet and central star. The result is a stratified shearing-box with a planet
potential (though with fixed boundaries rather than shear-periodic). Radiative
transfer is calculated with a formal solution of the transfer equation:

$\frac{dI}{d\tau}=S-I$ (1)

with a short characteristics method (Davis et al., 2012). Radiation is then
coupled to the hydrodynamics by computing $J=\int Id\Omega$ and adding as a
source term in the dimensionless energy equation:

$\frac{\partial E}{\partial t}+\nabla\cdot\left(Ev+pv\right)=-\rho
v\cdot\nabla\Phi+\rho\kappa\beta(J-S).$ (2)

We use gray, constant opacities with a Planck source function $S=(p/\rho)^{4}$
for the radiative transfer in these models. All simulations are run with
static mesh refinement and a maximum resolution of $256$ cells$/H_{0}$ until
at least time $t=400\Omega_{0}^{-1}$.

###  2.1 Model Parameters

The simulations adopt dimensionless units where disk scale height, orbital
frequency, and midplane gas density are set to unity
($H_{0}=\Omega_{0}=\rho_{0}=1$). A given model is then determined by the four
dimensionless parameters: $q_{t},\kappa,\beta,\epsilon/H_{0}$ which are
explained in more physical terms below. We use the subscript $0$ to refer to
characteristic dimensional values of the system which, by virtue of setting
$H_{0}=\Omega_{0}=\rho_{0}=1$, can be used to dimensionalize the dimensionless
variables presented for the remainder of the this work.

  * •

The planet core’s thermal mass $q_{t}$:

$q_{t}=\frac{M_{c}}{M_{\ast}}\left(\frac{H_{0}}{a_{0}}\right)^{-3}.$ (3)

This characterizes the core-to-star mass ratio in a way that quantifies the
strength of perturbation the planet exerts on the background protoplanetary
disk. Subthermal planets ($q_{t}\lesssim 1$) weakly perturb the disk whereas
superthemal planets ($q_{t}\gtrsim 1$) have masses and Hill radii large enough
to strongly perturb the surrounding disk.

  * •

The dimensionless opacity $\kappa$

$\kappa=\kappa_{0}\rho_{0}H_{0}.$ (4)

Because simulations are initialized as an isothermal stratified disk, the
dimensionless $\kappa$ can also be considered as the optical depth to the
midplane of the background disk: $\kappa\sim 0.8\tau_{\rm midplane}$

  * •

The dimensionless gas-radiation coupling $\beta$

$\beta=\frac{ac}{c_{0}}\left(\frac{\mu
m_{p}}{k}\right)\frac{T_{0}^{3}}{\rho_{0}}.$ (5)

This parameter encodes the thermal properties of the local disk and, along
with $\kappa$, couples the radiative transfer to the gas energy. Note that
$\beta$ can also be written as the ratio of energy fluxes characteristic to
the radiative and gas properties of the disk

$\beta=\frac{c}{c_{0}}\frac{aT_{0}^{4}}{p_{0}}\sim\frac{F_{\rm rad}}{F_{\rm
gas}},$ (6)

and can also be found in Jiang et al. (2012) as $\mathbb{CP}$ (see their Eq.
1).

  * •

The gravitational softening length $\epsilon/H_{0}$ defined by the adopted
planet potential:

$\Phi_{p}=-\frac{GM_{c}}{\sqrt{r_{0}^{2}+\epsilon^{2}}}.$ (7)

This is a parameter employed for numerical considerations. Ideally it would be
set to match the core radius of a forming planet. For planets like a proto-
Jupiter however, this scale lies orders of magnitude below the grid scale and
would be computational infeasible. Therefore we artificially set the parameter
to scale with the gravitational potential $\epsilon/H_{0}=q_{t}/10$ and study
the outer envelopes of such systems. While less physical, this choice at least
ensures that the models are fully resolved numerically.

For the models in this work we restrict the parameter space to models with
$\beta=1$ and test two different opacities with $\kappa=(1,100)$, and three
masses with $q_{t}=(0.5,1,2)$. We also isolate the role of background shear
and diagnose systematics of these models by running some simulations
designated -noshear. These models have rotational source terms omitted from
the momentum equation and a zero-velocity initial condition to simulate the
case of a core immersed in a static stratified medium. All models presented in
this work and their corresponding dimensionless parameters are summarized in
Table 1. Also provided in the table are possible choices of physical planet
parameters at $\sim 0.1$ and $\sim 5$au (assuming a solar-mass host star) that
coincide with the chosen dimensionless parameters. It is important to
recognize that the displayed physical parameters are not the only acceptable
choices and that quantities can be scaled up or down so long as they satisfy
equations 3-5 for a given set of dimensionless parameters.

Table 1: Table of Model Parameters

| Dimensionless Values | Consistent $\sim 5$au Values | Consistent $\sim 0.1$au Values
---|---|---|---
Name | $q_{t}$ | $\beta$ | $\kappa$ | $\epsilon/H_{0}$ | $\Sigma_{0}$ | $T_{0}$ | $\kappa_{0}$ | $M_{c}$ | $\Sigma_{0}$ | $T_{0}$ | $\kappa_{0}$ | $M_{c}$
|  |  |  |  | (g/cm${}^{2})$ | (K) | (cm2/g) | ($M_{\oplus}$) | (g/cm2) | (K) | (cm2/g) | ($M_{\oplus}$)
b1k1 | 0.5 | 1 | 1 | 0.05 | 250 | 65 | 0.01 | 8 | $1.5\times 10^{4}$ | 1800 | $0.017$ | 3
b1k1-noshear | - | - | - | - | - | - | - | - | - | - | - | -
b1k1-q1 | 1 | - | - | 0.1 | - | - | - | 16 | - | - | - | 6
b1k1-q2 | 2 | - | - | 0.2 | - | - | - | 32 | - | - | - | 12
b1k100 | 0.5 | - | 100 | 0.05 | - | - | 1 | 8 | - | - | 1.7 | 3
b1k100-noshear | - | - | - | - | - | - | - | - | - | - | - | -

Note. — cells marked with - are assumed to have the value of the preceding row

###  2.2 Flow Structure

To provide an example of the setup and value of these models, two snapshots
taken at the end of the fiducial simulations b1k1 and b1k100 are shown in Fig.
1.

Figure 1: Top-down snapshots of simulations b1k1 and b1k100 taken at
$t=400\Omega_{0}^{-1}$. In green is a volume rendering of passive scalars but
on a greatly exaggerated scale. In reality, the spiral streamers are orders of
magnitude lower density than the core material. The scalars are superimposed
upon a slice of the simulation midplane where temperature and velocity
directions are displayed.

The top panels display the full simulation width while the bottom panel zoom-
in to the scale of the planet’s Bondi radius $R_{b}$. Passive scalars shown in
green allow for a visualization of the recycling of circumplanetary material
along horseshoe orbits back to the protoplanetary disk. While the optical
properties of the volume rendered scalars are the same between both displayed
models, the scaling is exaggerated relative to scalar density to highlight the
recycling process. In fairness, the tenuous streams of recycled material are
substantially lower density than the saturated inner core. Also shown is the
midplane temperature structure and flow field which appears similar to the
results obtained in other works.

These simulations, even with more sophisticated thermodynamics, retain a
qualitatively similar flow field to other $\beta$ cooling models (Kurokawa &
Tanigawa, 2018) or even some isothermal ones (Ormel et al., 2015; Fung et al.,
2015; Bailey et al., 2021). In the midplane, at larger lateral distance the
flow reduces to retrograde oriented shearing streamlines of the background
disk. Along $x=0$, the models recover the typical family of horseshoe
streamlines. Closer to the planet are recycling flows sourced from horseshoes
at higher altitude. These outflows through which material is being returned to
the background disk are markedly cooler than the background equilibrium and
lie just interior to the heated spiral arms. Though this cooling is far more
obvious in the $\kappa=1$ case.

The circumplanetary flow patterns at least qualitatively also match those
found in other radiation-hydrodynamics models of atmospheric recycling
(Moldenhauer et al., 2021, 2022). Akin to other works, we fill the planet’s
Bondi sphere with passive scalars to visualize the recycling. Because our
simulations are not in thermodynamic equilibrium however, the passive scalars
are also tracking the accretion process, lessening their overall effectiveness
as a recycling diagnostic. Upon injecting the passive scalars ten orbits after
simulation start, material outside of $R_{b}/2$ in all simulations is recycled
within several orbits due to the short recycling timescales there. After
$400\Omega_{0}^{-1}$, as is shown in Figures 1 & 2, passive scalars remain
only within the innermost $R_{b}/5$.

Figure 2: Vertical slices showing the evolution of scalar concentrations for
models b1k1 and b1k100. Corresponding toroidal and poloidal (vectors)
velocities at $t=400\Omega_{0}^{-1}$ are plotted in the bottom panel.

We emphasize however that because the models are not in thermodynamic
equilibrium $R_{b}/5$ is not representative of a characteristic recycling
radius. While it may appear that passive scalars in the region between
$R_{b}/5<r<R_{b}/2$ have disappeared between $t=100\Omega_{0}^{-1}$ and
$400\Omega_{0}^{-1}$ perhaps due to recycling, there is also ongoing
accretion. Rather than being lost to recycling, said scalars have been
accreted radially into the core. Indeed, out of the scalar mass within
$R_{b}/5$ at simulation end, roughly one-fourth to one-third is additionally
accreted scalar material with the rest being unrecycled mass remaining from
the initial injection. We also find that while scalar material is
predominantly lost to the background disk via the $x$ and $y$-directions,
accreted scalar fluxes are largest in the vertical direction – reminiscent of
the pattern of meridional circulation. Between the different opacity models,
the passive scalar field appears qualitatively similar even though the
optically thicker case exhibits small scale turbulent velocities. The
optically thinner case does show stronger retention of passive scalars near
the midplane, a feature also appearing in the Moldenhauer et al. (2022)
models.

Apart from the midplane, the models exhibit flow aspects significantly
different from existing cases however – aspects which we believe ultimately
impact the implications for planet formation. Viewing the vertical structure
of the velocity field in Fig. 2, we see similar toroidal rotation patterns
independent of the opacity. The circumplanetary region, sourced from high
altitude prograde oriented streamlines, retains a prograde rotating atmosphere
that is highly sub-Keplerian. The poloidal velocities however transition from
strong polar inflow and are dramatically damped within the inner $R_{b}/2$.
This appears reminiscent of the work of Kurokawa & Tanigawa (2018) where
envelope models that include cooling, establish a strong entropy gradient in
the atmosphere. This entropy gradient results in an additional buoyant force
upon a fluid element impinging upon the inner envelope and creates a barrier
to recycling. This in in contrast to both adiabatic and isothermal models
which do not have cooling and thus do not set up strong entropy gradients. In
these models, recycling appears much more capable in accessing the inner
envelope which would act to interfere with the quasi-static contraction
process. Though we use a more realistic cooling by doing the radiative
transfer, the buoyancy barrier phenomenon put forth by Kurokawa & Tanigawa
(2018) appears consistent with the models presented here.

##  3 Mass Exchange

Recycling flows tend to exchange high entropy disk material with the cooled
low entropy gas within the planet atmosphere. This effectively heats the outer
regions of the planetary envelope to an equlibrium where mass can no longer be
accreted. Within this context, it is natural to examine the properties of mass
accretion in our models and study locations where accretion has been inhibited
to determine the influence recycling flows have on envelope growth.

In Figure 3 we plot the mass interior to several radii for the fiducial b1k1
and b1k100 models.

Figure 3: Evolution of the envelope mass as a function of time for models b1k1
and b1k100. Three choices of envelope radius ($r=R_{b}$, $R_{b}/2$, $R_{b}/5$)
are shown and compared against the mass in an equivalent -noshear model.

At early times, on the dynamical timescale, the mass increase occurs as the
envelope reaches an approximate hydrostatic equilibrium in response to the
newly introduced planet potential. On longer timescales, radiative cooling in
the interior drives the additional accretion of mass. Notably, we see that the
envelope masses measured are still increasing even at the end of the fiducial
runtime $t=400\Omega_{0}^{-1}$. This contrasts the simulations of Moldenhauer
et al. (2022) where thermodynamic equilibrium and steady state are achieved
after $\sim 300\Omega_{0}^{-1}$. We also emphasize that the relative mass
change in logspace at the largest scale (blue curves) is small but the
relative mass change at the small scale (green curves) is large, suggesting
that significant mass accretion is occurring in the interior envelope.

While Fig. 3 demonstrates that envelope masses are increasing for the duration
of the simulations, it is not sufficient to determine the extent of the
envelopes. After all, the mass increases for all curves could be a reflection
of only the innermost regions accreting. To address this, we present Fig. 4 in
which integrated net mass fluxes at a given distance from the center of the
simulation domain (i.e. planet center) are mapped out as a function of time.

Figure 4: Top: integrated mass fluxes through a surface located at $r$ for
models b1k1 and b1k100. Middle: $\dot{M}$ of a shell of single grid cell
thickness located at $r$. Bottom: The corresponding $\dot{M}_{\rm
shell}/M_{\rm shell}$ at $r$. The geometry in all cases is a box of half-width
$r$ to match the cartesian simulation domain.

Moving outwards from the planet center we construct concentric boxes aligned
on cell interfaces and calculate the net mass fluxes through each box. These
integrated fluxes are plotted in the upper panels of Fig. 4 with the ordinate
giving the half-width of the box (in units of $H_{0}$) used to perform the
surface integral. Locations marked blue possess a net inward flux of mass at
corresponding time while those marked red show an outward flux. And by the
divergence theorem, these fluxes are equal albeit opposite to the net
accretion rate $\dot{M}$ taking place within the constructed box.

We find that the fluxes generally end up pointing inward everywhere, a finding
reflected by the positive derivatives of the mass curves in Fig. 3. The full
flux map of Fig. 4 however, also contains gradients in the $r/H_{0}$
direction, suggesting variations in the mass accretion rate for different
distances from the planet core. For the $\kappa=1$ model, the fluxes appear
nearly constant at large distances from the planet, suggesting no net mass
accretion in those regions. For more interior regions where gradients in the
flux appear steeper there is greater associated mass accretion. The
$\kappa=100$ model appears to show more extended gradients in flux, but as we
show in further analysis, these are more an artifact of the b1k100 fluxes
existing in the linear part of the log-linear color normalization than a true
difference between the models. Similarly the near zero flux at $r\sim 0.01$
appears to be an artifact from our post-processing of the fluxes that arises
when close to the grid scale ($\sim 1/256$).

Because density changes are related to fluxes through the continuity equation
$\dot{\rho}=-\nabla\cdot F$, the local accretion rates can be better presented
and quantified if instead of considering the fluxes, flux gradients are
considered. This is reflected in the middle panels of Fig. 4 where plotted are
the differences in fluxes between the adjacent cell interfaces plotted in the
upper panels. By the continuity equation, these plotted flux differences are
equivalent to the mass accretion rates $\dot{M}_{\rm shell}$ in concentric
shells of single cell width thickness. To make the measurement independent of
resolution, we also supply $\dot{M}_{\rm shell}/M_{\rm shell}$ in the lower
panels of Fig. 4. Physically this would represent some inverse mass accretion
timescale for a given radial location. The middle and lower panels of Fig. 4
show that by the end of simulation runtime, there is an inner core where
$\dot{M}$ is solidly red and still increasing. Exterior, there is a recycled
region consistent with zero accretion. This outer recycled region is not
quiescent but shows temporal variation between local net mass accretion and
loss. The variation appears to be due to unsteady fluid motion in the envelope
also reflected in the turbulent looking outflow streams in Fig. 1.

The main advantage of Fig. 4 is that the transition region from accreting to
non-accreting can be quantified. To define this transition radius for a given
simulation snapshot we take the minimum radius at which a shell switches to
mass loss, effectively taking a contour on Fig. 4 for $\dot{M}_{\rm shell}=0$.
With this definition, we find that by the end of simulation the transition
radius levels out to a value $r_{t}\sim H_{0}/5\sim 0.4R_{b}$ for both
simulations depicted in Fig. 4. Repeating this procedure for our other models,
we find that the transition radius scales linearly with $q_{t}$ and therefore
linearly with the planet mass. These transition radii are plotted accordingly
in Fig. 5.

Figure 5: Box half-widths where $\dot{M}_{\rm shell}$ vanishes (blue) or at
least becomes less than $10^{-6}$ (orange). Values for $r/H_{0}$ come from
models b1k1, b1k1-q1, b1k1-q2 averaged over the last $100\Omega_{0}^{-1}$ with
errorbars giving the standard deviation over the final $100\Omega_{0}^{-1}$
and therefore some estimate of the time variability. Black dashed lines are
drawn as references for linear scaling $r/H_{0}\propto q_{t}$

To give some account of the variability, we average the value of transition
radius over the final $100\Omega_{0}^{-1}$ and present the standard deviation
as error bars. In addition to points for $\dot{M}_{\rm shell}=0$, points for
the radius where $\dot{M}_{\rm shell}=10^{-6}$ are also plotted in Fig. 5.
While less physically relevant, the $\dot{M}_{\rm shell}=10^{-6}$ boundary is
where the color-scale in Fig. 4 switches from log to linear and therefore can
be readily discerned by eye in Fig. 4. The $\dot{M}_{\rm shell}=10^{-6}$
transition location shows a similarly approximate linear scaling with the
thermal mass demonstrating the insensitivity of the results to a precise
definition of this transition radius.

##  4 Energetics

In addition to transferring mass, atmospheric recycling can transfer energy
between the forming envelope and background disk. Because accretion is
governed by the cooling timescale for the bulk of protoplanet lifetimes, one
can also hope to study the envelope accretion process by looking at the
envelope energetics. In this section we study the energy budget of our
simulations by considering three energy sources/sinks and testing systematics
against our -noshear models. The energy sources considered come from a rewrite
of the energy conservation equation (2). To do so, one can split the potential
into planetary and stellar components $\Phi=\Phi_{p}+\Phi_{s}$ where
$\Phi_{p}=q_{t}/\sqrt{r^{2}+\epsilon^{2}}$ and $\Phi_{s}=(3x^{2}+z^{2})/2$.
The energy variable can then be redefined to include the gravitational
potential energy of the planet $E\rightarrow p/(\gamma-1)+\rho
v^{2}/2+\rho\Phi_{p}$. With this definition, the energy appears more in
keeping with that of standard 1D models (but with also the kinetic energy
included). Using the fact that $\partial\Phi/\partial t=0$, the full energy
equation (2) can be rewritten as

$\dot{E}=\dot{E}_{\rm adv}+\dot{E}_{\rm rad}+\dot{E}_{\rm work}$ (8)

with the three component terms described as follows.

$\dot{E}_{\rm adv}=-\nabla\cdot(Ev+pv)$ (9)

is an advective energy term that tends to be positive and act as an energy
source, heating the envelope and slowing accretion. Being purely the
divergence of an energy flux, the advective energy is simply the portion of
the total energy that obeys energy conservation simply by considering the
continuity of energy moving into or out of a given volume. This advective
transport can be due to both accretion with accreted material adding its own
energy to the envelope or recycling removing low entropy material and
replacing it with high entropy disk gas. In the -noshear simulations where
recycling is minimal, this term ought to be dominated by the accretion
component, giving also some estimate of the relative importance of accretion
vs. recycling in the fiducial runs.

$\dot{E}_{\rm rad}=\rho\kappa\beta(J-S)$ (10)

is an energy sink term due to radiative cooling given by our solution of the
time-independent transfer equation.

$\dot{E}_{\rm work}=-\rho v\cdot\nabla\Phi_{s}$ (11)

is an energy term acquired by doing work against the vertical stellar and the
background shear components of the potential. Because these components of the
potential scale as $z^{2}$ and $x^{2}$ respectively, this term tends to be
most relevant at large $\sim H_{0}$ scales.

In the previous section, and in keeping with the finite volume methods of
simulation, we studied integrated mass measurements in concentric volumes
centered on the planet. Similarly here we consider the energy budget of these
planet envelopes by integrating each of the three energy terms $\dot{E}_{\rm
adv}$, $\dot{E}_{\rm rad}$, $\dot{E}_{\rm work}$ over equivalent control
volumes. In Figures 6 & 7 we plot these energy terms after integrating over
control volumes of half-width $R_{b}/5$, $R_{b}/2$ and $R_{b}$. Figure 6 shows
the energy budgets for the b1k1 and b1k1-noshear models while Figure 7
contains the results for the optically thicker k100 equivalents. As a check,
we compare the sum of the components $\dot{E}_{\rm adv}+\dot{E}_{\rm
rad}+\dot{E}_{\rm work}$ calculated by post-processing the primitive variables
against the true $\dot{E}$ measured during the simulation. In all cases, the
$\dot{E}$ measured from each of these calculations agree sufficiently well.

Figure 6: Energy sources and sinks integrated over cubes of half-width
$R_{b}/5$ (left), $R_{b}/2$ (center), and $R_{b}$ (right). Dotted lines are
negative in sign. Top panels are taken from the fiducial k1b1 model and bottom
panels from the equivalent -noshear model.  Figure 7: Same as Figure 6 but for
the b1k100 models.

From the rightmost panels of Figures 6 & 7 we see that on a scale comparable
to the canonical planet envelope $R_{b}$, the energy budget is predominantly
determined from a balance of heating from the advective term (orange curves)
and radiative cooling (green curves). The difference between the two is the
subsequent net envelope cooling (blue curves) which regulates accretion.
Judging by the similarity between $\dot{E}$ curves irrespective of the choice
of $r=R_{b}/5$, $R_{b}/2$, $R_{b}$ is being dominated by cooling occurring at
small $r$. This is consistent with the idea that recycling is acting to bring
the outer envelope to steady state while the inner envelope continues to
radiative cool and accrete. Furthermore, as one goes to smaller radius, the
difference between advection and cooling is larger. This suggests a point deep
enough in the envelope at which radiative cooling fully dominate over
advection and a Kelvin-Helmholtz like contraction can be recovered.

To be fair, the leftmost panels showing the energy budget in the very interior
envelope do present a non-zero advective term. Though the advective flux term
associated with recycling tends to be smaller than the net $\dot{E}$ by 1-2
orders of magnitude by the end of the simulation, there remains the question
whether this would persist over longer timescales. However, it appears that
this $\dot{E}_{\rm adv}$ on the interior scales cannot be even be attributed
to recycling if we consider the -noshear runs. After all, the advective energy
term does not solely include energy advected due to recycling but also energy
advected due to accretion. In -noshear runs where we turn off the shearing
terms and have a roughly spherical accretion of material, the role of
recycling is expected to be limited and the $\dot{E}_{\rm adv}$ ought to be
predominantly due to the accretion. If the inner envelope energy budgets are
compared to the lower panels showing the equivalent -noshear runs, we see the
magnitude of $\dot{E}_{\rm adv}$ in the fiducial runs is comparable to or even
lower than the $\dot{E}_{\rm adv}$. This would suggest that the $\dot{E}_{\rm
adv}$ in the fiducial runs is not actually heating due to recycling at all but
rather heating due to the accretion of high energy gas from larger radius.
This leads us to conclude that at least at $R_{b}/5$, the energetics show no
sign of being effected by recycling dynamics. This would be consistent with
the results of Section 3, that suggest a recycling radius at $\sim 0.4R_{b}$.

An additional interesting feature of this energy budget analysis is that on
the larger scales, the magnitude of the radiative cooling term (i.e. the
luminosity), grows by one to two orders of magnitude compared to the actual
energy liberated by Kelvin-Helmholtz contraction. In a standard quasi-static
picture one would expect these two to be equal – any observed luminosity would
be a direct result of the contraction. In these models, the recycling
occurring in the outer envelope is continually tapping into the fresh hot
reservoir of disk gas, converting that energy to radiation, and adding to the
luminosity from KH contraction. The result is a significantly higher
luminosity than otherwise expected from a forming protoplanet. This could
potentially make embedded protoplanets more observable than previously
thought, although separating the signal of KH contraction from that of the
recycling would be another confounding factor in understanding the true
accretion of the planet. Stated another way: were an embedded protoplanet
observed, the true KH contraction could be significantly slower than inferred
by the simple conversion of observed luminosity to liberated gravitational
binding energy.

##  5 discussion

###  5.1 Connection with 1D Models

Inherent to the 3D simulation of planet envelopes is the problem that these
models remain prohibitively expensive to run for the end-to-end formation of
the planet. Nevertheless, we can develop 1D models with parameterized 3D
effects and simulate for the entire planet formation timescale. In the
process, it is also possible to mitigate some of the shortcomings of the 3D
models such as the assumption of a softened gravitational potential and the
lack of self-gravity.

The results of the preceding sections indicate that some interior portion of
the envelope models continues to cool and contract, suggesting that a 1D
quasi-static treatment could be appropriate to model the long-term evolution.
In Figures 8 & 9 we present snapshots of the spherically averaged mass,
temperature, entropy, and luminosity profiles derived from our 3D models. In
the quasi-static problem, the envelope structure at a given time is entirely
determined by the corresponding luminosity profile. This appears to also be
the case in our 3D simulations and can be seen in Figures 8 & 9 where also
plotted are the envelope profiles retrieved by applying the hydrostatic
equations to the measured luminosity profiles. The dashed hydrostatic profiles
are produced using only the spherically averaged 1D luminosity profile and 1D
structure equations and match the averaged 3D results very well. This provides
some support to the idea that a quasi-static 1D framework can be applied to
adequately model the evolution. Furthermore, if one attempts to reproduce our
averaged model profiles with the standard 1D methods (Pollack et al., 1996;
Piso & Youdin, 2014), the fit (not pictured) is reasonable (further suggesting
that these models are participating in a sort of 1D Kelvin-Helmholtz
contraction). That being said, 1D quasi-static models can be improved with
some minor effort, informed by the 3D envelope energetics, to better match our
3D results. In the following we describe the process and rationale for
including these multi-dimensional hydrodynamical effects in 1D models before
applying them to the long-term growth of a core under Jupiter-like conditions
as well as conditions at 0.1 AU characteristic of close-in super-Earth/mini-
Neptunes.

Figure 8: Spherically averaged mass, temperature, entropy and luminosity
profiles for the b1k1 models (solid). Also shown (dashed) are the mass,
temperature, and entropy profiles that would be obtained from applying an
assumption of hydrostatic equilibrium to the measured luminosity profiles.
Figure 9: Same as Figure 8 but for the optically thicker b1k100 models.

For a first modification to 1D evolutionary models we consider the results of
our mass and energy budgeting. It seems that material exterior to $\gtrsim
0.4R_{b}$ is recycled efficiently, producing a net accretion rate consistent
with zero in that region. On the other hand, material interior to this radius
continues to cool and accrete. This is reflected in the bottom panels of Fig.
2 as the location interior to which the poloidal velocity is damped to zero,
in Fig. 4 as the location interior to which the envelope continues to accrete,
and in the leftmost panels of Figs. 6 & 7 where the radiative cooling
dominates the envelope energetics. This prompts us to effectively place the
canonical outer boundary of the cooling and accreting 1D envelope at
$0.4R_{b}$ (with the scaling here being informed by the results of Fig. 5).

Moving the outer boundary of the canonical envelope from its standard location
of $R_{\rm out}=\min\left(R_{b},R_{h}\right)$, to $0.4R_{b}$ prompts an
additional modification, as one must now set a value for the pressure and
temperature at $0.4R_{b}$ to fully determine the boundary condition on the 1D
model. One option is to say the envelope exterior to $0.4R_{b}$ is isentropic,
reflecting the idea that the recycled outer envelope is continually being
replenished with high-entropy disk material. Pressure and temperature can then
be integrated along this isentrope from the background disk to derive the
required values at $0.4R_{b}$ e.g. Ali-Dib et al. (2020). This results in what
we will refer to as a two-layer recycling envelope model – an outer isentropic
recycling region and an inner Kelvin-Helmholtz contracting core. While this is
not a terrible approximation, it can be seen in both Figures 8 & 9 that the
entropy in the 3D models is not constant exterior to $0.4R_{b}$ (i.e.
$0.2H_{0}$). Rather, a recycled fluid element does undergo some amount of
cooling in its journey from disk to planet and back. While mass is conserved
under advection, energy is certainly not as it can be lost through radiative
processes. This is reminiscent of the results of Zhu et al. (2021) where the
thermal timescale appears to be shorter than the recycling timescale resulting
in an envelope where mass is recycled exterior to $0.1R_{h}$ yet the thermal
structure remains unaltered. Here, the models similarly cool exterior to the
fiducial $0.4R_{b}$ core even though mass can be rapidly replenished there.
The end result is outer steady-state envelope region that is constantly cooled
and recycled in equilibrium, with the innermost regions being the most cooled.
This latter point is perhaps somewhat counter-intuitive as the envelope
interior is higher density and therefore has higher optical depth and slower
cooling. However, the recycling in the inner envelope also occurs more slowly,
with a given fluid element potentially making many orbits around the planet
and consequently extending the time available to cool.

While one could hope to describe this cooled but recycled region by equating
the recycling and cooling times for given fluid elements, the recycling flow
pattern is very much 3-dimensional and also dependent on the opacity itself
making a simple analytic description for the energetics of this equilibrium
challenging. In Zhu et al. (2021) for example, three estimates for the radial
thermal timescale are provided, all of which would give a different 1D
parameterized boundary if equated to the recycling timescale here. One
effective but simple parameterization we have found to model this recycled but
cooled region in 1D is to subdivide it into two regions – an outer part where
the recycling time is shorter than the cooling time, is made to be isentropic
and connected to the background nebula as in the two-layer model. The inner
part, where recycling is slower than cooling (but still steady-state), is made
to evolve for a few $10^{2}-10^{3}\Omega_{0}^{-1}$ and taken as a fixed
recycling-cooling equilibrium profile. The resultant equilibrium profile
appears fairly insensitive to the choice of the evolution time but from a
physical point of view should be on the order of a characteristic recycling
time. The only free parameter in this method is the radius at which the
recycled region switches from isentropic to the cooled-recycled equilibrium.
This parameter can be estimated from the 3D simulations and for $\kappa=1$
models is taken as $H_{0}$ and for $\kappa=100$ models is taken as $0.6R_{b}$.
Everything taken together results in a 1D three-layer recycling model as
depicted in 10 which matches the 3D spherically averaged profiles better than
classical quasi-static calculations and two-layer recycling models.

Figure 10: Diagram representing the structure of our three-layer recycling 1D
models.

We recap and summarize the steps needed to reproduce our three-layer recycling
models as follows.

  1. 1.

Integrate the nebular pressure and temperature, assuming isentropy, from the
background disk111We adopt $H_{0}$ as nebular boundary condition radius. For
$\kappa=1$ models, this means there is formally no isentropic outer layer.
However, the results are insensitive to this choice and $1.5H_{0}$ could have
been chosen for equivalent results with a truly three-layer model. to the
outer radius of the recycling-cooling equilibrium layer (estimated from 3D
simulations as $H_{0}$ for $\kappa=1$, and $0.6R_{b}$ for $\kappa=100$).

  2. 2.

Evolve the quasi-static model for $\sim 10^{2}-10^{3}\Omega_{0}^{-1}$ to get
the steady-state profile for the recycled-cooled equilibrium region.

  3. 3.

Take the pressure and temperature at $R_{\rm out}=0.4R_{b}$ as an outer
boundary condition to evolve the Kelvin-Helmholtz contraction of the inner
envelope under the standard quasi-static framework.

With this 3D motivated methodology, we study the long-term evolution of
planetary cores at Jupiter and close-in super-Earth distances and compare to
existing calculations. To integrate the quasi-static equations for our cores
we use the method of Piso & Youdin (2014) which integrates the usual structure
equations:

$\frac{dm}{dr}=4\pi r^{2}\rho$ (12) $\frac{dp}{dr}=-\rho\nabla\Phi_{p}$ (13)
$\frac{dT}{dr}=\nabla\frac{T}{p}\frac{dp}{dr}$ (14)
$\nabla=\min\left(\nabla_{\rm ad},\frac{3\kappa pL}{64\pi Gm\sigma
T^{4}}\right)$ (15)

but makes the assumption of constant luminosity thus becoming an eigenproblem
for the eigenvalue $L$ given an envelope mass. The sequence of time-
independent structure models are then stitched together using the global
energy equation

$L=\frac{dE}{dt}$ (16)

to determine the temporal evolution.

With this methodology we explore the long-term growth of cores subject to
different initial conditions and different recycling parameterizations. For
each core we run a standard 1D model with no recycling parameterization
($R_{\rm out}=\min(R_{b},R_{h})$), a two-layer recycling parameterization
($R_{\rm out}=0.4R_{b}$ with the boundary condition determined by an
isentropic layer extending to $H_{0}$), and our three-layer recycling
parameterization. All 1D models remove the gravitational softening treatment
and instead include self-gravity and are integrated from an inner boundary
placed at the core radius assuming core density $\rho_{c}=3$ g/cm3. All models
also adopt $\gamma=1.4$ and $\mu=2.3$, similar to Ali-Dib et al. (2020).

Under Jupiter-like conditions ($M_{\ast}=M_{\odot}$, $T_{0}=65$ K,
$\Sigma_{0}=250$ g/cm2, $a=5$ AU), we display the accretion history of an
$8M_{\oplus}$ core at $10^{-2}$ and $1$ cm2/g (corresponding to dimensionless
$\kappa=1$, $\kappa=100$) in Fig. 11.

Figure 11: Evolution of an $8M_{\oplus}$ core at 5 AU. The $\kappa_{0}=1$
cm2/g opacity is characteristic of Jupiter-like conditions while
$\kappa_{0}=10^{-2}$ cm2/g is only included in analogy to the 3D simulations.
Solid lines are standard quasi-static models where $R_{\rm
out}=\min(R_{b},R_{h})$, dashed lines are our three-layer recycling models,
and dotted lines are an equivalent two-layer recycling model.

The opacity at $1$ cm2/g is characteristic of Jupiter-like conditions and
reaches a phase of runaway accretion after $1$ Myr consistent with existing
models. The low $10^{-2}$ cm2/g opacity is more unphysical and shown more for
the sake of completeness to apply our b1k1 simulations. Under these
conditions, the evolution is relatively unchanged by the inclusion of a two or
three-layer recycling parameterization, changing growth times by less than a
factor of two. This is due to the fact that the envelope energetics are
dominated by a high mass convective core interior to $0.4R_{b}$ consistent
with other works (Ali-Dib et al., 2020) that employ a two-layer recycling
framework and find no significant effects at $5$ AU.

In the inner disk however, recycling can act much more strongly to limit
envelope growth. In Fig. 12, we provide the accretion history of a set of
models subject to hypothetical conditions ($M_{\ast}=M_{\odot}$, $T_{0}=2150$
K, $\Sigma_{0}=2.5\times 10^{4}$ g/cm2) at $0.1$ AU.

Figure 12: Evolution of an $4.3M_{\oplus}$ core (left) and $6M_{\oplus}$ core
(right) at 0.1 AU. Presentation is similar to that of Fig. 11. Disk conditions
at $0.1$ AU for these models are chosen as $T_{0}=2150$ K and
$\Sigma_{0}=2.5\times 10^{4}$ g/cm2 so that $\beta=1$. Several of the low-mass
curves fail to runaway which is a consequence of the more approximate method
and the imposed constant luminosity but are included for the sake of
completion.

Again we adopt constant opacities corresponding to dimensionless opacities
$\kappa=1,100$ which under the adopted disk parameters is $10^{-4}$ and
$10^{-2}$ cm2/g. While the $10^{-4}$ cm2/g opacity might appear low relative
to the assumed background temperature and density, the dust properties of
these protoplanet envelopes are poorly understood and indeed models of grain
settling in these environments (Ormel, 2014; Mordasini, 2014) predict
significantly reduced opacities.

When modeling a $4.3M_{\oplus}$ core under these conditions and no recycling
(left panel of Fig. 12), the low opacity envelope easily runs away to gas
giant within a disk lifetime. With a three-layer recycling model this
timescale is lengthened by a factor $\sim 2$. The other recycled models for a
$4.3M_{\oplus}$ core reach a breaking point before the runaway phase sets in.
This appears to be a consequence of the more approximate method and imposed
constant luminosity constraint that could be addressed in the future with a
coupled time-dependent treatment e.g. Bodenheimer & Pollack (1986); Pollack et
al. (1996); Ikoma et al. (2000). In the meantime, we run instead a higher mass
$6M_{\oplus}$ core which can robustly sustain a hydrostatic constant
luminosity atmosphere (shown in the right panel of Fig. 12). The low-opacity
calculation similarly runs away within a disk lifetime when no recycling is
included. A two-layer recycling treatment however lengthens the growth time by
nearly two orders of magnitude. At the same time, a three-layer recycling
model only lengthens it by one order of magnitude. The appears to be a fairly
generic hierarchy in that in all cases a two-layer recycling will act to
lengthen the growth timescale for core and a three-layer recycling will only
act to push it back to a more intermediate value. In this recycling framework,
given a fixed recycling radius, a no-recycling model gives a lower bound on
the KH contraction timescale while a two-layer model produces an upper bound.

In most cases examined here, there is no significant difference between the
two and three-layer treatment, but for lower opacities in the inner disk the
distinction could be important like in the $10^{-4}$ cm2/g, $6M_{\oplus}$
case. While these results are more pessimistic about the ubiquity of recycling
halting the growth of all super-Earths than say Moldenhauer et al. (2022),
recycling only the outer $0.6R_{b}$ can still slow the accretion process by a
significant amount depending on the formation conditions. In the future, a
full population synthesis would do well to determine the entire scope of this
process as well as the distinction between a two/three-layer treatment. Until
that time, here there is at least some proof-of-concept.

###  5.2 Comparison with Other Works

Several existing works (Moldenhauer et al., 2021, 2022) have investigated the
radiation-hydrodynamics of protoplanetary envelopes in a manner similar to
this work but obtained results different from those here. After several
hundreds of dynamical times, the models presented in those works appear to
reach a steady state and stop accreting. This would lead to the conclusion
that super-Earth envelopes similarly be entirely recycled thus stalled before
hitting the stage of runaway accretion. This would naturally explain the
envelope mass fractions at the level of a few percent as a similar mass is
present in the model envelopes once a steady state is reached. Our models do
not appear to reach a steady state but rather the inner $\sim 0.4R_{b}$ show
quite good agreement with a 1D quasi-static Kelvin-Helmholtz contraction. We
have not been able to isolate the cause of this difference between our results
and previous work but here examine some differences and comment on them.

A major difference between the two works stems from the parameter space of the
models studied. In our dimensionless units, the Moldenhauer et al. (2021)
models are run with $q_{t}=0.3$, $\kappa=20$, and most notably
$\beta=4.5\times 10^{-3}$ primarily due to a cooler protoplanetary disk model.
This is several orders of magnitude below $\beta=1$ used here, and has a far
slower cooling rate. In light of the results of Kurokawa & Tanigawa (2018), it
is entirely possible that $\beta=1$ conditions result in a quasi-static
interior with $\beta\ll 1$ being fully recycled. Said work showed that
adiabatic thermodynamics with $\beta=0$ (note this is our definition of
$\beta$ and differs from the referenced work) are uninhibited by the entropy
gradients and associated buoyancy barrier that arises in models with non-zero
$\beta$. While we have tested some models with an orders of magnitude smaller
$\beta$ (Bailey et al., 2023), they appear essentially indistinguishable from
an equivalent adiabatic simulation. In such a case it is difficult to
determine if the cooling is truly resolved or if the radiative source term is
being dominated by numerical errors. Hence we do not include those models here
and restrict this paper only to models that are confidently converged where
the energy budget analysis of Sec. 4 can be applied reliably.

Another difference between these models and others comes from the solution of
the radiative transfer. While other works adopt the popular method of flux-
limited diffusion, we calculate the formal solution to the time-independent
transfer equation by discretizing along short characteristics (Davis et al.,
2012). Unlike M1 or diffusion approximations, by solving the transfer equation
directly, the short characteristics method properly treats any anisotropies in
the radiation field. However, in light of the results of Bailey et al. (2023)
which suggest only minor differences between FLD and these methods, we do not
believe this accounts for the discrepancy.

Another potential cause lies in the choice of coordinate system and associated
boundary conditions. While Moldenhauer et al. (2021, 2022) adopt a spherical
polar system centered on the planet, our models utilize a Cartesian one. While
we don’t believe the coordinate system to have such a substantial effect on
the problem, we imagine the associated boundary conditions could have some
effect. In spherical polar one may adopt a reflecting boundary at the
planetary surface, a good proxy for the central rocky core of the planet. In
Cartesian coordinates, the models are precluded from having a boundary at the
planet’s location and must adopt some gravitational softening length instead.
In having no physical boundary, there is an additional reservoir of material
interior to the physical core radius is able to cool and accrete – potentially
lengthening the time to recycle. This would not entirely halt recycling
however as material can still flow interior to the softening length. Given
sufficiently long time it is possible that non-boundary models could recover
the full steady state seen in spherical polar. If this were the case, it would
actually provide a good mechanism for bifurcating the planet formation process
into forming super-Earths at small semi-major axis and gas giants at large
semi-major axis. This is because for a given core mass, the relative size
difference between core radius and envelope radius grows as one goes further
out in the disk. In this way, a core with no boundary and its reservoir of
interior mass would be akin to a core moved further out in the disk. If the
non-boundary model then showed slower recycling it would be akin to cores
further out in the disk also exhibiting slower recycling. At Jupiter
distances, perhaps this recycling would become so slow that runaway accretion
would be recovered with the bifurcation of planet outcomes being related
mainly to semi-major axis of formation. In the future some effort could be
made to compare a diffusion method between coordinate systems or attempt to
model a free reflecting sphere in Cartesian to test this.

A final hypothesis for the discrepancy in model outcomes comes from resolution
effects. As seen in Bailey et al. (2023), resolution can have a significant
impact on the inferred accretion rates in these types of models with under-
resolution resulting in a false steady state. In this work, at the cost of
long simulation runtime and physically small softening length, we have focused
on obtaining converged results, running at a high enough spatial resolution to
obtained converged accretion rates informed by convergence testing in our
previous work. We emphasize that higher-order methods are essential for this
kind of modeling as piecewise linear methods at the resolutions employed here
would be insufficient to resolve the long term cooling we have measured.

###  5.3 Summary

We have presented 3D models of envelope growth evaluating the role of mass and
energy transport between envelope and protoplanetary disk for several choices
of planet mass and opacity. For the first time we have conducted an energy
budget analysis to determine how robustly energy from the background disk is
transported into the interior envelope through recycling. By looking at the
mass and energy budgets of the 3D models, we have determined that while there
is obvious recycling forcing the envelope to steady state out to $\sim
0.4R_{b}$, this effect does not persist down to the inner envelope where a
Kelvin-Helmholtz contraction appears to be recovered. While the 3D models make
a number of simplifying assumptions including a large softening length,
constant opacity, etc., this allows us to ensure the measured accretion rates
are converged. With lessons learned from these simplified 3D models, we have
extended the analysis back to more realistic cases in 1D with the development
of a three-layer recycling model. This three-layer model accounts for
differences in the way mass and energy are influenced by the recycling process
rather than assuming isentropy. With these 1D models we conclude that the
hydrodynamic effects would have had relatively little impact in determining
the formation process of planets like Jupiter. For close-in super-Earths/mini-
Neptunes however, recycling could have a more significant effect. While our
models suggest that only the outer $0.4R_{b}$ is recycled, growth times can
still be delayed by one to two orders of magnitude depending on the core
conditions and the recycling parameterization. In general, our preferred
three-layer recycling models tend to delay the growth rates of protoplanet
envelopes though not as much as an equivalent two-layer model. In the future,
we hope to push these 3D models further both in time and in terms of realism
to continue to test whether and under what conditions these models reach a
true steady state.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ali-Dib et al. (2020) Ali-Dib, M., Cumming, A., & Lin, D. N. C. 2020, MNRAS, 494, 2440, doi: 10.1093/mnras/staa914
  * Bailey et al. (2021) Bailey, A., Stone, J. M., & Fung, J. 2021, ApJ, 915, 113, doi: 10.3847/1538-4357/ac033b
  * Bailey et al. (2023) Bailey, A. P., Stone, J. M., & Fung, J. 2023
  * Bodenheimer & Pollack (1986) Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391, doi: 10.1016/0019-1035(86)90122-3
  * Cimerman et al. (2017) Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662, doi: 10.1093/mnras/stx1924
  * Davis et al. (2012) Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9, doi: 10.1088/0067-0049/199/1/9
  * Dong & Zhu (2013) Dong, S., & Zhu, Z. 2013, ApJ, 778, 53, doi: 10.1088/0004-637X/778/1/53
  * Dullemond & Turolla (2000) Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187, doi: 10.48550/arXiv.astro-ph/0003456
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102, doi: 10.1086/428383
  * Fung et al. (2015) Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101, doi: 10.1088/0004-637X/811/2/101
  * Fung & Lee (2018) Fung, J., & Lee, E. J. 2018, ApJ, 859, 126, doi: 10.3847/1538-4357/aabaf7
  * Ginzburg & Chiang (2019) Ginzburg, S., & Chiang, E. 2019, MNRAS, 487, 681, doi: 10.1093/mnras/stz1322
  * Ginzburg et al. (2016) Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29, doi: 10.3847/0004-637X/825/1/29
  * Guillot et al. (1997) Guillot, T., Gautier, D., & Hubbard, W. B. 1997, Icarus, 130, 534, doi: 10.1006/icar.1997.5812
  * Ida & Lin (2004) Ida, S., & Lin, D. N. C. 2004, ApJ, 616, 567, doi: 10.1086/424830
  * Ikoma & Hori (2012) Ikoma, M., & Hori, Y. 2012, ApJ, 753, 66, doi: 10.1088/0004-637X/753/1/66
  * Ikoma et al. (2000) Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013, doi: 10.1086/309050
  * Jiang et al. (2012) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14, doi: 10.1088/0067-0049/199/1/14
  * Kurokawa & Tanigawa (2018) Kurokawa, H., & Tanigawa, T. 2018, MNRAS, 479, 635, doi: 10.1093/mnras/sty1498
  * Lee (2019) Lee, E. J. 2019, ApJ, 878, 36, doi: 10.3847/1538-4357/ab1b40
  * Lee & Chiang (2015) Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41, doi: 10.1088/0004-637X/811/1/41
  * Lee et al. (2014) Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95, doi: 10.1088/0004-637X/797/2/95
  * Militzer et al. (2022) Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, , 3, 185, doi: 10.3847/PSJ/ac7ec8
  * Mizuno (1980) Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544, doi: 10.1143/PTP.64.544
  * Moldenhauer et al. (2021) Moldenhauer, T. W., Kuiper, R., Kley, W., & Ormel, C. W. 2021, A&A, 646, L11, doi: 10.1051/0004-6361/202040220
  * Moldenhauer et al. (2022) —. 2022, A&A, 661, A142, doi: 10.1051/0004-6361/202141955
  * Mordasini (2014) Mordasini, C. 2014, A&A, 572, A118, doi: 10.1051/0004-6361/201423702
  * Mordasini (2018) —. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 143, doi: 10.1007/978-3-319-55333-7_143
  * Ormel (2014) Ormel, C. W. 2014, ApJ, 789, L18, doi: 10.1088/2041-8205/789/1/L18
  * Ormel et al. (2015) Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512, doi: 10.1093/mnras/stu2704
  * Petigura et al. (2018) Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, AJ, 155, 89, doi: 10.3847/1538-3881/aaa54c
  * Piso & Youdin (2014) Piso, A.-M. A., & Youdin, A. N. 2014, ApJ, 786, 21, doi: 10.1088/0004-637X/786/1/21
  * Piso et al. (2015) Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82, doi: 10.1088/0004-637X/800/2/82
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Zhu et al. (2021) Zhu, Z., Jiang, Y.-F., Baehr, H., et al. 2021, MNRAS, 508, 453, doi: 10.1093/mnras/stab2517