# Comparison of model selection techniques for seafloor scattering statistics

###### Abstract

In quantitative analysis of seafloor imagery, it is common to model the
collection of individual pixel intensities scattered by the seafloor as a
random variable with a given statistical distribution. There is a considerable
literature on statistical models for seafloor scattering, mostly focused on
areas with statistically homogeneous properties (i.e. exhibiting spatial
stationarity). For more complex seafloors, the pixel intensity distribution is
more appropriately modeled using a mixture of simple distributions. For very
complex seafloors, fitting 3 or more mixture components makes physical sense,
but the statistical model becomes much more complex in these cases. Therefore,
picking the number of components of the mixture model is a decision that must
be made, using a priori information, or using a data driven approach. However,
this information is time consuming to collect, and depends on the skill and
experience of the human. Therefore, a data-driven approach is advantageous to
use, and is explored in this work. Criteria for choosing a model always need
to balance the trade-off for the best fit for the data on the one hand and the
model complexity on the other hand. In this work, we compare several
statistical model selection criteria, e.g., the Bayesian information
criterion. Examples are given for SAS data collected by an autonomous
underwater vehicle in a rocky environment off the coast of Bergen, Norway
using data from the HISAS-1032 synthetic aperture sonar system.

Derek R. Olson | Naval Postgraduate School (NPS), Monterey, California, USA
---|---
Marc Geilhufe | Norwegian Defence Research Establishment (FFI), Kjeller, Norway

##  1 Introduction

Acoustic measurements of seafloor backscattering are a source of unwanted
sound in seafloor object detection 1, 2, 3, but also provide a rich set of
information regarding the seafloor properties and structure 4, 5, 6, 7. The
intensity in a sonar image (i.e. a spatial map of measured backscattering) is
typically characterized by a random process 7. There are a variety of metrics,
or features, that can be used to describe this random process, including the
autocorrelation function, power spectrum 8, wavelet decomposition 9, gray-
level co-occurance matrix 10, 6, the mean intensity (scattering cross section)
11, 12 and in general, the intensity probability density function (pdf) 13,
14, 15.

It was found that for complex scattering environments (such as rocky
seafloors, and man-made structures), a mixture pdf was most appropriate16, 14,
which was justified by the non-stationary character of the acoustic data. Each
sample of the data was modeled as being drawn from a finite number of
distributions, e.g. either from the seafloor or man-made structure, or from
horizontal or vertical facets.

In general, the number of components that make up a non-stationary sonar image
is unknown, and must be selected prior to choosing a model and estimating the
parameters. The more model parameters are used (i.e. more components, or a
more complex statistical model for each component), the better the data will
be fit, but the parameters may lose meaning. In this work, we explore the use
of several model selection techniques based on Bayesian statistics, primarily
the Bayesian information criterion (BIC) and Akaike information criterion
(AIC). These techniques penalize more complex models in different ways. We
also use the log-likelihood (LL) to characterize the model-data fit.

This paper is organized as follows. A description of the sonar data used in
this work and example images are given in Section 2. The background
statistical modeling and model selection techniques are given in Section 3.
Results are presented and discussed in Section 4. Conclusions are given in
Section 5.

##  2 Data

The sonar measurements used in this work are synthetic aperture sonar (SAS)
images collected off the coast of Bergen, Norway by the Norwegian Defence
Research establishment (FFI). The platform used for these measurements is the
HUGIN-HUS autonomous underwater vehicle (AUV), using a HISAS-1032
interferometric SAS. This sonar system has a center frequency of 100 kHz and a
bandwidth of 30 kHz. The beamformed data is oversampled on a grid with
2$\times$2 cm resolution. Data that is used to fit mixture models is decimated
by a factor of 6 in each dimension to reduce the computational complexity of
parameter estimation, while reducing the correlation between samples due to
the point spread function of the sonar system.

Figure 1: An example SAS image, plotted as a function of along track distance
on the horizontal axis and across track distance on the vertical axis. The
color scale is in decibels referenced to an arbitrary pixel pressure, since
the system is uncalibrated.

An example image is shown in Fig. 1. The image consists of an exposed rock
outcrop, with sedimented areas in between. To show the detailed environmental
structure, two tiles are plotted in Fig. 2, both of which are 600$\times$600
pixels, or 12 m per side.

These tiles show that the rock structure consists of a low intensity uniform
scatterer that varies continuously due to undulations in the rock structure.
These continuous variations in intensity are punctuated by bright and dark
lines due to fractures, and steps created by glacial erosion. These features
are distinguishable due to their different intensities, and the SAS system
likely has a high enough resolution that discrete scatterers in the
environment are distinguishable. Therefore a mixture model is appropriate for
modeling the pdf of the ensemble consisting of the pixels from each tile.

Figure 2: Two tiles from Fig. 1, plotted on a decibel scale with 40 dB of
dynamic range. Each image is 600$\times$600 pixels.

##  3 Background

In this section, the basic definitions, statistical models, and model
selection definitions are given. We assume data is given in the form of an
$N_{1}\times N_{2}$ array of intensity values, where the intensity samples are
statistically uncorrelated with each other. The total number of samples is
$N=N_{1}N_{2}$. The random variable for intensity is denoted $I$, and the
amplitude random variable is $A=\sqrt{I}$. We denote the probability of this
variable by $p(a)$, with $a$ being a member of the population.

We use two different distributions to build up a mixture model: a Rayleigh
distribution for amplitude, and a K-distribution. A Rayleigh model for
amplitude pdf has the form

$\displaystyle
p_{R}(a|\lambda_{0})=2a/(\lambda_{0})e^{-a^{2}/\lambda_{0}},a\geq 0,$ (1)

where $\lambda_{0}$ is the mean square value of the pdf. The expected value of
the intensity is $E[a^{2}]=\lambda_{0}$.

A K-distribution for the scattered field amplitude has the form

$\displaystyle
p_{K}(a|\lambda,\alpha)=4/(\sqrt{\lambda}\Gamma(\alpha))(a/\sqrt{\lambda})^{\alpha}K_{\alpha-1}(2a/\sqrt{\lambda}),a>0,$
(2)

where $\Gamma(\cdot)$ is the gamma function, $\lambda$ is the scale parameter
of the K distribution, and $\alpha$ is the shape parameter. The expected value
of the intensity for this model is $E[I]=\alpha\lambda=\sigma_{K}$, where
$\sigma_{K}$ is the average intensity. In the results below, the K
distribution is parameterized using the pair $(\sigma_{K},\alpha)$ rather than
the shape parameter. When used in a mixture distribution, the parameters have
subscripts to denote which K component the parameters correspond to.

Mixture models are formed by a weighted sum of individual pdf components. The
physical meaning of this type of model is that every measurement in a
population, or sample, can be identified with one of the $M$ components. The
weights of the distributions, $w_{m}$ are normalized such that $\sum w_{m}=1$,
and therefore the weights can be interpreted as the fraction of pixels
corresponding to each component.

The form of the mixture models used here is

$\displaystyle
p(a|\theta)=w_{0}p_{R}(a|\lambda_{0})+\sum\limits_{m=1}^{M-1}p_{K}(a|\sigma_{Km},\alpha_{m})$
(3)

where $M$ is the number of mixture components, $\theta$ is a vector of length
$k$, consisting of the parameters of the model consisting of $w_{m}$,
$\lambda$ and the average powers and shape parameters of the K distributions.
The parameters of this mixture model are found using the expectation-
maximization (EM) algorithm 17. This method maximizes a slightly altered
version of the log-likelihood for each component, but asymptotically maximizes
the likelihood function for mixture models 17.

Since the number of components that constitutes the environment is in general
unknown, model selection techniques are used to pick the number of K
distribution components. As the number of components, $M$, increases, it is
better able to match the pdf of the measurement, but at the cost of more
uncertainty per parameter 14. The likelihood function, $\ell(\theta|a)$, is a
common metric for model-data fit. It is defined for N independent samples by,

$\displaystyle\ell(\theta|a)=\prod\limits_{n=1}^{N}p(a_{n}|\theta)$ (4)

where $a_{n}$ is the $n-th$ member of the ensemble, and again $\theta$ is the
parameter vector. The parameter

vector $\hat{\theta}$ that maximizes $\ell$ is called the maximum likelihood
estimate. It is common to work with the log-likelihood,
$\mathcal{L}=\log(\ell)$, which is given as

$\displaystyle\mathcal{L}(\theta|a)=\sum\limits_{n=1}^{N}\log(p(a_{n}|\theta))$
(5)

As stated before, more complex models typically result in a more uncertain
maximum likelihood function, so $\mathcal{L}$ cannot be used as a basis on
which a decision about $M$ can be made. There exist various model selection
techniques, some of which are based on Bayesian concepts18. Two simple metrics
are based on the log-likelihood function, but with an additional penalty that
depends on the number of parameters, $k$. The Bayesian information criterion
(BIC) is defined as

$\displaystyle BIC=-2\mathcal{L}(\hat{\theta})+\log(N)k$ (6)

with smaller $BIC$ preferred. This criterion results from an asymptotic (i.e.
large N) Gaussian approximation of the posterior probability density (ppd)
function of $\theta$ given the data. The penalty to the log-likelihood is
therefore a function of the number of parameters, which is due to the ppd
becoming narrower as $N$ becomes asymptotically large. Another model selection
criterion is the Akaike information criterion (AIC), due to Akaike 19, and
also described by Gelman 18. It is given by the formula

$\displaystyle AIC=-2\mathcal{L}(\hat{\theta})+2k$ (7)

This criterion is simpler than the BIC and penalizes the log-likelihood
independently of the number of samples.

Both the AIC and BIC are based on the log-likelihood function evaluated at
$\hat{\theta}$, and are therefore quite sensitive to the numerical estimate.
Additionally, point estimates do not contain any information about the ppd as
a whole. Other information criteria such as the deviance information criterion
and the Wantanabe-Akaike information criterion18 are not studied here, but we
consider them as fruitful areas for future work. Another possibility is to use
trans-dimensional Monte-Carlo methods to estimate the model with the highest
posterior probability 20.

##  4 Results and Discussion

R-K mixture models were fit to the amplitude data in Figs. 2 (a) and (b),
using between $M=2$ and $M=5$ components. This means that the maximum number
of K-distribution components was 4. The model-data fit is shown graphically in
Fig. 3 for both tiles in terms of the log of the probability of false alarm
(PFA), also called the excedence distribution function (EDF). The PFA is
related to the cumulative distribution function (CDF) through $PFA=1-CDF$, and
is a common method of presenting sonar reverberation statistics 21. The data
in Tile 1 shows a slight “knee” near a normalized amplitude of 2, 6, and 8.
These changes in slope of the log-PFA indicate different components that make
up the model. The log-PFA of Tile 2 has more pronounced knees in the curve,
near the normalized amplitudes of 2, 4, and 12, although the last one is more
uncertain due to the finite sample size being more evident at high amplitudes
(i.e. the PFA curve becomes more stair-case like, rather than a smooth curve).

Model-data fits for both tiles are poor for both R-K1 and R-K2. This behavior
is likely due to the fact that there are not enough components to fit the
data. R-K3 and R-K4 fit the data much better, but are almost the same. This
behavior indicates that a 5-component model does not provide significantly
better fit than the 4-component model. From a visual standpoint, it makes
sense to use a 4-component model for both of these data sets. This hypothesis
will be compared to the results of the more formal model selection techniques.

Model section criteria, the log-likelihood (LL), BIC, and AIC are shown for
both tiles in Table 1. The maximum LL value for Tile 1 is a tie between the 4-
and 5-component models, and the maximum LL for Tile 2 is for the 4-component
model. Intuitively, the more complex model should have a higher likelihood
function, but here, this may be an issue with the numerical parameter
estimates. Additionally, as discussed in Section 3, one cannot base a model
selection decision purely on the maximum of the likelihood function. Further
work on refining these estimate should be made. In terms of the AIC, the
smallest value occurs for the 4-component model for both datasets. However,
for the BIC, the smallest value occurs for the 3-component model for Tile 2,
and the 4-component model for Tile 1. The AIC penalizes model complexity only
slightly, whereas the BIC has a severe penalty for complexity for large number
of data samples. Since the number of data samples was about $10^{4}$, the
penalty for model complexity is about $0.5\log 10^{4}\approx 4.6$ times larger
for the BIC as it is for the AIC. We conclude that the 3-component model is
preferred by BIC for Tile 2 because of the much larger penalty for model
complexity.

Figure 3: The probability of false alarm (PFA) for the data from snippets 1 and 2, compared to the various mixture models explored here. | Tile 1 | Tile 2
---|---|---
Model | LL | AIC | BIC | LL | AIC | BIC
R-K1 | -4321 | 8652 | 8688 | -885 | 1780 | 1810
R-K2 | -4157 | 8331 | 8389 | -847 | 1711 | 1759
R-K3 | -4137 | 8297 | 8377 | -839 | 1700 | 1766
R-K4 | -4137 | 8301 | 8404 | -840 | 1709 | 1793
Table 1: Model selection results for both image tiles. The log-likelihood
(LL), Akaike information criterion (AIC), and the Bayesian information
criterion (BIC) are shown.

##  5 Conclusion

We presented a statistical model for SAS images of complex, non-stationary,
rocky seafloors. This model consisted of a Rayleigh distribution, plus an
unknown number of K distributions. The number of K-distributions was selected
using model selection techniques, including the AIC and BIC. For the AIC, the
4-component model was selected as the most appropriate for two image tiles
used here. For the BIC, different number of components were chosen for each
image tile. It is likely that a different number of components was chosen by
the BIC due to its more severe penalty for model complexity. Future work
should include more robust model selections techniques, such as the deviance
information criterion, which employs a Monte-Carlo Markov chain sampling of
the distribution parameters. This work can also be used to partition an image
into different scattering components, which can aid in estimates of image
complexity, and may also be the basis for quantitative seafloor remote sensing
of geological and/or biological parameters.

##  6 Acknowledgments

This research was supported by grant number N00014-23-WX01149 from the Office
of Naval Research under the Young Investigator Program. The authors thank
FFI’s HUGIN-HUS operator group for collecting the data used in this study.

## REFERENCES

  * 1 D. P. Williams, “Fast unsupervised seafloor characterization in sonar imagery using lacunarity,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 11, pp. 6022–6034, Nov 2015.
  * 2 R. Quinn, _Acoustic Remote Sensing in Maritime Archaeology_ , B. Ford, D. L. Hamilton, and A. Catsambis, Eds. Oxford University Press, sep 2012.
  * 3 A. P. Galusha, J. M. Keller, A. Zare, and G. Galusha, “A fast target detection algorithm for underwater synthetic aperture sonar imagery,” in _Detection and Sensing of Mines, Explosive Objects, and Obscured Targets XXIII_ , J. C. Isaacs and S. S. Bishop, Eds. SPIE, apr 2018.
  * 4 A. Zare, N. Young, D. Suen, T. Nabelek, A. Galusha, and J. Keller, “Possibilistic fuzzy local information c-means for sonar image segmentation,” in _2017 IEEE Symposium Series on Computational Intelligence (SSCI)_. IEEE, nov 2017\.
  * 5 J. Peeples, W. Xu, and A. Zare, “Histogram layers for texture analysis,” _IEEE Transactions on Artificial Intelligence_ , vol. 3, no. 4, pp. 541–552, aug 2022.
  * 6 P. Blondel and O. G. Sichi, “Textural analyses of multibeam sonar imagery from Stanton Banks, Northern Ireland continental shelf,” _Applied Acoustics_ , vol. 70, no. 10, pp. 1288–1297, oct 2009.
  * 7 W. K. Stewart, D. Chu, A. Malik, S. Lerner, and H. Singh, “Quantitative seafloor characterization using a bathymetric sidescan sonar,” _IEEE J. Ocean Eng._ , vol. 19, pp. 599–610, 1994.
  * 8 J. S. Weszka, C. R. Dyer, and A. Rosenfeld, “A comparative study of texture measures for terrain classification,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. SMC-6, no. 4, pp. 269–285, 1976.
  * 9 D. P. Williams, “Unsupervised seabed segmentation of synthetic aperture sonar imagery via wavelet features and spectral clustering,” in _2009 16th IEEE International Conference on Image Processing (ICIP)_. IEEE, nov 2009.
  * 10 R. M. Haralick, K. Shanmugam, and I. Dinstein, “Textural features for image classification,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. SMC-3, no. 6, pp. 610–621, nov 1973.
  * 11 D. R. Jackson and M. D. Richardson, _High-Frequency Seafloor Acoustics_. New York, NY: Springer, 2007\.
  * 12 D. R. Olson, A. P. Lyons, and T. O. Sæbø, “Measurements of high-frequency acoustic scattering from glacially eroded rock outcrops,” _J. Acoust. Soc. Am._ , vol. 139, no. 4, pp. 1833–1847, 2016.
  * 13 A. P. Lyons and D. A. Abraham, “Statistical characterization of high-frequency shallow-water seafloor backscatter,” _The Journal of the Acoustical Society of America_ , vol. 106, no. 3, pp. 1307–1315, sep 1999.
  * 14 D. R. Olson, A. P. Lyons, D. A. Abraham, and T. O. Sæbø, “Scattering statistics of rock outcrops: Model-data comparisons and bayesian inference using mixture distributions,” _J. Acoust. Soc. Am._ , vol. 145, no. 2, pp. 761–774, 2019.
  * 15 R. C. Gauss, J. M. Fialkowski, D. C. Calvo, R. Menis, D. R. Olson, and A. P. Lyons, “Moment-based method to statistically categorize rock outcrops based on their topographical features,” in _OCEANS 2015 - MTS/IEEE Washington_ , Oct 2015, pp. 1–5.
  * 16 D. Abraham, J. Gelb, and A. Oldag, “Background and clutter mixture distributions for active sonar statistics,” _IEEE J. Ocean Eng._ , vol. 36, no. 2, pp. 231–247, April 2011.
  * 17 A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 39, no. 1, pp. 1–22, sep 1977.
  * 18 A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, _Bayesian Data Analysis_. Chapman and Hall/CRC, nov 2013.
  * 19 H. Akaike, “A new look at the statistical model identification,” _IEEE Transactions on Automatic Control_ , vol. 19, no. 6, pp. 716–723, dec 1974.
  * 20 J. Dettmer, S. E. Dosso, and C. W. Holland, “Trans-dimensional geoacoustic inversion,” _The Journal of the Acoustical Society of America_ , vol. 128, no. 6, pp. 3393–3405, dec 2010.
  * 21 A. Lyons, S. Johnson, D. Abraham, and E. Pouliquen, “High-frequency scattered envelope statistics of patchy seafloors,” _IEEE Journal of Oceanic Engineering_ , vol. 34, no. 4, pp. 451–458, oct 2009.


11institutetext: Astrophysics Research Center (ARCO), Department of Natural
Sciences, The Open University of Israel, Raanana 4353701, Israel
22institutetext: Department of Astronomy, Tsinghua University, 30 Shuangqing
Rd, 100084 Beijing, China

# Rocky sub-Neptunes formed by pebble accretion:
Rain of rocks from polluted envelopes

Allona Vazan 11 Chris W. Ormel 22

(November 2022)

Sub-Neptune planets formed in the protoplanetary disk accreted hydrogen-helium
(H,He) envelopes. Planet formation models of sub-Neptunes formed by pebble
accretion result in small rocky cores surrounded by polluted H,He envelopes
where most of the rock (silicate) is in vapor form at the end of the formation
phase. This vapor is expected to condense and rain-out as the planet cools. In
this Letter we examine the timescale for the rainout and its effect on the
thermal evolution. We calculate the thermal and structural evolution of a
10$\,{\rm M}_{\oplus}\,$planet formed by pebble accretion, taking into account
material redistribution from silicate rainout (condensation and settling) and
from convective mixing. We find that the duration of the rainout in sub-
Neptunes is on ${\sim}$Gyr timescale and varies with envelope mass: planets
with envelopes below $\sim$0.75$\,{\rm M}_{\oplus}\,$rainout into a core-
envelope structure in less than $1\,\mathrm{Gyr}$, while planets in excess of
0.75$\,{\rm M}_{\oplus}\,$of H,He preserve some of their envelope pollution
for billions of years. The energy released by the rainout inflates the radius
with respect to planets that start out from a plain core-envelope structure.
This inflation would result in estimates of the H,He contents of observed
exoplanets based on the standard core-envelope structure to be too high. We
identify a number of planets in the exoplanet census where rainout may
operate, which would result in their H,He contents to be overestimated by up
to a factor two. Future accurate age measurements by the PLATO mission may
allow the identification of planets formed with polluted envelopes.

###### Key Words.:

Planets and satellites: formation – Planets and satellites: interiors –
Planets and satellites: composition –

##  1 Introduction

The radius valley (Fulton et al., 2017) splits the population of close-in
exoplanets into two groups: bare rocky objects, and larger radii planets. This
dichotomy may distinguish dry super-Earths from sub-Neptunes, which contains
lighter materials - volatiles and/or hydrogen and helium (hereafters H,He).
While the interior structure of super-Earth planets can be modeled to some
degree from the Earth’s structure, the structure of sub-Neptunes is yet
unclear. Based on the detected radii and masses, close-in sub-Neptunes contain
only a few percent of H,He / volatile. Planet formation models predict more
gas for planets in this mass range (Lee & Chiang, 2015; Mordasini, 2020; Ormel
et al., 2021), where mass loss processes after disk dissipation removed some
of the gas (Owen & Wu, 2013; Gupta & Schlichting, 2020; Rogers & Owen, 2021).
If so, farther-out sub-Neptune planets may have thicker envelopes than the
close-in planets we observe today.

Planet formation models of rocky sub-Neptune planets show that solid
sublimation in the growing gas envelope leads to massive pollution of the
envelope with silicate vapor once a planet exceeds about 2$\,{\rm
M}_{\oplus}\,$(Bodenheimer et al., 2018; Brouwers et al., 2018; Ormel et al.,
2021). The formation of sub-Neptunes by rocky pebble accretion results in a
typical structure of a small rocky core surrounded by a silicate vapor
composition gradient and a vapor-rich convective envelope on top of it (Ormel
et al., 2021). Progressive cooling of these planets lead to over-saturation of
the silicate vapors in the envelope and consequently condensation and settling
(rainout) of rock droplets to deeper undersaturated layers. Rainout sweeps the
silicate from the outer layers downwards and enhances late core growth
(Brouwers & Ormel, 2020). The Rainout phase terminates in a core-envelope
structure, where further cooling won’t change the interior structure. The
timescale for rainout depends on the thermal evolution of the planet. Brouwers
& Ormel (2020) estimated the timescales for rainout in super-Earths to be
billions of years, calculated for uniformly polluted envelopes, based on ideal
gas and adiabatic cooling scheme. Other recent works explored implications of
hydrogen-rock distribution in the interior of sub-Neptune planets, in which
the structure doesn’t evolve into a core-envelope structure (Markham et al.,
2022; Misener & Schlichting, 2022), or that rainout process wasn’t explicitly
modeled (Bodenheimer et al., 2018).

In this Letter we connect planet formation with evolution - we numerically
calculate the post-formation evolution of a 10$\,{\rm M}_{\oplus}\,$sub-
Neptune formed by rocky pebble accretion. We consider a simultaneous thermal
and structural evolution of the interior, which at formation has a silicate
composition gradient and a polluted envelope. In particular, we follow the
rainout process and its consequences for the radius evolution and mass-radius
relation. The 10$\,{\rm M}_{\oplus}\,$case was selected as an intermediate
case to emphasize the difference between interiors of smaller super-Earths and
larger Neptunes. A companion paper describes the evolution of 5-20$\,{\rm
M}_{\oplus}\,$planets formed by rocky pebble accretion under various
conditions (Vazan et al. - in prep.).

##  2 Method

In this section we describe the new methods used in this work which build upon
the interior evolution methods and results of previous works (Vazan et al.,
2015, 2018, 2022).

The initial model is the outcome of planet formation model of rocky pebble
accretion (Ormel et al., 2021). The formation model was calculated for a solid
accretion rate of $10^{-5}$$\,{\rm M}_{\oplus}\,$ yr-1 with both gas and
pebbles contributing to the opacity. At the end of the formation phase the
young 10$\,{\rm M}_{\oplus}\,$planet is embedded in the protoplanetary disk at
0.2AU, and its radius is the Hill radius (about 100$\,{\rm R}_{\oplus}\,$). It
contains 6.7$\,{\rm M}_{\oplus}\,$of metal (silicate rock) and 3.3$\,{\rm
M}_{\oplus}\,$of gas (H,He). More details on the planet formation model appear
in appendix A.

Input parameters are similar to what we used in the planet formation model:
the equations of state of H, He (Saumon et al., 1995), and $SiO_{2}$ rock
(Faik et al., 2018) are the same as in the formation work. Irradiation by the
parent star is for the same location as in the disk phase (i.e., no
migration). Since we don’t expect grains or pebbles to remain in the envelope
in the long term evolution (Brouwers et al., 2021; Movshovitz et al., 2010;
Ormel, 2014; Mordasini, 2014), the radiative opacity is that of a grain-free
solar metallicity atmosphere (Freedman et al., 2014).

The thermal evolution model is calculated for the entire interior - from
center to surface - on one mass grid, with no distinction between core and
envelope. The model solves the interior structure and evolution equations (eq.
1-5 in Vazan et al. (2015)) which allows for heat transport by convection,
radiation and conduction depends on local conditions in time (eq. 11-12 in
Vazan et al. (2015)). The redistribution of composition is by convective
mixing where the convection criterion is fulfilled, and by rainout in
oversaturated regions. Convective-mixing is calculated as described in the
appendix of Vazan et al. (2015). In this work we add the simulation of
composition redistribution by rainout, as explained in the next paragraph.

The Rainout model simulates the condensation and settling of silicate vapor
from an over-saturated layer into an undersaturated layer below, similar to
the approach of Stevenson et al. (2022) for modeling the evolution of Jupiter.
For consistency with the formation models (Ormel et al., 2021) we use the same
liquid-vapor curve of Kraus et al. (2012) (lower branch) to determine the
saturation in the interior:

$\rho_{\rm sat}=\exp\left(a-\frac{b}{T}+\left(\frac{T}{T_{\rm
cp}}\right)^{c}\right)$ (1)

where $a=9.0945$, $b=5.630\times 10^{4}$ K, $c=13.26$, and $T_{\rm cp}=5130$ K
is the critical point of $SiO_{2}$. As the planet cools along the evolution,
whenever the silicate mass fraction in a certain layer $Z_{i}$ exceeds the
saturation level $Z_{\rm i,sat}$, we remove the excess of silicate (above the
saturation level) to the layer below it. Numerically, if $Z_{i}>Z_{\rm i,sat}$
at time $t$ in layer $i$ the silicate mass fraction is:

$Z_{i}(t+1)=Z_{i}(t)-Z_{\rm i,sat}$ (2)

and the excess of rock $(Z_{i}-Z_{sat})m_{i}$ is removed from layer $i$ to the
layer below it, enhancing its metallicity. The rock droplets are assumed to
instantaneous fall between layers, as the droplets settling time is much
shorter than the Kelvin-Helmholtz time. The energy of the rainout - latent
heat of condensation and the settling energy - are included in the evolution
model via the equation of state and the interior structure equations
respectively. For simplicity, we ignore in this work the effect of convective
inhibition in condensation front on slowing the heat transport (Guillot, 1995;
Leconte et al., 2017; Markham et al., 2022).

Mass loss is included in the thermal evolution model, by removal of the mass
from the outermost layer of the planet at a constant rate. Since rainout of
rock droplets from oversaturated layer is fast, mass loss (if it occurs) is of
hydrogen and helium and not of silicates. The rate of mass loss by stellar XUV
depends on assumptions on atmospheric absorption and stellar properties, and
the range of possible outcomes is quite diverse (e.g., Lopez et al., 2012;
Mordasini, 2020; King & Wheatley, 2021; Kubyshkina & Fossati, 2022). We
therefore use constant mass loss rates to achieve mass loss fractions
amounting to 0-0.9 of the initial post-disk envelope mass within the first 0.1
Gyr, during which the stellar XUV is strong.

Figure 1: Silicate mass fraction (left) and temperature profile (right) of a
10$\,{\rm M}_{\oplus}\,$sub-Neptune with 33% H,He located at 0.2 AU. Profiles
are shown at different evolution phases. See appendix A for more details on
the initial (planet formation) profiles.

##  3 Results

###  3.1 Sub-Neptunes with massive envelope from formation

The planet formation phase ends with the dissipation of the disk and rapid
contraction of the outer diluted gas layers from 100$\,{\rm R}_{\oplus}\,$to
about a tenth of that size within a few million years. Despite the significant
radius contraction, the density and composition distribution in the deep
interior is almost not affected. After the rapid contraction, the long term
cooling of the outer metal-rich envelope of this planet is governed by large
scale convection, while the cooling of the deeper composition gradient is by
conduction and/or layered convection, which is less efficient. Consequently,
the metal-rich envelope cools faster and the silicate vapor becomes
oversaturated, while the deep interior remains very hot and thus
undersaturated. As a result, metal starts to condense and settle (rainout)
into deeper layers. As the planet cools the condensation front moves deeper in
the interior. Since the liquid-vapor curve is a strong function of
temperature, the layers above the condensation front are almost metal-free.

Figure 2: Silicate mass fraction (color) as a function of interior layers (y
axis) and time (x axis). Silicate mass fraction ranges between zero (gas only)
in blue and pure silicate in brown. The gradual distribution of silicate from
formation converges into a core-envelope structure after about 4.25 Gyr. The
solid lines signify $Z=0.98$ and $Z=0.02$ enrichment levels. Versions of this
figure for radius and pressure layers instead of mass are presented in
appendix B.

We find that the effect of rainout on the structure of rocky sub-Neptunes
formed by pebble accretion is significant. In figure 1 we show the silicate
(Z) mass fraction distribution (left) and temperature profile (right) of the
10$\,{\rm M}_{\oplus}\,$planet located at 0.2AU. The initial state, based on
the results of the formation model (dashed) and after 0.1, 1, and 5 Gyr
(solid) are shown. Although the planet’s interior at the end of the formation
phase contains most of its silicate as vapor in the polluted envelope, after 5
Gyr all the silicate vapor has condensed and settled and the planet has
transitioned to a core-envelope structure.

Interestingly, the redistribution of silicate in the planet interior is found
to take place on billion years timescale, as is shown in figure 2. The early
fast contraction increases the temperate at the outer layers in the first 10
Myrs, and therefore significant rainout of silicate starts only after about 20
Myr. The rainout continues as the planet cools and ends when all the silicates
are located in a rocky core, at age of 4.25 Gyr for this 10$\,{\rm
M}_{\oplus}\,$sub-Neptune. The simulation shown in figure 2 is for planet with
its initial gas content of 33% in mass (no mass loss). Since the H,He amount
may decrease after disk dissipation, we will examine in section 3.2 how the
rainout timescale changes with envelope mass.

The rainout releases energy from the condensation process (latent heat) and
gravitational settling. Moreover, the settling erases parts of the original
composition gradient, resulting in the release of formation energy that was
stored there, a consequence of low heat transport efficiency in regions
characterized by composition gradients. These energy sources cause the planet
radius to expand for over the duration of the rainout period. When the rainout
terminates the three energy sources vanish and the efficient (convective)
cooling results in faster radius contraction. In the left panel of figure 3 we
present a comparison between the radius of a raining-out planet (red) formed
by pebble accretion to a similar planet with static structure, i.e., no
rainout111The model of rocky plant with polluted envelope and no rainout is
unrealistic, and is shown here in order to emphasis the rainout effect on
radius evolution. (dashed black) and to a similar planet with core-envelope
structure from formation (solid black). As is shown, the effect of rainout on
radius evolution is significant. The luminosity of the planets is shown in the
right panels, indicating significant luminosity increase by the rainout.

Figure 3: Radius (left) and luminosity (right) evolution of the 10$\,{\rm
M}_{\oplus}\,$planet shown in figures 1-2 (red). The evolution of same planets
without rainout (dashed black) and with core-envelope structure from formation
(solid black) are shown for comparison. The green vertical line signifies the
end of the rainout phase. Radius inflation at that point is 7.8% in comparison
to the core-envelope case.

We define the rainout timescale as the time from formation until all the
silicates have settled and a core-envelope structure has been achieved. We
define maximum radius inflation as the radius change at the end of the rainout
phase, in comparison to a similar (mass and composition) planet formed with
core-envelope structure at the same age. In Fig. 3 the end of the rainout
phase, i.e., when the planet reaches a core-envelope structure, is illustrated
by the green vertical dashed line. The rainout timescale is $t=4.25$ Gyr and
the maximum radius inflation is $\Delta R/R=0.078$. As can be seen, the slope
of the radius evolution changes since the planet become fully adiabatic and
cools more efficiently from that point.

The radius inflation of the planet in Fig. 3 is retained in the post-rainout
phase. The energy release from massive envelopes is slow and the radius
evolution is therefore delayed in comparison to planets formed with core-
envelope structure. In other words, planets that formed with polluted massive
envelopes would look ”younger” than planets that formed by core-envelope,
caused by the later rainout heat release. The radius of the raining-out planet
in Fig. 3 is 2.5 Gyr younger along the post-rainout evolution, than a similar
planet that formed with a core-envelope structure. Alternatively, if age is
known this larger radius might be inferred as a higher H,He content for such
observed planet, as will be discussed in the next section.

###  3.2 Sub-Neptunes with lower mass envelopes

The duration of the rainout phase depends on the planet’s cooling efficiency,
hence on the blanketing effect of the envelope. Thicker envelopes will retain
interior heat for longer. In some cases the initial envelope (or part of it)
is lost after the disk dissipates by XUV radiation from the star. The sub-
Neptune shown in figure 2 contains 33% of H,He from formation, but a large
fraction of this envelope can be removed in close-in sub-Neptunes, up to
complete loss (e.g., Owen & Wu, 2013; Rogers & Owen, 2021).

We next explore the rainout process and its duration in simulations that
include mass loss, which we express in terms of the remaining envelope mass.
We repeat the thermal evolution calculation for the same initial 10$\,{\rm
M}_{\oplus}\,$planet with different mass loss rates during the first 0.1 Gyr
of their evolution. In figure 4 we present by the green line the rainout
timescale as a function of the remaining H,He mass. As expected, the smaller
the gas mass (the greater the mass loss) the faster the rainout.
Interestingly, sub-Neptunes with with H,He amounts exceeding 0.75$\,{\rm
M}_{\oplus}\,$ are in the rainout phase for billions of years. They are thus
characterized by inflated radii and gradual / polluted interiors. Gas-poor
sub-Neptunes with less than 0.75$\,{\rm M}_{\oplus}\,$of H,He ”converge” to
core-envelope structure within less than 1 Gyr, resulting in fast (adiabatic)
cooling and contraction in the post-rainout phase. Yet, planets with envelope
masses larger than 0.3$\,{\rm M}_{\oplus}\,$would look ”younger” by their late
released heat content for Gyrs, while planets with less than 0.3$\,{\rm
M}_{\oplus}\,$of H,He cool fast enough that their radius evolution coincides
with the radius evolution of a similar planet that was born with a core-
envelope structure in 1-2 Gyr. In that sense, rocky planets with less than
0.3$\,{\rm M}_{\oplus}\,$H,He that formed with polluted envelopes cannot be
distinguished from planets that formed with a core-envelope structure after
$\sim$1 Gyr.

Figure 4: Green: time from formation until convergence to a core-envelope
structure (rainout timescale) as a function of envelope mass. Trend is shown
for sub-Neptune planets that contain 6.7$\,{\rm M}_{\oplus}\,$of silicates and
their gas (H,He) mass is determined by the mass loss rate, from 3.3$\,{\rm
M}_{\oplus}\,$(original model, no mass loss) down to 0.33$\,{\rm
M}_{\oplus}\,$($3\,M_{\oplus}$ mass loss). Blue: maximum radius inflation by
rainout in comparison to the core-envelope structure model at the rainout
timescale. Curves are polynomial fits for the evolution data points.

Faster rainout in planets with low mass envelopes causes stronger radius
inflation, as similar amount of energy is being released in shorter time.
Thus, the $\sim$8% radius inflation of a planet with massive (3.3$\,{\rm
M}_{\oplus}\,$) envelope at age of 4.25 Gyr becomes 23% radius inflation in a
planet with 0.33$\,{\rm M}_{\oplus}\,$envelope at 0.55 Gyr. In figure 4 we
present (dashed blue) the trend of maximum radius inflation with remaining
envelope mass. As can be seen, the earlier the rainout occurs the larger the
radius inflation is. The early but stronger radius inflation in planets with
low H,He mass might enhance mass loss by photoevaporation, which depends on
the planetary radius, if it overlaps with the strong XUV phase. The longer and
moderate radius inflation in planets with significant H,He envelopes can be
inferred, when using core-envelope model, as a larger H,He content.

##  4 Implications for exoplanets

The envelope (H,He) mass fraction of observed planets are usually derived from
mass-radius relations, assuming a core-envelope structure (e.g., Lopez &
Fortney, 2014; Chen & Rogers, 2016; Baraffe et al., 2008). However, when
inflation by rainout operates the H,He fractions inferred for these planets in
this way would be overestimated. To understand the significance of rainout, we
examine a selection of planets in the exoplanet census with well defined
masses, radii, and ages. Since the radii of these planets are more sensitive
to their H,He fraction than to their planetary mass (Lopez & Fortney, 2014),
we extrapolate our results to slightly larger and smaller planets between
7-12$\,{\rm M}_{\oplus}\,$. For each exoplanet we derived the core-envelope
H,He fraction based on Lopez & Fortney (2014) for standard (solar) opacity. We
then calculated the expected radius inflation by rainout at its age, and used
it to correct the radius and obtain a new H,He fraction, taking into account
the rainout effect. Finally, we run models with the new H,He fraction to
verify their consistency with the observed parameters. As can be seen in Table
1, the actual H,He content can be half of the estimated by the standard model,
when these planets formed with polluted envelopes.

| mass | radius | period | age | H,He | H,He
---|---|---|---|---|---|---
| ($\,{\rm M}_{\oplus}\,$) | ($\,{\rm R}_{\oplus}\,$) | (day) | (Gyr) | standard | w/ rainout
Kepler-36 c | 7.13 | 3.6 | 16.2 | 6.9$\pm 0.372$ | 7.8% | 6%
Kepler-11 e | 7.95 | 4.1 | 32 | $8.5^{+1.1}_{-1.4}$ | 15% | 11%
TOI-1136 d | 7.95 | 4.53 | 12.5 | 0.7$\pm 0.1$ | 13% | 6.5%
TOI-1136 f | 8.3 | 3.8 | 26.3 | 0.7$\pm 0.1$ | 9% | 4.5%
TOI-1422 b | 8.9 | 3.88 | 13 | $5.1^{+3.9}_{-3.1}$ | 10% | 6%
KOI-142 b | 9.54 | 3.37 | 10.9 | $2.4^{+1.2}_{-0.77}$ | 7% | 4%
K2-314 d | 10.2 | 3.86 | 35.7 | 9$\pm 0.6$ | 13% | 10%
K2-19 c | 10.8 | 4.76 | 11.9 | 8+ | 20% | 14%
Kepler-79 b | 10.9 | 3.4 | 13.5 | $3.4^{+0.6}_{-0.91}$ | 6.6% | 4.5%
Kepler-30 b | 11.4 | 3.82 | 29.3 | 2$\pm 0.8$ | 11% | 6%
Table 1: Parameters of observed exoplanets (exoplanet.eu) that we used as
input parameters in the model, and the resulting H,He fraction estimates for
standard core-envelope models and for planets formed with polluted envelopes
and rainout their silicates afterwards. The selected planets are all the
exoplanets with masses 7-12$\,{\rm M}_{\oplus}\,$, radii above 3.3$\,{\rm
R}_{\oplus}\,$, and estimates for stellar age (marked in green in Figure 5).
Error bars were omitted from all observation values and we calculate each
model for a single set of values.

In Figure 5 we show the mass-radius relation of observed exoplanets, and mark
(in green) the planets that we model and are potentially affected by the
rainout, either by post-rainout radius inflation or because they are still
raining-out their polluted envelopes. Planets with radii smaller than
3.3$\,{\rm R}_{\oplus}\,$have less than 5% H,He, thus the rainout effect after
1 Gyr is negligible. Sub-Neptunes with radii above 3.3$\,{\rm
R}_{\oplus}\,$may have underestimated metal content if rainout effect on the
radius is ignored, as is shown in Table 1. Planets above 3.3$\,{\rm
R}_{\oplus}\,$that aren’t marked in green don’t have well defined age estimate
and therefore cannot be interpreted by our model.

Figure 5: Mass-radius relation of observed exoplanets. Planets in green are
potentially inflated by past or ongoing rainout, in comparison to planets that
formed with core-envelope structures. The rainout effect on the inferred H,He
mass of the marked exoplanets appear in Table 1. The sub-Neptune planets with
radii above 3.3$\,{\rm R}_{\oplus}\,$that aren’t marked in green lack an age
measurement and thus are not included in our analysis. Planets with radii
below 3.3$\,{\rm M}_{\oplus}\,$are consistent with low ($<5$%) H,He and thus
rainout has negligible effect on radius at Gyrs age.

Rainout takes place on ${\sim}$Gyr timescale, making the age of the planet an
important parameter in determining its planetary structure. As can be seen
from Table 1, younger sub-Neptunes are prone to larger rainout effects. Since
sub-Neptunes rain out on billion year timescale, which is the age of most of
the observed exoplanets, knowing their age will help us to better constrain
the interior structure of individual observed planets. As of today, many sub-
Neptune planets don’t have well defined age measurement. Future measurements
of stellar age will allow us to suggest corrected H,He fractions for these
planets, for the case they formed with polluted envelopes. Moreover, large
enough statistical database of radius and age will allow us to examine trends
in radius evolution which may hint towards formation mechanisms of sub-
Neptunes. For example, planets formed with polluted envelopes will feature a
spread in radius during their early evolution which is larger than expected by
standard adiabatic cooling. Planets may display a typical radius diversity for
a specific range of H,He mass content. In this context, the expected data from
the PLATO mission (Rauer et al., 2014) on stellar age is expected to drive a
significant improvement in our understanding of planets evolution in general,
and specifically of planets formed with polluted envelopes and on the rainout
process.

Another observational implication is the possible overlap of the rainout phase
with the period of high energy irradiation by the young star. Since mass loss
by photoevaporation strongly depends on planetary radius (third power, Owen &
Wu (2016)), planets with low-mass envelopes and thus early rainout and large
radius inflation, are prone to enhanced ”sedimentation-powered” mass loss. As
we find, a 7$\,{\rm M}_{\oplus}\,$sub-Neptune planet with 5% H,He content will
suffer radius inflation of up to 23%, within the first 0.5 Gyr, which might
double the expected mass loss by photoevaporation at early time. Thus, planets
with less than 0.3$\,{\rm M}_{\oplus}\,$may experience ”sedimentation-powered”
mass loss of their primordial hydrogen atmospheres. A detailed calculation of
the effect of radius inflation on mass loss will be discussed in a companion
study (Vazan et al. - in prep.)

##  5 Conclusions

In this work we find that rainout of high-Z vapor significantly affects the
evolution of sub-Neptunes born with polluted envelopes as it can operate on
evolutionary ($\sim$Gyr) timescales. Timescales we find in this work are
nevertheless shorter than in Brouwers & Ormel (2020), which estimated 10 Gyr
rainout time for super-Earth planets with a few percent H,He. The longer
timescale found in their work is probably the result of relying on the ideal
gas equation of state, which usually suffers from too high temperatures for
the same conditions (see Ormel et al., 2021, for equation of state
comparison).

The key findings from this study are as follows :

  1. 1.

Sub-Neptune planets formed with polluted envelopes experience late growth of a
rocky core due to rainout of high-$Z$ vapor. The timescale for rainout depends
on envelope mass.

  2. 2.

The interior structure of a rocky sub-Neptune formed by pebble accretion
depends on its envelope (H,He) mass and age. Planets with less than
0.75$\,{\rm M}_{\oplus}\,$H,He evolve towards a core-envelope structure by 1
Gyr, while planets with more massive envelopes or younger planets may still be
raining-out when we observe them, and thus have interiors with composition
gradients and/or polluted envelopes.

  3. 3.

Rocky planets with less than 0.3$\,{\rm M}_{\oplus}\,$H,He that formed with
polluted envelopes cannot be distinguished from planets that formed with a
core-envelope structure after about 1 Gyr, due to short rainout phase and
efficient cooling.

  4. 4.

The Rainout process causes radius inflation by release of gravitational energy
(settling), latent heat (condensation), and formation energy (composition
gradient erosion). The radius inflation is shorter and larger in planets with
low mass envelopes: a few percent for sub-Neptunes with massive gas envelope,
and larger (tens of percent) but on much shorter timescale in planets with
lower gas mass fraction.

  5. 5.

Planets that formed with polluted envelopes would look ”younger” than planets
that formed with core-envelope structure, as a result of the later heat
release. Alternatively, if age is known, the inflated radii might be
interpreted as higher H,He mass in these planets, which can be up to twice of
the actual value.

The model in this Letter is for rocky sub-Neptunes, formed as 10$\,{\rm
M}_{\oplus}\,$planet, with and without envelope mass loss. Significantly more
(less) massive planets have stronger (weaker) gravity fields, which will
affect the results quantitatively. In a companion paper (Vazan et al. - in
prep.) we execute parameter space study for a wider range of masses and
formation parameters.

## Acknowledgments

We thank the referee for constructive comments that improved the manuscript.
A.V. acknowledges support by ISF grants 770/21 and 773/21. C.W.O. acknowledge
support by the National Natural Science Foundation of China (grant no.
12250610189). Figures were plotted using Matplotlib (Hunter, 2007).

## References

  * Baraffe et al. (2008) Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
  * Bodenheimer et al. (2018) Bodenheimer, P., Stevenson, D. J., Lissauer, J. J., & D’Angelo, G. 2018, ApJ, 868, 138
  * Brouwers & Ormel (2020) Brouwers, M. G. & Ormel, C. W. 2020, A&A, 634, A15
  * Brouwers et al. (2021) Brouwers, M. G., Ormel, C. W., Bonsor, A., & Vazan, A. 2021, A&A, 653, A103
  * Brouwers et al. (2018) Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65
  * Chen & Rogers (2016) Chen, H. & Rogers, L. A. 2016, ApJ, 831, 180
  * Faik et al. (2018) Faik, S., Tauschwitz, A., & Iosilevskiy, I. 2018, Computer Physics Communications, 227, 117
  * Freedman et al. (2014) Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109
  * Guillot (1995) Guillot, T. 1995, Science, 269, 1697
  * Gupta & Schlichting (2020) Gupta, A. & Schlichting, H. E. 2020, MNRAS, 493, 792
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * King & Wheatley (2021) King, G. W. & Wheatley, P. J. 2021, MNRAS, 501, L28
  * Kraus et al. (2012) Kraus, R. G., Stewart, S. T., Swift, D. C., et al. 2012, Journal of Geophysical Research (Planets), 117, E09009
  * Kubyshkina & Fossati (2022) Kubyshkina, D. & Fossati, L. 2022, A&A, 668, A178
  * Leconte et al. (2017) Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, A&A, 598, A98
  * Lee & Chiang (2015) Lee, E. J. & Chiang, E. 2015, ApJ, 811, 41
  * Lopez & Fortney (2014) Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1
  * Lopez et al. (2012) Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59
  * Markham et al. (2022) Markham, S., Guillot, T., & Stevenson, D. 2022, A&A, 665, A12
  * Misener & Schlichting (2022) Misener, W. & Schlichting, H. E. 2022, MNRAS, 514, 6025
  * Mordasini (2014) Mordasini, C. 2014, A&A, 572, A118
  * Mordasini (2020) Mordasini, C. 2020, A&A, 638, A52
  * Movshovitz et al. (2010) Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616
  * Ormel (2014) Ormel, C. W. 2014, ApJ, 789, L18
  * Ormel et al. (2021) Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175
  * Owen & Wu (2013) Owen, J. E. & Wu, Y. 2013, ApJ, 775, 105
  * Owen & Wu (2016) Owen, J. E. & Wu, Y. 2016, ApJ, 817, 107
  * Rafikov (2006) Rafikov, R. R. 2006, ApJ, 648, 666
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249
  * Rogers & Owen (2021) Rogers, J. G. & Owen, J. E. 2021, MNRAS, 503, 1526
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Stevenson et al. (2022) Stevenson, D. J., Bodenheimer, P., Lissauer, J. J., & D’Angelo, G. 2022, PSJ, 3, 74
  * Vazan et al. (2015) Vazan, A., Helled, R., Kovetz, A., & Podolak, M. 2015, ApJ, 803, 32
  * Vazan et al. (2013) Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283
  * Vazan et al. (2018) Vazan, A., Ormel, C. W., Noack, L., & Dominik, C. 2018, ApJ, 869, 163
  * Vazan et al. (2022) Vazan, A., Sari, R., & Kessel, R. 2022, ApJ, 926, 150

##  Appendix A Planet formation model

The input model we use for the evolution is the outcome of the planet
formation model of Ormel et al. (2021), which calculates the growth of the
core and envelope of an embedded, pebble-accreting planet. The envelope is
assumed to be in hydrostatic equilibrium on dynamical timescales and in
pressure equilibrium with the surrounding disk gas. At each time step a steady
solution to the stellar structure equations is found. A key element in the
model is the relaxation of the standard assumption that accreting solids end
up in the core.

The formation calculation assumes constant solid accretion rate of
$10^{-5}$$\,{\rm M}_{\oplus}\,$ yr-1. The effects of pebble accretion on gas
accretion are included via acceleration of gas accretion through accumulation
of a vapor-rich atmosphere, and by contribution to the opacity, which reduces
the intake of nebular gas. The formation model uses tabular equation of states
of H,He and $SiO_{2}$ and their mixtures, based on Vazan et al. (2013). The
vaporization is determined by the silicate liquid-vapor (saturation) curve of
Kraus et al. (2012).

In figure 6 we show the initial model we use in this work. The model is
identical to what appears in the middle-bottom panel of figure 13 in Ormel et
al. (2021). It is the outcome of the growth of a planet embedded in the
protoplanetary disk at 0.2 AU, starting from an embryo of 0.06$\,{\rm
M}_{\oplus}\,$and proceeding to the point where 10$\,{\rm M}_{\oplus}\,$mass
is reached after about 5 Myr. The solid core in the planet formation model is
not part of the interior grid and is modeled by constant density and arbitrary
temperature. These input density and temperature profiles converge to
isentropic profiles after the initial procedure in the evolution model, in
which realistic equation of state and structure equations are applied also to
the core.

The end of the planet formation phase is chosen at the state where the planet
mass reached 10$\,{\rm M}_{\oplus}\,$. As expected, if the simulation is
allowed to continue in gas-rich vicinity the planet would reach the crossover
point and start its runaway gas accretion. Nevertheless, pre-runaway
evolutionary timescales typically slow down (e.g., Rafikov 2006; Lee & Chiang
2015; Ormel et al. 2021, figure 6), making the formation of such 10$\,{\rm
M}_{\oplus}\,$planet with a significant H,He envelope conceivable.

Figure 6: The initial model we use for the evolution calculations, resulting
from planet formation model of Ormel et al. (2021). Metallicity (black) and
temperature (red) profiles for a planet formed by rocky pebble accretion at
0.2 AU. The temperature in the classical core region that is not modeled in
the formation model, is fitted by isentropic profile in the evolution model,
see text for details.

##  Appendix B Rainout vs. radius and pressure

The data shown in figure 2 as a function of the code’s mass grid (y-axis) may
give the impression that rainout takes place in the outer atmosphere. We
therefore plotted the same data also as a function of radius (figure 7) and
pressure (figure 8) in the y-axis. Solid lines indicate values of $Z=0.02$ and
$Z=0.98$. As can be seen in the two figures, the envelope pollution
(evaporation front) is located in the deep interior at pressure $>0.1$ GPa, at
about half of the total radius. While rainout takes place deep in the
envelope, the outer envelope (and atmosphere) are left almost silicate-free,
and thus prone to easier inflation by the rainout heat release.

Figure 7: Same data of figure 2, here as a function of radius instead of mass
(y axis). Silicate mass fraction (color) ranges between zero (gas only) in
blue and pure silicate in brown. White area is empty, due to radius
contraction.

Figure 8: Same data of figure 2 and 7, here as a function of pressure (y
axis). Silicate mass fraction (color) ranges between zero (gas only) in blue
and pure silicate in brown. White areas are empty, due to pressure increase in
contraction.