# Constraining the effect of convective inhibition on the thermal evolution of
Uranus and Neptune

Steve<EMAIL_ADDRESS>California Institute of Technology, Dept.
of Geological and Planetary Sciences  Dave Stevenson  California Institute of
Technology, Dept. of Geological and Planetary Sciences

###### Abstract

The internal heat flows of both Uranus and Neptune remain major outstanding
problems in planetary science. Uranus’ surprisingly cold effective temperature
is inconsistent with adiabatic thermal evolution models, while Neptune’s
substantial internal heat flow is twice its received insolation. In this work
we constrain the magnitude of influence condensation, including latent heat
and inhibition of convection, can have on the thermal evolution of these
bodies. We find that while the effect can be significant, it is insufficient
to solve the Uranus faintness problem on its own. Self-consistently
considering the effects of both latent heat release and stable stratification,
methane condensation can speed up the cool down time of Uranus and Neptune by
no more than 15%, assuming 5% molar methane abundance. Water condensation
works in the opposite direction; water condensation can slow down the cool
down timescale of Uranus and Neptune by no more than 15% assuming 12% molar
water abundance. We also constrain the meteorological implications of
convective inhibition. We demonstrate that sufficiently abundant condensates
will relax to a state of radiative-convective equilibrium requiring finite
activation energy to disrupt. We also comment on the importance of considering
convective inhibition when modeling planetary interiors.

##  1 Introduction

Giant planet atmospheres are primarily heated by a combination of sunlight and
internal heat leftover from formation. All giant planets except Uranus are
observed to emit more infrared radiation into space than the absorbed
sunlight, by approximately a factor of two. Jupiter’s present day luminosity
can be approximately explained by convective cooling from an initially hot
state over the age of the solar system [Hubbard, 1977][Hubbard et al., 1999].
To accurately reproduce Saturn’s present day state, one may need to account
for additional heating by the settling of helium rain from the envelope into
the interior [Hubbard et al., 1999][Stevenson, 1983]. However, luminosity is a
crude indicator of thermal evolution since planets can store heat internally
and may have internal heat sources (e.g. differentiation). The present
luminosities of Uranus and Neptune are not well understood because even their
basic structures, including composition, internal structure, and thermal
transport properties, are not well understood.

Measurements of the ice giants’ electromagnetic emission to space began in the
1960s [Kellermann and Pauliny-Toth, 1966], with high quality far infrared
measurements constraining the effective temperatures beginning in the 1970s
[Fazio et al., 1976][Loewenstein et al., 1977][Stier et al., 1978]. These
early observations concluded that Uranus appeared approximately in equilibrium
with its received sunlight, while Neptune emitted more than twice the
radiation it received. These observations were corroborated by higher quality
analysis after the Voyager 2 flybys [Pearl et al., 1990][Pearl and Conrath,
1991]. The 1$\sigma$ upper limit for Uranus’ energy balance (the ratio between
its emitted and absorbed thermal flux) is 1.14. The lower limit is below
unity, indicating the results are consistent with zero internal heat flow.
However, we know the heat flow cannot be zero because Uranus has a magnetic
field. Moreover, the higher microwave temperatures at long wavelengths (e.g.
[Gulkis et al., 1983]) are compatible with heat flow from depth. Uranus must
be convective at depth.

Theoretical attempts to explain these observations began promptly. It was
immediately clear that the ice giants could not have the same thermal
histories as the gas giants. Early studies concluded that, if these planets
cool convectively like the gas giants, they must have formed at a temperature
not much warmer than their current states [Hubbard, 1978][Hubbard and
MacFarlane, 1980], a highly unlikely interpretation because the energy of
accretion $\sim GM^{2}/R$ far exceeds their current heat content for any
plausible assumption of structure. Alternative theories suggested a large
fraction of gravitational heat of formation remains trapped in the interior,
but by some mechanism cannot escape to space [Podolak et al., 1991]. More
recent studies suggest that there is no problem for Neptune [Fortney et al.,
2011][Linder et al., 2019], or even that Neptune’s present luminosity is
higher than expected [Nettelmann et al., 2016][Scheibe et al., 2019]. Uranus’
very low internal heat flux, sometimes known as the faintness problem [Helled
et al., 2020], remains largely unsolved, although it has been suggested that
the problem can be solved by modeling thin layers of static stability near
phase boundaries [Nettelmann et al., 2016]. Today it is largely accepted that
the adiabatic assumption for the interior is probably inappropriate for Uranus
and Neptune [Helled et al., 2020].

In this work we present a mechanism that inhibits convection near the methane
cloud level, thereby trapping internal heat beneath the clouds. This mechanism
has already been theorized and discussed e.g. [Leconte et al., 2017] [Friedson
and Gonzales, 2017] [Guillot, 2005], but the effect of methane on the ice
giants has not yet been explicitly quantified and worked into a thermal
evolutionary model. In hydrogen atmospheres, sufficiently abundant condensible
species can shut off convection near the cloud level [Guillot, 1995][Guillot,
2005]. By “sufficiently abundant” we mean greater than an analytically
calculable critical mole fraction $q_{\text{crit}}$. This value is about 1.4%
for methane and 1.2% for water under the relevant conditions in Uranus and
Neptune. Recent theoretical study confirms this effect is also stable against
double diffusive convection in a saturated medium in the fast precipitation
limit, indicating radiation would be the only remaining efficient thermal
transport mechanism [Leconte et al., 2017] [Friedson and Gonzales, 2017].
Methane is certainly sufficiently abundant for convective inhibition to occur
[Helled et al., 2020]. The long-term survival of the configuration against
entrainment is still a subject of research; the configuration may be
intermittently eroded, destroyed, and reformed [Friedson and Gonzales, 2017].

In Section 2 we begin by laying the heuristic groundwork and providing
analytic order of magnitude estimates of the effect of convective inhibition.
Then in Section 3 we outline a more detailed atmospheric model. We then use
this model to quantify the effects of condensation on the difference between
the planet’s observed effective temperature and its internal entropy. We also
comment on the meteorological implications arising from convective inhibition
on Uranus, Neptune and Saturn, as well as its importance for interior
modeling. In Section 4 we constrain the importance of both methane and water
condensation on the planets’ thermal evolution. Finally in Section 5 we make
recommendations for future missions to the ice giants, and comment on
additional applications of this mechanism to the thermal histories of
exoplanets, especially super-Earths.

##  2 Intuition and analytic approximations

A hydrogen atmosphere becomes stable against convection at a critical value of
the condensate mole fraction $q_{\text{crit}}$ that depends on temperature and
the properties of the condensate and the gas mixture [Guillot, 1995] [Guillot,
2005] [Li and Ingersoll, 2015] [Leconte et al., 2017] [Friedson and Gonzales,
2017]. The mechanism is as follows: consider an isobaric open system hydrogen
gas parcel saturated with a vapor species of higher molecular weight. Assume
there exists a finite reservoir of liquid condensate in equilibrium with the
saturated parcel, outside but in contact with the system. If the parcel is
relatively cool, the effect of the condensate will be a small correction, and
the parcel will approximately behave like an ideal gas such that density
decreases as temperature increases. However, as temperature increases at fixed
pressure, the mixing ratio of the condensate likewise increases. Because the
condensate vapor is heavier than the dry air, there comes a crossover where
the Arrhenius relationship in temperature governing vapor pressure saturation
overcomes the linear relationship in temperature governing mean spacing
between molecules in a gas. After this crossover point for the system outlined
above, increasing temperature actually increases the density of the parcel.
For this reason, a hydrogen atmosphere with sufficiently abundant condensate
($q_{\text{max}}>q_{\text{crit}}$) with an internal heat source will not
convect. It does not convect because the warmer underlying gas is Ledoux
stable to convection for any temperature gradient if saturated at all levels.
In our notation,

$q_{\text{crit}}^{-1}=\left(\frac{\mu_{c}L}{RT}-1\right)(\epsilon-1)$ (1)

where $R$ is the ideal gas constant, $L$ is the latent heat of vaporization,
$\mu_{c}$ is the molecular weight of the condensate, and
$\epsilon\equiv\mu_{c}/\mu_{d}$ is the ratio of the condensate molecular mass
to that of dry air. Note that our definition of $q$ is the molar mixing ratio,
and is different from the quantity defined as $q$ in [Leconte et al., 2017].

If we neglect radiative transfer of heat, we can analytically approximate the
temperature difference between the top and bottom of the stable layer. The
bottom of the upper saturated level satisfies $q=q_{\text{crit}}$, while the
deep well-mixed convective atmosphere has uniform composition of the vapor
satisfying $q=q_{\text{max}}$. In the limit of no thermal transport, the
stable layer will reduce to a stable interface, with an unsaturated convective
level beneath a saturated convective level. Thus, we can approximate the top
of the deep well mixed convective layer and the bottom of the upper saturated
convective layer to be at the same pressure level $p$. We define the
temperature at the bottom of the saturated convective layer to be $T_{1}$, and
the temperature at at the top of the well-mixed convective layer to be
$T_{2}$. In this case, if the saturation vapor pressure $p_{s}$ satisfies the
Clausius-Clapeyron relationship where $q=p_{s}/p$, we can immediately write
down

$\frac{1}{T_{2}}-\frac{1}{T_{1}}=\frac{R}{\mu_{c}L}\log\left(\frac{q_{\text{crit}}}{q_{\text{max}}}\right)$
(2)

Equation 2 turns out to be a good approximation for very deep clouds (pressure
of order a hundred bars, e.g. the water cloud level in the contemporary ice
giants), where the atmosphere is relatively opaque. This approximation is less
accurate when the atmosphere is less opaque and radiative transfer is more
efficient, such as the methane cloud deck of the contemporary ice giants, or
the water cloud deck earlier in their thermal histories. Equation 2 will also
always be an upper bound to the difference between $T_{1}$ and $T_{2}$ because
it neglects thermal transport. In practice the difference between the pressure
levels at the top and bottom of the stable layer play an important role. In
Section 3, we solve the problem of radiative convective equilibrium
explicitly.

The goal of this analysis is to quantify the effect of condensation on a
planet’s effective temperature. For a given internal entropy, there will be
three important theoretical effective temperatures we discuss throughout this
paper: $T_{e}$, $T_{\text{int}}$, and $T_{ab}$. We define these quantities
here, referring to Figure 1 to illustrate a sample atmospheric profile that
corresponds to its respective temperature. The effective temperature $T_{e}$
(corresponding to the solid red curve on Figure 1) accounts for convective
inhibition, and should correspond to the observed effective temperature of the
planet from the outside. The internal effective temperature $T_{\text{int}}$
(faded solid blue curve on Figure 1) accounts for latent heat but not
convective inhibition. Finally the adiabatic effective temperature $T_{ab}$
(black curve on Figure 1) is the effective temperature the planet would have
if no condensation occurred at all, i.e. if the whole troposphere were dry
adiabatic.

Figure 1: Sample atmospheric profile useful for intuition. Blue curves show
wet adiabatic (wet ad) temperature (solid) and virtual temperature (dashed)
profiles. Red curves show our radiative-convective equilibrium (RCE)
solutions, and the black curve shows the dry adiabatic reference solution, all
corresponding to the same interior entropy. Magenta and purple curves
correspond to saturated abundance of methane for RCE and wet ad respectively,
while the corresponding dashed curves show the local critical mixing ratio
$q_{\text{crit}}$.

We now approximate how this temperature difference between $T_{1}$ and $T_{2}$
changes the effective temperature of the planet. Consider an initially
convective atmosphere (for example, mixed by a cosmic ladle) that has
effective temperature $T_{\text{int}}$. If $q_{\text{max}}>q_{\text{crit}}$,
then part of this atmosphere will be stable to convection. The atmosphere will
cool from the top, but cannot carry that heat out convectively, causing the
upper layer to relax onto a cooler adiabat until radiative-convective
equilibrium is reached. We wish to estimate the difference between the initial
effective temperature $T_{\text{int}}$ and the final effective temperature
$T_{e}$ in equilibrium. Entropy increment scales as $dS\propto dT/T$ in an
isobaric environment, according to the first law of thermodynamics, and
adiabatic processes are isentropic. Therefore if the difference between
$T_{1}$ and $T_{2}$ is not too large, we can approximate

$\frac{T_{e}-T_{\text{int}}}{T_{\text{int}}}\sim\frac{T_{1}-T_{2}}{T_{2}}\sim\frac{RT_{2}}{\mu_{c}L}\log\left(\frac{q_{\text{crit}}}{q_{\text{max}}}\right)$
(3)

where $T_{e}$ is the observed effective temperature of the planet after
accounting for convective inhibition. This approximation neglects the non-
adiabaticity due to latent heat release, which we address in the following
paragraph.

In reality we expect two effects arising from condensation: the tendency
toward sub-adiabatic wet pseudo-adiabaticity arising from latent heat, and the
tendency toward super-adiabatic stable stratification arising from convective
inhibition. These two effects will be in opposite directions; while latent
heat will tend to produce a warmer effective temperature, convective
inhibition will tend to produce a cooler effective temperature with the same
interior entropy. The former effect has been studied in detail [Kurosaki and
Ikoma, 2017], while the latter effect is the subject of this work. We can
estimate the magnitude of the latent heat effect using the definition of
equivalent potential temperature that is conserved along a moist adiabat.

$\theta_{e}(p,T)=\theta(p,T)\exp\left[\frac{\epsilon
p_{s}(T)}{c_{p}T(p-p_{s}(T))}\right]$ (4)

where $\theta\equiv T(p_{0}/p)^{\nabla_{ab}}$ is the potential temperature,
$\nabla_{ab}\equiv\frac{\gamma-1}{\gamma}$ is the adiabatic gradient, and
$\gamma$ is the Grüneisen parameter. Using this, we can estimate the
temperature change due to moist adiabaticity as
$\frac{T_{e}-T_{\text{int}}}{T_{\text{int}}}\sim\exp\left[\frac{\nabla_{ab}\mu_{c}q_{\text{max}}}{RT_{2}}\right]$.
Comparing this to Equation 3 demonstrates these two quantities should not
scale in the same way. However, they are comparable in order of magnitude
under the conditions of interest. Therefore both effects must be accounted for
explicitly in order to fully understand the effect of condensation on thermal
evolution of planets with polluted hydrogen atmospheres. We perform this
calculation in Section 3.

The importance of these effects on thermal evolution are as follows. Because
potential temperature relates linearly to a reference temperature, the
temperature at all pressures will scale linearly with the temperature at some
reference pressure. Likewise, the effective temperature of a planet scales
linearly with a reference temperature in the adiabatic region, assuming
constant opacity (this is actually a poor assumption, and a fully complete
model must include opacity variations due to condensation explicitly. See
Section 5 for further details). Therefore, it is possible to model thermal
evolution by assuming a planet’s effective temperature is linearly related to
its internal heat content. The purpose of calculating the difference between
the observed effective temperature $T_{e}$ and the adiabatic equivalent
effective temperature $T_{ab}$ is to explicitly quantify the non-linearity
arising from condensation, so that thermal evolution can be modeled self-
consistently. We carry this out in Section 4.

##  3 Atmospheric model

We model a radiative-convective equilibrium atmosphere using a two stream gray
opacity approximation for thermal radiative transfer. We seek to uniquely
define the apparent effective temperature $T_{e}$ as a function of a planet’s
internal equivalent effective temperature $T_{\text{int}}$ and condensate
abundance $q_{\text{max}}$. Planetary and physical properties, such as surface
gravity and the physical properties of the gas mixture, are considered to be
fixed.

We also assume the planet is subject to intermittent moist convective events
that overcome the potential barrier of the stable layer. These could occur due
to instabilities caused by entrainment over long timescales [Friedson and
Gonzales, 2017], rare impact events, or strong updrafts from the interior. The
equilibrium configuration then is reached by gradual cooling, with the upper
layer relaxing onto a moist pseudo-adiabat set by a different potential
temperature than the adiabat that sets the interior. The stable layer
meanwhile will have a super-adiabatic temperature gradient set by thermal
radiative equilibrium, see Figure 2.

In this section we explore the effect of methane abundances varying between
2-5%. The most commonly cited number for the deep mixing ratio of methane in
Uranus and Neptune are 2.3% and 2% respectively, because these are the nominal
values in the first published work on the atmospheric structure of these
planets derived from radio refractivity data from Voyager [Lindal et al.,
1987] [Lindal, 1992]. However, these early works provide solutions to the
refractivity data using assumptions of methane abundance between 1% and 4%,
finding all these solutions to be theoretically compatible with the
observations. Subsequent analysis from ground based and Hubble observations
have likewise found a range of acceptable values for both planets ranging
between roughly 2%-4% for both planets, as well as latitudinal variation in
methane mixing ratio [Baines et al., 1995][Rages et al., 1991] [Baines and
Hayden Smith, 1990] [Karkoschka and Tomasko, 2011]. In this work we are
interested in understanding how thermal evolution is affected by methane
condensation. Since the exact mixing ratio is not precisely constrained and we
are interested in this question broadly, we take sample values for methane
concentration between 2-5% to understand how the effect changes with methane
abundance. Unsurprisingly, the effect becomes monotonically more important as
concentration increases within this range, as shown in Figures 3 and 6.

###  3.1 Defining the boundaries of the stable layer

Figure 2: A schematic sketch of the model (not to scale). The red curves
represent temperature increasing to the right, while black horizontal lines
demarcate important pressure levels. The convective layers are labeled
$\theta_{e}(T_{e})$ and $\theta(T_{\text{int}})$ respectively, to indicate
that their temperature structure is dictated by setting $\theta_{e}$ and
$\theta$ to a constant value uniquely determined by $T_{e}$ and
$T_{\text{int}}$.

In this subsection we quantify important pressure boundaries we need to define
the radiative transfer model in the following subsection. We are interested in
the case where condensate is sufficiently abundant to inhibit convection,
$q_{\text{max}}>q_{\text{crit}}$. We consider an atmosphere where optical
depth unity in the IR is at lower pressure than the level at which convection
is expected. At deeper levels (higher pressure) we assume there is a region of
rapidly varying condensate mixing ratio in the vapor phase; this region can be
convectively stable as discussed in Section 2. Deeper still, below the
conventionally defined cloud deck, the condensate mixing ratio is a constant
because the vapor pressure is always less than the saturated vapor pressure at
that temperature. We refer to this below as the “bulk mixing ratio” though it
is strictly only applicable to whatever deep, well mixed layer lies beneath
the clouds and says nothing about the actual methane abundance at far deeper
levels (i.e., the methane abundance of the planet as a whole). Accordingly,
our atmosphere has (from the top downward) a radiative layer (the
stratosphere) a convective layer, another radiative layer (called the ‘stable
layer’ below) and a deep convective layer.

Assuming a two stream approximation with collimated light beams, the thermal
structure of an atmosphere in radiative equilibrium as a function of optical
depth is $T(\tau)=T_{e}(\tau+1/2)^{1/4}$. We assume the IR opacity to be
dominated by pressure-induced opacity of hydrogen collisions that
approximately obeys $\kappa\sim\kappa_{0}(p/p_{0})$, where
$\kappa_{0}=10^{-2}$g-1cm2 and $p_{0}=1$ bar. Assuming a different
$\kappa_{0}$ does not significantly affect our findings, despite altering some
of the model details. An atmosphere in radiative equilibrium becomes unstable
to convection at the point where its lapse rate becomes superadiabatic. Under
our assumptions, this places the pressure level of the radiative-convective
boundary

$p_{\text{rc}}^{2}=\frac{4gp_{0}R_{d}}{(c_{p}-2R_{d})\kappa_{0}}$ (5)

where $R_{d}$ is the specific gas constant of dry air, and $c_{p}$ is its
constant pressure heat capacity. Note for ideal gasses
$R_{d}/c_{p}=\nabla_{ab}$. We use a dry adiabatic lapse rate to set the
radiative-convective boundary because we assume the effect of moist
adiabaticity is small in this relatively cold part of the atmosphere. Beneath
the boundary, we assume the atmosphere to be moist adiabatic. A moist
adiabatic atmosphere conserves the equivalent potential temperature, Equation
4. We set the moist adiabatic equivalent potential temperature using the
temperature and pressure at the radiative-convective boundary. By doing this,
we define a unique moist adiabat for a given effective temperatures.

In equilibrium, the mole mixing ratio is set by the condensate’s saturated
vapor pressure $q(p,T)=p_{s}(T)/p(T)$. When we define a moist adiabat, the
temperature is uniquely defined at every pressure level. Therefore we can
solve for the level $p_{1}$ where the atmosphere becomes stable to convection
by solving $q(p_{1})=q_{\text{crit}}(T(p_{1}))$, where $q_{\text{crit}}$ is
defined in Equation 1.

Similarly, we can solve for the bottom of the stable layer by solving
$q(p_{2})=q_{\text{max}}$, where $q_{\text{max}}$ is the bulk abundance of the
condensate species. This is set using the pseudoadiabat corresponding to a
planet with effective temperature that neglects convective inhibition. By
fixing the uninhibited effective temperature $T_{\text{int}}$ that defines
$p_{2}$, we can solve for the corresponding effective temperature $T_{e}$ that
satisfied radiative convective equilibrium, i.e.
$F_{\uparrow}(p_{2})-F_{\downarrow}(p_{2})=F_{\uparrow}(p_{1})-F_{\downarrow}(p_{1})=\sigma(T_{e}^{4}-fT_{0}^{4})$
where $f$ is the fraction of sunlight absorbed above the stable layer, and
$T_{0}$ is the equilibrium effective temperature with the sun in the absence
of internal heat (see Figure 2).

Figure 2 contains apparent temperature discontinuities, which exist in the
model. Of course, temperature discontinuities are not stable in natural media,
as conductive heat transport will be infinite. Additionally, a temperature
discontinuity–even if stabilized by a compositonal difference–will lead to
negligible temperature differences due to thin thermal boundary layer
convection in an inviscid fluid. The temperature discontinuities in Figure 2
therefore do not actually represent discontinuities in nature, but steep
temperature gradients. We provide the following order of magnitude analysis to
determine how important these steep quasi-discontinuities are in the context
of the model. The relaxation timescale for thermal diffusion in a medium with
thermal conductivity $k_{t}$ is $\tau_{c}\sim\frac{\rho c_{p}D^{2}}{k_{t}}$,
where $D$ is the vertical length scale of relevance. The timescale for
radiative relaxation is the radiative time constant
$\tau_{\text{rad}}\sim\frac{c_{p}}{8\sigma T^{3}\kappa}$. We solve for the
thickness of the conductive layer by equating the two timescales solving for
$D$. Using appropriate parameters for hydrogen around 1 bar ($k_{t}\sim
10^{4}$g cm s-3 K-1, $\rho\kappa\sim 10^{-6}$) we find the length scale to be
of order 10-100 meters, small compared to atmospheric length scales (i.e., the
scale height). Then using Fourier’s Law, we find heat conduction to be of
order $10^{-2}$erg cm-2 s-1 for discontinuities of order 1K (scaling linearly
with the size of the temperature discontinuity), about four orders of
magnitude smaller than $\sigma T_{e}^{4}$ and therefore not included in the
model. Thus the discontinuities in Figure 2 are really there in the model, but
are physically understood to be steep temperature gradients nevertheless
unimportant for the purposes of calculating total heat flow.

We must also comment on entrainment by convection outside the stable layer.
Entrainment will tend to erode and thin the stable layer over time [Friedson
and Gonzales, 2017]. In general for water, the erosion timescale is greater
than the cooling timescale, indicating that the equilibrium configuration
should exist at some times. However as the stable layer is eroded and becomes
thinner, heat transport across the stable layer will be enhanced, reducing the
difference between $T_{e}$ and $T_{\text{int}}$. Eventually the thinning
stable region will reduce to an interface that may be stable, unstable, or
conditionally stable. We acknowledge that these complications are confounding
factors for our model, and therefore our results that neglect entrainment
erosion of the stable layer should be thought of as an upper bound on the
magnitude of the effect on $\Delta T\equiv T_{e}-T_{\text{int}}$ and on
evolution.

###  3.2 Radiative transfer across the stable layer

Figure 2 is a useful visual reference for this section. In order to compute
the radiative-convective equilibrium solution, we first solve for the
equilibrium heat flow for a system specifying the boundaries of the stable
layer $p_{1}$ and $p_{2}$, along with their corresponding temperatures $T_{1}$
and $T_{2}$. The temperature structure above $p_{1}$ is moist adiabatic, while
the temperature structure below $p_{2}$ is dry adiabatic. In equilibrium
$\frac{dT}{dt}\propto\frac{dF}{dp}=0$. Then using the two stream approximation
of radiative transfer

$\frac{dF}{dp_{\uparrow}}=\frac{dF}{dp_{\downarrow}}=\frac{\kappa}{2g}(F_{\uparrow}-F_{\downarrow})$
(6)

with appropriate boundary conditions, we can analytically solve for the upward
and downward heat flux at every level in the stable layer. Of interest for our
problem is the net heat flow from the deep/stable layers to the shallow layer
that is convectively coupled to the photosphere. This is

$F_{\uparrow}(p_{1})-F_{\downarrow}(p_{1})=\frac{4F_{2}gp_{0}+F_{1}(p_{2}^{2}-p_{1}^{2})\kappa_{0}}{4gp_{0}+(p_{2}^{2}-p_{1}^{2})\kappa_{0}}-F_{1}$
(7)

where $F_{1}=F_{\downarrow}(p_{1})$ and $F_{2}=F_{\uparrow}(p_{2})$ are the
boundary conditions. We can solve for these boundary conditions using the
uniquely determined temperature structures above and below the stable layer,
i.e.

$F_{1}=\frac{\kappa_{0}\sigma}{g}\int_{0}^{p_{1}}\exp\left[\frac{\kappa_{0}(p^{2}-p_{1}^{2})}{2gp_{0}}\right]T(p)^{4}(p/p_{0})dp$
(8) $F_{2}=\frac{\kappa_{0}\sigma
T_{2}^{4}}{g}\int_{p_{2}}^{\infty}\exp\left[\frac{\kappa_{0}(p_{2}^{2}-p^{2})}{2gp_{0}}\right]\left(p/p_{2}\right)^{4R/c_{p}}(p/p_{0})dp$
(9)

Using this process, for a given
$(T_{\text{int}},q_{\text{max}})\implies(p_{2},T_{2})$ we solve the above
nonlinear equation to obtain $T_{e}\implies(p_{1},T_{1})$ using the condition
$F_{\uparrow}(p_{1})-F_{\downarrow}(p_{1})=\sigma(T_{e}^{4}-fT_{0}^{4})$. This
defines the observed effective temperature as a function of the internal
effective temperature $T_{e}(T_{\text{int}})$, shown as the dashed-dotted
curves in Figure 4(a). In order to relate the observed effective temperature
$T_{e}$ to the adiabatic effective temperature $T_{ab}$, one must additionally
consider the effect of latent heat, shown as solid curves in Figure 4(a). Then
the net effect, $T_{e}(T_{ab})$ is shown as dashed curves in Figure 3, where
$\Delta T=T_{e}-T_{ab}$. Figure 1 shows a sample temperature/pressure profile
to illustrate the contributions from each source. Similar results are shown
for water in Figure 4. These figures show $\Delta T$ under two assumptions: an
initially dry adiabatic atmospheric profile (solid curves) where moist
adiabaticity is not considered, and an initially moist adiabatic upper layer
(dashed curves), as described in detail in the text. In the dry adiabatic
case, $T_{\text{int}}=T_{ab}$. One may notice that $\Delta T=T_{e}-T_{ab}$
arising from the dry adiabatic case is substantially smaller than
$T_{e}-T_{\text{int}}$ in the wet adiabatic case (dashed-dotted curves in
Figure 4(a)). There are two reasons for this. First, the dry adiabatic lapse
rate is steeper than the wet pseudo-adiabatic lapse rate, which improves the
efficiency of radiative transfer. Second, the super-adiabatic stable region
makes less of a difference when the initial profile was already comparatively
steeper than the moist adiabatic case. For this reason, the results are less
sensitive to assumptions about the initial temperature profile than one might
initially expect.

(a) Effective temperature changes arising from different sources. Latent heat
(solid), convective inhibition (dashed-dotted), and the net effect (dashed).

(b) Net temperature differences, for an initially dry adiabatic atmosphere
(solid) and an initially moist adiabatic upper layer (dashed)

Figure 3: $\Delta T(T_{ab})$ for different envelope abundances of methane
$q_{\text{max}}$ between 2-5%, where $\Delta T\equiv T_{e}-T_{ab}$. The dashed
(net $\Delta T$) curves are the same in (a) and (b).

(a) Effective temperature changes arising from different sources. Latent heat
(solid), convective inhibition (dashed-dotted), and the net effect (dashed).

(b) Net temperature differences, for an initially dry adiabatic atmosphere
(solid) and an initially moist adiabatic upper layer (dashed)

Figure 4: $\Delta T(T_{ab})$ for different envelope abundances of water
$q_{\text{max}}$ between 3-12%, where $\Delta T\equiv T_{e}-T_{ab}$. The
dashed (net $\Delta T$) curves are the same in (a) and (b).

Finally we must consider the deposition of sunlight. The primary absorber of
sunlight in Uranus’ troposphere is methane vapor [Marley and McKay, 1999],
while the primary absorber of infrared light is hydrogen collisions. At 1 bar
with 2% methane for example, accounting for methane absorption plus Rayleigh
scattering vs. thermal absorption by hydrogen, the ratio between visible to
thermal infrared opacity is approximately $\kappa_{\odot}/\kappa_{t}<10^{-2}$
averaging over a broad band in wavelength. This value of course is not unique;
there are windows at certain wavelengths, the methane abundance changes
rapidly with depth, and this simple calculation neglects absorption by haze
and cloud particles. Nevertheless $\kappa_{\odot}/\kappa_{t}$ is sufficiently
small that we can plausibly argue that sunlight penetrates significantly
beyond the 1 bar level, and most of the sunlight is absorbed deeper in the
atmosphere. This can be parameterized by simply arguing some fraction $f$ of
sunlight is absorbed above the cloud level, and $1-f$ is absorbed below. In
principle this procedure accommodates any value of $f$, but for our purposes
for simplicity we approximate using the limiting cases $f\rightarrow 0$ for
the shallow ($\sim 1$bar) methane condensation level, while $f\rightarrow 1$
for deep ($\sim 100$bar) water clouds. The difference for the stable layer is
that in radiative equilibrium, the flux through the stable layer must balance
with $\sigma(T_{e}^{4}-fT_{0}^{4})$ as explained above.

###  3.3 Meteorological implications

The equilibrium solution will be stable if the virtual potential temperature
$\theta_{v}$ is monotonically decreasing with increasing pressure between
$p_{1}$ and $p_{2}$.

$\theta_{v}(p,T)=T_{v}(p,T)\left(\frac{p_{0}}{p}\right)^{\nabla_{ab}}$ (10)
$T_{v}(p,T)=T\left(1-q(1-\epsilon)\right)^{-1}$ (11)

In the equilibrium cases discussed here, this condition is always satisfied
(see dashed red curve in Figure 1). This result should be fully general; it
does not depend on, for example, the choice of atmospheric opacity. In order
to trigger a convective instability, the upper layer would need to become more
dense than the deep layer (note in this section we use the word “dense” in the
virtual potential temperature sense, i.e. accounting implicitly for adiabatic
expansion/compression). In order to accomplish this by cooling, the upper
layer would need to cool such that the bottom of the upper layer satisfies
$q<q_{\text{crit}}$, so that further cooling makes the gas mixture more dense
rather than less dense. However, as soon as $q$ becomes infinitesimally less
than $q_{\text{crit}}$, it will be more dense than the material in the stable
layer directly beneath it, even while remaining less dense than the well mixed
gas in the deep layer. This will cause a small convective instability wherein
a portion of the stable layer is eroded into the upper layer, causing the
stable layer to thin and restoring the upper layer to satisfy
$q=q_{\text{crit}}$ at a new pressure level. This will happen in all cases
where the stable layer is of finite thickness. Therefore, in order to trigger
a convective instability with the deep layer, the stable layer must vanish
completely, reducing to a compositional discontinuity between the upper and
deep layers. Such a scenario in the limit of heat transport by thermal
conduction and an inviscid fluid results in an infinitesimal thermal boundary
layer with an infinitesimal temperature discontinuity. In the optically thin
radiative transfer case, the situation is somewhat more subtle, because rather
than heat flux diverging to infinity, it converges on a finite value (in the
two stream approximation: the upward heat flux from below minus the downward
heat flux from above). Nevertheless, if there is a significant temperature
difference between the upper and deep layer (as would be required in order to
trigger a convective instability) this finite value is orders of magnitude
larger than the luminosity of the planets in question, for any appreciable
temperature discontinuity. This result contrasts with previous findings [Li
and Ingersoll, 2015], which posited the system may behave as a relaxation
oscillator when the upper layer cools sufficiently to become over-dense and
trigger a convective instability. This previous study did not explicitly
account for radiative transfer across the stable layer, instead dynamically
cooling from above and treating the stable layer as a perfect insulator. In
our case, we find instead that the profile relaxes into a state of global
radiative-convective equilibrium; eventually the upper layer stops cooling as
the heat it loses to space balances with the heat radiated across the stable
layer.

After solving for the equilibrium atmospheric profile, it is possible to
compute the convective available potential energy (CAPE), as well as the
activation energy required to disrupt the equilibrium stable layer. Both
quantities can be computed in the same way, using

$E=\int_{p_{a}}^{p_{b}}R_{d}(T_{v,f}-T_{v,i})\frac{dp}{p}$ (12)

where $T_{v,f}$ is the final virtual temperature profile (for example, a moist
adiabat) , $T_{v,i}$ is the initial virtual temperature profile (e.g. in
radiative-convective equilibrium with a stable layer), and $p_{a}$ and $p_{b}$
are bounding pressures that satisfy $T_{v,i}=T_{v,f}$. The blue and red dashed
curves in Figure 1 provide a visual example for how this calculation can be
done. This computed activation energy in units of energy per mass, can be
converted into updraft velocities required to disrupt the stable layer. In the
case of Saturn, if $q_{\text{crit}}\sim q_{\text{max}}$ as suggested in [Li
and Ingersoll, 2015], then the activation energy required to disrupt this
stable equilibrium is quite small, requiring updraft velocities of only a few
meters per second. Therefore the basic premise of [Li and Ingersoll, 2015] can
still be valid, although it may require some additional mechanism to jump-
start the process, for example a strong updraft or entrainment erosion as
described in [Friedson and Gonzales, 2017]. In the case of Uranus and Neptune
with methane mixing ratios further from the critical value, the activation
energy necessary to initiate a convective instability is correspondingly
larger, requiring updraft velocities of order tens of meters per second, far
larger than expected convective velocities. However, if such a disruption were
able to occur by some anomalous updraft, the resulting storm would be
extremely energetic, with CAPE exceeding 10 J/g, substantially larger than
even the most extreme weather events on Earth.

###  3.4 Interior implications

Figure 5: Apparent effective temperature $T_{e}$ as a function of internal
equivalent effective temperature $T_{\text{int}}$ considering both methane and
water. Colors correspond to water abundance, increasing downward from 0% to
12% water. The curve style corresponds to methane abundance, increasing
downward from 0% to 5% abundance.

Interior models of giant planets generally assume a 1-bar equivalent
temperature that sets the internal entropy of the planet. This is usually done
using the atmospheric temperature, and corrected for any non-adiabatic
behavior in the atmosphere. Non-adiabatic models for Uranus’ and Neptune’s
interiors and thermal evolution have recently been carried out (e.g.
[Nettelmann et al., 2016], [Vazan and Helled, 2020], [Scheibe et al., 2021]),
finding self-consistent solutions to Uranus’ contemporary heat flow. Our
results are still relevant to many of these models. The so called thermal
boundary layers at depth from [Nettelmann et al., 2016] and [Vazan and Helled,
2020] retain an adiabatic convective envelope set by a 1-bar equivalent
temperature. Furthermore, the U-1 and U-2 models from [Vazan and Helled, 2020]
involve a convective envelope of homogeneous composition, whose temperature
profile is also set using a 1-bar equivalent temperature. Although our results
make no direct statement about the behavior of the deep interior of the planet
that may include extended regions of static stability, it still applies to the
thermal evolution of these upper envelopes. Because Uranus and Neptune possess
magnetic fields, and the latter possesses a substantial internal heat flow, we
deduce that both must be convective at depth, and that this adiabatic
description is probably relevant for at least some fraction of the interior of
Uranus and Neptune. The primary source of non-adiabatic behavior in the
atmosphere is condensation, usually accounted for by the sub-adiabatic
gradients caused by latent heat of condensation. In Figure 5, we show the
apparent effective temperature $T_{e}$ against the internal equivalent
effective temperature $T_{\text{int}}$ for various envelope condensate
abundances. If convective inhibition does occur, it will largely cancel out
the effect of moist adiabaticity. Therefore, using a dry adiabat to guess the
1-bar equivalent temperature is a better approximation than accounting for
latent heat alone but neglecting convective inhibition. A better approach
could be to use the analytic scaling relationships in Section 2 to estimate
the magnitude of these effects. The best approach would be to explicitly model
the effect of convective inhibition, using methods from this work or [Leconte
et al., 2017].

##  4 Evolutionary model

We present an adiabatic thermal evolution model of Uranus and Neptune. As
discussed, treating the interior as adiabatic is probably inappropriate for
the ice giants [Helled et al., 2020]. Nevertheless it provides a convenient
framework to understand the effect of methane condensation on these planets’
thermal histories in the absence of an accepted interior model. For an
adiabatic model, we assume the total heat content of the planet’s interior to
be a linear function of its adiabatic equivalent effective temperature
$T_{ab}$.

$\int_{M_{\text{min}}}^{M}c_{p}Tdm=A\bar{c_{p}}MT_{ab}$ (13)

One way of thinking about this equation is to imagine a small set of layers,
or possibly even one layer, in the form of concentric shells, each of which is
isentropic and homogeneous but of different composition to neighboring layers,
with negligible thermal boundary layers between them as would be fluid
dynamically expected for a low viscosity system. Beneath this set of shells
there could be a region, possibly a substantial fraction of the planet, where
there is a compositional gradient and therefore inefficient convective
transport. This deeper region would not contribute to $A$ or to the resulting
thermal evolution of the planet because it stores primordial heat. In reality,
there would be non-zero thermal diffusion from a stably stratified interior
portion of the planet and its convective envelope, slowing down planetary
cooling. However, if the diffusion timescale for the planet is longer than the
age of the solar system, then this contribution would be small. It is not
guaranteed that the diffusion timescale is in fact longer than the age of the
solar system. It depends on the (unknown) thermal transport properties of the
(unknown) compositional constituents of the ice giants’ interiors, and is
further complicated by the possibility of thermal transport by double
diffusive convection. Therefore we acknowledge this description of the
interior evolution is imperfect, but it does at least approximately describe a
wide variety of possible interior behaviors, and provides a convenient
framework to self-consistently assess the relative influence of convective
inhibition on thermal evolution while remaining agnostic about the details of
the ice giants’ interior structures. A full description of Uranus’ and
Neptune’s thermal evolutions would require a detailed interior model.
Nevertheless under our assumptions, the parameter $A$ is approximately
constant through time because the Gruneisen parameter is rather insensitive to
temperature, and its value is set by the fraction of the total mass that is
fully convective. Some fraction of the planets must be convective in order to
generate their observed magnetic fields.

The rate of cooling depends on the apparent effective temperature $T_{e}$,
while the heat content of the bulk of the interior is linearly related to the
adiabatic equivalent effetive temperature $T_{ab}$. The equation governing
thermal evolution is then

$4\pi R^{2}\sigma(T_{e}^{4}-T_{0}^{4})=-A\bar{c_{p}}M\frac{dT_{ab}}{dt}$ (14)

where $T_{0}$ is the equilibrium effective temperature of the planet with
sunlight if there were no internal heat source. The radius $R$ is treated as
constant because we are concerned with most of the evolution where the body is
degenerate, not any early very hot phase. The steady increase in solar
luminosity (i.e. time variation of $T_{0}$) is ignored. In order to solve the
thermal evolution equation, we need an explicit relationship between $T_{e}$
and $T_{ab}$. This is done using the method from Section 3, with the results
for methane in Figure 3 and for water in Figure 4. In general the relationship
between $T_{e}$ and $T_{ab}$ depends on the condensate bulk interior abundance
$q_{\text{max}}$. Without condensation, the relationship is $T_{e}=T_{ab}$. In
this case, Equation 14 can be non-dimensionalized into the canonical form

$\frac{x^{4}-x_{0}^{4}}{1-x_{0}^{4}}=-\tau_{K}\frac{dx}{dt}$ (15)

where $x\equiv T_{e}/T_{e}^{(\emptyset)}$ where $T_{e}^{(\emptyset)}$ is the
apparent effective temperature today, and $x_{0}=T_{0}/T_{e}^{(\emptyset)}$.
The Kelvin timescale $\tau_{K}=\frac{A\bar{c_{p}}MT_{e}^{(\emptyset)}}{4\pi
R^{2}\sigma}(T_{e}^{(\emptyset)4}-T_{0}^{4})^{-1}$ scales how long it takes to
cool to $T_{e}^{(\emptyset)}$ from an initial arbitrarily hot state. In the
asymptotic case $T_{e}^{(\emptyset)}\gg T_{0}$ (not true for Uranus!), this is
about $\tau_{K}/4$, but can be a different fraction of $\tau_{K}$ in general.

Accounting for condensation, Equation 14 becomes

$\frac{dx}{dx_{ab}}\frac{x^{4}-x_{0}^{4}}{1-x_{0}^{4}}=-\tau_{K}\frac{dx}{dt}$
(16)

where $x_{ab}=T_{ab}/T_{e}^{(\emptyset)}$. The difference between Equations 15
and 16 then straightforwardly demonstrates the effect of convective inhibition
by condensation on the planets’ thermal evolution: it alters the rate of
cooling by a factor of $\frac{dx}{dx_{ab}}$. As we will see, this factor can
be greater than or less than unity. This means the effect can either speed up
or slow down the rate of change of the planets’ apparent effective temperature
at different points in its thermal history. This is especially important for
understanding the results for water.

The fact that Equation 16 retains the Kelvin timescale $\tau_{K}$ makes this
formulation especially convenient. This allows us to directly compare the
fraction of that timescale that a given evolutionary model takes to cool from
arbitrarily hot bodies to their current temperatures for different assumptions
of the condensate abundance $q_{\text{max}}$. Leaving the the effect in terms
of the Kelvin timescale allows our results to be roughly independent of
accurate interior models, because $\tau_{K}$ implicitly encodes an arbitrary
interior model. The results for methane are shown in Figure 6, and for water
in Figure 7. The results for methane are relatively straightforward; for the
early stages of Uranus and Neptune’s thermal histories, the effect of methane
is unimportant, because the atmosphere is warm enough that methane does not
condense anywhere. As the atmosphere cools, methane begins to condense, at
first in the stratosphere above the radiative convective boundary. As cooling
continues, convective inhibition begins to extend the radiative-convective
boundary downward, as superadiabatic gradients can be stable. At this point,
the measured effective temperature $T_{e}$ departs from its adiabatic
equivalent $T_{ab}$, causing the effective temperature to drop faster than the
interior is cooling. As cooling continues, the layered system described in
Section 3 emerges, and perhaps persists today [Guillot, 1995].

(a) Uranus thermal evolution

(b) Neptune thermal evolution

Figure 6: Thermal evolution model for Uranus and Neptune, with different
colored curves representing different methane abundances. The x-axis is the
time before the present day, scaled to the cooldown time in the dry adiabatic
(no condensation) case. The y-axis is $x$ from Equation 16. Line styles and
colors are identical to Figure 3.

The case of water, shown in Figure 7, the behavior is more subtle. In this
case, thermal evolution is actually slowed down compared to the adiabatic
case. This contrasts with previous findings [Kurosaki and Ikoma, 2017], which
considered the effect of moist adiabaticity (i.e., latent heat) but did not
quantify the effects of convective inhibition. If we consider latent heat
only, we obtain results in good agreement with this previous study. In their
case, thermal evolution is sped up, because the atmosphere initially remains
warm while the interior cools. This allows the planet to lose heat efficiently
when condensation first occurs, speeding up evolution. Our findings
demonstrate that the effect of convective inhibition overwhelms the effect of
moist adiabaticity, so that our story is the opposite. Early on, as
condensation occurs, we find the atmosphere cools faster than the interior,
reducing cooling efficiency. Later on, $\frac{dx}{dx_{ab}}$ from Equation 16
becomes less than unity, as demonstrated by the negative slope at low
temperatures in Figure 4. Therefore in the case of water, condensation early
in the ice giants’ thermal histories caused the effective temperature to drop
faster than the internal temperature, analogous to what happened with methane
condensation more recently. However, this temporary speedup of
$\frac{dT_{e}}{dt}$ coincides with a loss of luminosity, slowing down the rate
at which the interior loses heat. Then, over subsequent evolution, the
interior cools inefficiently, and in recent history the effective temperature
changes slower than the internal temperature. The net effect is a cooldown
time that is longer than the dry adiabatic case.

(a) Uranus thermal evolution

(b) Neptune thermal evolution

Figure 7: Thermal evolution model for Uranus and Neptune, with different
colored curves representing different water abundances. Axes are identical to
Figure 6. Line styles and colors are identical to Figure 4.

Because the water and methane cloud decks are well separated, the
superposition of the two effects is straightforward. Immediately beneath the
methane cloud deck, the behavior can be accurately modeled as a dry adiabat,
because the water mixing ratio is so small at these relatively low
temperatures. Therefore the results of modeling the whole atmosphere with both
cloud decks explicitly is virtually identical to using $\Delta T_{tot}=\Delta
T_{CH_{4}}+\Delta T_{H_{2}O}$ from Figures 3 and 4.

##  5 Discussion

Provided $q_{\text{max}}>q_{\text{crit}}$, when the planet cools to a
temperature low enough for condensation, convection can be interrupted. At
this point the apparent effective temperature departs from the internal
equivalent effective temperature, by the mechanism described in Section 3. We
can solve for the equilibrium configuration to derive the apparent effective
temperature $T_{e}$ as a function of $q_{\text{max}}$ and $T_{\text{int}}$.
These results are shown as the dashed-dotted lines for methane in Figure 3,
and Figure 4 for water. Interior modelers should bear in mind that the
internal 1-bar equivalent temperature may depart from simple adiabatic
extrapolation of the troposphere by nearly a significant factor (see Section
3.4 and Figure 5).

This behavior leads to the evolutionary behavior observed in Figure 6. Before
condensation occurs, the planet cools normally. Upon the onset cloud
formation, the apparent effective temperature drops rapidly. However, upon
reaching the minimum, the apparent effective temperature actually begins to
decrease more slowly than the fully adiabatic case. This effect is present in
Figure 3 but is more apparent in Figure 4. The net effect for the ice giants
is a net speedup of thermal evolution for methane, and a net slowdown for
water. The magnitude of this speedup or slow down can be no more than 15% in
either case, assuming 5% methane molar abundance or 12% water molar abundance.
Both effects can occur simultaneously, and superimpose straightforwardly
because their cloud decks are well separated.

As Neptune continues to cool, methane will begin to behave similarly to water,
exhibiting a local minimum in $\Delta T(T_{ab})$. This local minimum is seen
most clearly for water in Figure 4 around $T_{ab}=$180K, but can also be seen
for the dashed curves in Figure 3 around 75K, and would be present in the
solid curves if the x-axis extended to lower temperatures. As the planet
continues to cool below this local minimum, the slope of $\Delta T(T_{ab})$
becomes negative, and the rate of change of the thermal state of the
atmosphere slows. Consider the implications of this for methane clouds near
the 1-bar level. This state persists for longer from a thermal evolution
perspective than an arbitrary/random thermal state. That is, this state is a
local minimum in $\Delta T(T_{ab})$, meaning the planet reaches this state
faster than it would if it were cooling adiabatically, and leaves this state
more slowly than it would if it were cooling adiabatically. Therefore these
planets will spend a longer portion of their thermal histories in the state
where the cloud level is $\sim 1-10$bars than they would in a thermal
evolution model that does not consider convective inhibition by condensation.
Perhaps this consideration renders the surprising similarity of Uranus and
Neptunes’ atmospheres’ shallow temperature structures despite their vast
difference in insolation somewhat less improbable than it first appears.

We must consider whether this atmospheric structure is compatible with
existing data, especially Voyager radio refractivity data. The current data
has been shown to be consistent with many different models, including
subadiabatic, adiabatic, moist adiabatic, and superadiabatic temperature
gradients [Helled et al., 2020]. The data has also been shown to be compatible
with a wide range of temperature structures and methane abundances [Lindal et
al., 1987] [Lindal, 1992]. The data itself shows a layer of rapidly varying
refractivity near the condensation level, generally interpreted to be methane
clouds [Lindal et al., 1987] [Lindal, 1992][Marley and McKay, 1999]. Another
interpretation of the same data supports a layer of superadiabatic temperature
lapse rate in the cloud-forming regions of these planets [Guillot, 1995]. In
general, our understanding of the thermal structure of the ice giant
atmospheres is incomplete, as the results from Voyager 2 refractivity data are
model dependent, with a particular degeneracy between assumed methane
enrichment and temperature structure. In order to disentangle these variables
and have a more confident understanding of these planets’ atmospheres’ thermal
structures, we must return with a mission. It should be a priority for a
future mission to independently measure methane abundance and temperature,
perhaps with entry probes or a well designed microwave radiometer experiment.

These general findings do not consider the long term stability of stable
layers in the atmosphere. As long as the stability timescale is greater than
the relaxation timescale for a stable layer, the results should approximately
reflect reality. However, the stability timescale is poorly constrained
[Friedson and Gonzales, 2017]. If it is sufficiently short, this could further
complicate the dynamics. If that condition is satisfied, then even if stable
layers are intermittently interrupted by massive internal plumes, large meteor
impacts, or instability due to long term erosion by entrainment, they will
reform again on geologically short timescales ($\sim 100$yr). Therefore the
thermal evolution will be governed primarily by the equilibrium state, and not
possible intermittent periods of enhanced activity. Intermediate states where
the equlibrium configuration is thinned over time but not totally destroyed by
entrainment erosion would in general reduce the magnitude of $\Delta
T(T_{\text{int}})$, so the findings in this paper should be considered an
upper bound. Furthermore, we use a highly simplified thermal evolution model,
not considering changes in planetary radius or explicitly accounting for the
effects of non-adiabaticity at depth. Seasonal variations in insolation were
not included in the model, as these variations average out over geologic time.
However, seasonal variations have been shown to create local temperature
variations of order 10K [Orton et al., 2007], comparable to the magnitude of
the effect of convective inhibition by methane. The possible dynamical and
evolutionary consequences could be the subject of future work. Our atmospheric
model also did not explicitly include the condensate opacities, and may
therefore not capture possible feedback mechanisms. We discuss further the
possible effects of opacity variation due to condensation in the following
paragraphs. For these reasons, this work should be considered exploratory, and
further work is needed in order to more confidently establish the thermal
histories of the ice giants while accounting for convective inhibition.

Here we must include a discussion about the effects of opacity variation due
to condensation, which are not considered in this model but which are
certainly important for a fully complete understanding of Uranus and Neptune’s
thermal states and thermal histories, and has been considered explicitly by
prior works, e.g. [Kurosaki and Ikoma, 2017]. The variation of opacity affects
our results in two important ways: first, by changing the radiative-convective
boundary as vapor condenses out of the stratosphere; and second, by affecting
radiative transfer within and across the layer of stable stratification caused
by convective inhibition.

We begin by discussing the stratospheric effect of opacity variations due to
condensation. Water and methane are both more opaque than hydrogen in the
thermal infrared, therefore as the planet cools and these volatiles begin to
condense and rain out of the stratosphere, the stratosphere becomes more
transparent and the radiative-convective boundary deepens. At fixed effective
temperature, the temperature at the radiative-convective boundary is
relatively unchanged, therefore decreasing the opacity of the stratosphere has
the net effect of decreasing the entropy of the troposphere at fixed effective
temperature. Therefore during this stage, the temperature of the troposphere
is cooling faster than the effective temperature of the planet as the
stratosphere extends downward. Convective inhibition only begins to become
relevant when the stratosphere has cooled sufficiently such that the
radiative-convective boundary has a lower vapor mixing ratio than the bulk
abundance. By the time this occurs, the bulk of the stratosphere is cooler
than the radiative convective boundary by approximately a factor of $2^{1/4}$,
and therefore relatively dry due to the highly sensitive dependence of
saturation vapor pressure on temperature. So the important stratospheric
effect due to opacity variation we have just described qualitatively,
essentially predates the onset of convective inhibition. This allows us to
neglect these dynamics in our context, although we caution the reader that a
fully realistic consideration of the effects of condensation must also include
the effects of opacity variations, which are important.

The second effect of opacity variation is on the radiative transfer across the
stable layer. To estimate the importance of this effect, we modifed our method
so that we increased the opacity of the deep layer and and stable layer (see
Figure 2) by an order of magnitude. This changes our results for contemporary
methane clouds by no more than 2%, and for contemporary deep water clouds by
$<$0.01%. The effect for water clouds is larger earlier in its evolution, but
is always $<$2%. In either case, the direction of this consideration is to
increase the magnitude of $\Delta T$. The reason this matters less for deep
clouds is because the opacity is already very large, and the effect due to
deep clouds is accurately approximated by Equation 2, which assumes a high
opacity limit. For shallower clouds where the details of thermal transport are
more relevant, it affects the results, but only as a relatively small
correction even assuming a very large order of magnitude change in opacity.

If there are indeed layers of static stability in the troposphere or deep
atmosphere of Uranus and/or Neptune, then they should support gravity waves.
Whether we expect gravity waves to be excited, what their general
characteristic would be, and whether they could be detected from space (for
example using an Doppler imager) is a subject worthy of future theoretical
consideration.

Whatever the uncertainties about the specifics, the basic physical mechanism
is likely to be important in the ice giants because of their highly enriched
atmospheres. There may be additional stable layers, for example a silicate
cloud level beneath the water cloud level, or a sulfide/ammonia cloud level.
We focus on only two in this work to demonstrate the general principle without
getting bogged down in largely unconstrained assumptions about the envelope
enrichment in each species. However, the intuition we build here for methane
and water can be straightforwardly applied to other cloud levels using exactly
the same method. This method is also likely to be applicable to the majority
of exoplanets, ranging from super-Earths and water worlds with hydrogen
envelopes, to metal-enriched gas giants. It is clear from this work that
thermal evolution and internal thermal structure may be profoundly influenced
by convective inhibition by condensation. Any complete model of thermal
evolution or internal structure is advised to consider convective inhibition.

#### Acknowledgments

This work has been funded by the NASA FINESST program. This work has greatly
benefited from the ISSI Ice Giants Science Team meetings of 2019 and 2020 in
Bern. We would also like to thank the anonymous reviewers for their invaluable
feedback that greatly strengthened this work.

## References

  * Baines and Hayden Smith, 1990 Baines, K. and Hayden Smith, W. (1990).  Icarus, 85:65–108.
  * Baines et al., 1995 Baines, K., Mickelson, M., Larson, L., and Ferguson, D. (1995).  Icarus, 114:328–340.
  * Fazio et al., 1976 Fazio, G., Traub, W., Wright, E., Low, F., and Trafton, L. (1976).  ApJ, 209:633–657.
  * Fortney et al., 2011 Fortney, J., Ikoma, M., Nettelmann, N., Guillot, T., and Marley, M. (2011).  ApJ, 729(32).
  * Friedson and Gonzales, 2017 Friedson, A. and Gonzales, E. (2017).  Icarus, 297:160–178.
  * Guillot, 1995 Guillot, T. (1995).  Science, 269:1697–1699.
  * Guillot, 2005 Guillot, T. (2005).  AREPS, 33:493–530.
  * Gulkis et al., 1983 Gulkis, S., Olsen, E., Klein, M., and Thompson, T. (1983).  Science, 221(4609):453–455.
  * Helled et al., 2020 Helled, R., Nettelmann, N., and Guillot, T. (2020).  SSRv, 216(38).
  * Hubbard, 1977 Hubbard, W. (1977).  Icarus, 30(2):305–310.
  * Hubbard, 1978 Hubbard, W. (1978).  Icarus, 35:177–181.
  * Hubbard et al., 1999 Hubbard, W., Guillot, T., Marley, M., Burrows, A., Lunine, J., and Saumon, D. (1999).  P&SS, 47:1175–1182.
  * Hubbard and MacFarlane, 1980 Hubbard, W. and MacFarlane, J. (1980).  JGR, 85(B1).
  * Karkoschka and Tomasko, 2011 Karkoschka, E. and Tomasko, M. (2011).  Icarus, 211:780–797.
  * Kellermann and Pauliny-Toth, 1966 Kellermann, K. and Pauliny-Toth, I. (1966).  Observations of the radio emission of uranus, neptune, and other planets at 1.9cm.  ApJ, 145(3):954–956.
  * Kurosaki and Ikoma, 2017 Kurosaki, K. and Ikoma, M. (2017).  ApJ, 153(260):9.
  * Leconte et al., 2017 Leconte, J., Selsis, F., Hersant, F., and Guillot, T. (2017).  A&A, 598.
  * Li and Ingersoll, 2015 Li, C. and Ingersoll, A. (2015).  NatGe, 8:398–403.
  * Lindal, 1992 Lindal, G. (1992).  ApJ, 103(3).
  * Lindal et al., 1987 Lindal, G., Lyon, J., Sweetnam, D., Eshleman, V., et al. (1987).  JGRA, 92(A13):14987–15001.
  * Linder et al., 2019 Linder, E., Mordasimi, C., Molliere, P., Dominik-Marleau, G., Malik, M., et al. (2019).  A&A, 623(A85).
  * Loewenstein et al., 1977 Loewenstein, R., Harper, D., and Moseley, H. (1977).  ApJ, 218:L145–L146.
  * Marley and McKay, 1999 Marley, M. and McKay, C. (1999).  Icarus, 138(2):268–286.
  * Nettelmann et al., 2016 Nettelmann, N., Wang, K., Fortney, J., Hamel, S., Yellamilli, S., et al. (2016).  Icarus, 275:107–116.
  * Orton et al., 2007 Orton, G., Encrenaz, T., Leyrat, C., Puetter, R., and Friedson, A. (2007).  A&A, 473(L5-L8).
  * Pearl and Conrath, 1991 Pearl, J. and Conrath, B. (1991).  JGR, 96:18921–18930.
  * Pearl et al., 1990 Pearl, J., Conrath, B., Hanel, R., Pirraglia, J., and Coustenis, A. (1990).  Icarus, 84:12–28.
  * Podolak et al., 1991 Podolak, M., Hubbard, W., and Stevenson, D. (1991).  Models of uranus’ interior and magnetic field.  In Uranus, pages 29–61. University of Arizona Press, Tucson, AZ.
  * Rages et al., 1991 Rages, K., Pollack, J., Tomasko, M., and Doose, L. (1991).  Icarus, 89:359–376.
  * Scheibe et al., 2019 Scheibe, L., Nettelmann, N., and Redmer, R. (2019).  A&A, 632(A70).
  * Scheibe et al., 2021 Scheibe, L., Nettelmann, N., and Redmer, R. (2021).  A&A, in review.
  * Stevenson, 1983 Stevenson, D. (1983).  JGR, 88(B3):2445–2455.
  * Stier et al., 1978 Stier, M., Traub, W., Fazio, G., Wright, E., and Low, F. (1978).  ApJ, 226:347–349.
  * Vazan and Helled, 2020 Vazan, A. and Helled, R. (2020).  A&A, 633(A50).


IDS Collaboration

# On the nature of yrast states in neutron-rich polonium isotopes

R. Lică<EMAIL_ADDRESS>Horia Hulubei National Institute for R&D in
Physics and Nuclear Engineering, RO-077125 Bucharest, Romania  A.N. Andreyev
School of Physics, Engineering and Technology, University of York, York, YO10
5DD, United Kingdom  H. Naïdja  Université Constantine 1, Laboratoire de
Physique Mathématique et Physique Subatomique, Constantine 25000, Algeria  A.
Blazhev  Institut für Kernphysik, Universität zu Köln, D-50937 Köln, Germany
P. Van Duppen  KU Leuven, Instituut voor Kern- en Stralingsfysica, B-3001
Leuven, Belgium  B. Andel  Department of Nuclear Physics and Biophysics,
Comenius University in Bratislava, 84248 Bratislava, Slovakia  A. Algora
Instituto de Física Corpuscular, CSIC - Universidad de Valencia, E-46980,
Valencia, Spain  Institute of Nuclear Research (ATOMKI), P.O.Box 51, H-4001
Debrecen, Hungary  S. Antalic  Department of Nuclear Physics and Biophysics,
Comenius University in Bratislava, 84248 Bratislava, Slovakia  J. Benito
Grupo de Fìsica Nuclear, EMFTEL & IPARCOS, Universidad Complutense de Madrid,
28040 Madrid, SPAIN  G. Benzoni  Istituto Nazionale di Fisica Nucleare,
Sezione di Milano, I-20133 Milano, Italy  T. Berry  Department of Physics,
University of Surrey, Guildford GU2 7XH, United Kingdom  M. J. G. Borge
Instituto de Estructura de la Materia, CSIC, Serrano 113 bis, E-28006 Madrid,
Spain  C. Costache  Horia Hulubei National Institute for R&D in Physics and
Nuclear Engineering, RO-077125 Bucharest, Romania  J. G. Cubiss  School of
Physics, Engineering and Technology, University of York, York, YO10 5DD,
United Kingdom  H. De Witte  KU Leuven, Instituut voor Kern- en
Stralingsfysica, B-3001 Leuven, Belgium  L. M. Fraile  Grupo de Fìsica
Nuclear, EMFTEL & IPARCOS, Universidad Complutense de Madrid, 28040 Madrid,
SPAIN  H. O. U. Fynbo  Department of Physics and Astronomy, Aarhus University,
DK-8000 Aarhus C, Denmark  P. T. Greenlees  University of Jyväskylä,
Department of Physics, Accelerator laboratory, P.O. Box 35(YFL) FI-40014
University of Jyväskylä, Finland  L. J. Harkness-Brennan  Department of
Physics, Oliver Lodge Laboratory, University of Liverpool, Liverpool L69 7ZE,
United Kingdom  M. Huyse  KU Leuven, Instituut voor Kern- en Stralingsfysica,
B-3001 Leuven, Belgium  A. Illana  Instituto Nazionale di Fisica Nucleare,
Laboratori Nazionali di Legnaro, I-35020 Legnaro, Italy  J. Jolie  Institut
für Kernphysik, Universität zu Köln, D-50937 Köln, Germany  D. S. Judson
Department of Physics, Oliver Lodge Laboratory, University of Liverpool,
Liverpool L69 7ZE, United Kingdom  J. Konki  University of Jyväskylä,
Department of Physics, Accelerator laboratory, P.O. Box 35(YFL) FI-40014
University of Jyväskylä, Finland  I. Lazarus  STFC Daresbury, Daresbury,
Warrington WA4 4AD, United Kingdom  M. Madurga  Dept. of Physics and
Astronomy, University of Tennessee, Knoxville, Tennessee 37996, US  N.
Marginean  Horia Hulubei National Institute for R&D in Physics and Nuclear
Engineering, RO-077125 Bucharest, Romania  R. Marginean  Horia Hulubei
National Institute for R&D in Physics and Nuclear Engineering, RO-077125
Bucharest, Romania  C. Mihai  Horia Hulubei National Institute for R&D in
Physics and Nuclear Engineering, RO-077125 Bucharest, Romania  R. E. Mihai
Horia Hulubei National Institute for R&D in Physics and Nuclear Engineering,
RO-077125 Bucharest, Romania  P. Mosat  Department of Nuclear Physics and
Biophysics, Comenius University in Bratislava, 84248 Bratislava, Slovakia  J.
R. Murias  Grupo de Fìsica Nuclear, EMFTEL & IPARCOS, Universidad Complutense
de Madrid, 28040 Madrid, SPAIN  Institut Laue-Langevin, CS 20156, 38042
Grenoble Cedex 9, France  E. Nacher  Instituto de Física Corpuscular, CSIC -
Universidad de Valencia, E-46980, Valencia, Spain  A. Negret  Horia Hulubei
National Institute for R&D in Physics and Nuclear Engineering, RO-077125
Bucharest, Romania  R. D. Page  Department of Physics, Oliver Lodge
Laboratory, University of Liverpool, Liverpool L69 7ZE, United Kingdom  A.
Perea  Instituto de Estructura de la Materia, CSIC, Serrano 113 bis, E-28006
Madrid, Spain  V. Pucknell  STFC Daresbury, Daresbury, Warrington WA4 4AD,
United Kingdom  P. Rahkila  University of Jyväskylä, Department of Physics,
Accelerator laboratory, P.O. Box 35(YFL) FI-40014 University of Jyväskylä,
Finland  K. Rezynkina  KU Leuven, Instituut voor Kern- en Stralingsfysica,
B-3001 Leuven, Belgium  Université de Strasbourg, CNRS, IPHC UMR7178, F-67000,
Strasbourg, France  V. Sánchez-Tembleque  Grupo de Fìsica Nuclear, EMFTEL &
IPARCOS, Universidad Complutense de Madrid, 28040 Madrid, SPAIN  K. Schomacker
Institut für Kernphysik, Universität zu Köln, D-50937 Köln, Germany  M.
Stryjczyk  KU Leuven, Instituut voor Kern- en Stralingsfysica, B-3001 Leuven,
Belgium  University of Jyväskylä, Department of Physics, Accelerator
laboratory, P.O. Box 35(YFL) FI-40014 University of Jyväskylä, Finland  C.
Sürder  Institut für Kernphysik, Technische Universität Darmstadt, D-64289
Darmstadt, Germany  O. Tengblad  Instituto de Estructura de la Materia, CSIC,
Serrano 113 bis, E-28006 Madrid, Spain  V. Vedia  Grupo de Fìsica Nuclear,
EMFTEL & IPARCOS, Universidad Complutense de Madrid, 28040 Madrid, SPAIN  N.
Warr  Institut für Kernphysik, Universität zu Köln, D-50937 Köln, Germany

###### Abstract

Polonium isotopes having two protons above the shell closure at $Z=82$ show a
wide variety of low-lying high-spin isomeric states across the whole chain.
The structure of neutron-deficient isotopes up to 210Po ($N=126$) is well
established as they are easily produced through various methods. However,
there is not much information available for the neutron-rich counterparts for
which only selective techniques can be used for their production. We report on
the first fast-timing measurements of yrast states up to the 8+ level in
214,216,218Po isotopes produced in the $\beta^{-}$ decay of 214,216,218Bi at
ISOLDE, CERN. In particular, our new half-life value of 607(14) ps for the
8${}_{1}^{+}$ state in 214Po is nearly 20 times shorter than the value
available in literature and comparable with the newly measured half-lives of
409(16) and 628(25) ps for the corresponding 8${}_{1}^{+}$ states in
216,218Po, respectively. The measured $B(E2;8_{1}^{+}\to 6_{1}^{+})$
transition probability values follow an increasing trend relative to isotope
mass, reaching a maximum for 216Po. The increase contradicts the previous
claims of isomerism for the $8^{+}$ yrast states in neutron-rich 214Po and
beyond. Together with the other measured yrast transitions, the $B(E2)$ values
provide a crucial test of the different theoretical approaches describing the
underlying configurations of the yrast band. The new experimental results are
compared to shell-model calculations using the KHPE and H208 effective
interactions and their pairing modified versions, showing an increase in
configuration mixing when moving towards the heavier isotopes.

††preprint: APS/123-QED

Nuclear isomers are excited meta-stable states with half-lives of the order of
nanoseconds and longer. The isomerism is usually caused by large differences
between initial and final states in physical properties such as total angular
momentum, its projection on the symmetry axis ($K$-isomers), shape, seniority
or the low energy difference between the states [1, 2]. By determining the
decay energy, branching ratio of de-exciting transitions and half-life of the
state in question, one can deduce the reduced transition probability (a
measure of the decay strength) and compare it to the single-particle estimate
(which assumes that the decay involves a single nucleon). The ratio of the
experimental and single-particle values defines a hindrance (or enhancement)
factor which is related to the underlying decay mechanism. The occurrence of
isomeric states in the vicinity of doubly-magic shell closures represents one
of the important benchmarks for testing the predictive power of the nuclear
shell model using different residual nucleon-nucleon interaction [3].

In the region around the doubly-magic ${}^{208}_{82}$Pb126, low-energy nuclear
structure is often dominated by a relatively high-$j$ single-particle orbital
(with $j\geq$ 7/2). Seniority ($\nu$) represents the number of particles that
are not paired (i.e., they are not coupled to angular momentum $J=0$) and can
be regarded as a good quantum number. Seniority isomers arise due to the
selection rules associated with the conservation of $\nu$ and are encountered
in semi-magic nuclei because electric quadrupole ($E2$) transitions between
$\nu$ = 2 states ($e.g.$ $8_{1}^{+}\to 6_{1}^{+}$) are small when the valence
shell is close to half-filled [4].

Polonium isotopes, having two protons in the $h_{9/2}$ orbital above the
$Z=82$ shell closure, are especially suitable for testing the seniority scheme
across the long chain of isotopes, also crossing the $N=126$ shell closure. A
large bulk of data exists in the literature on excited states and reduced
transition probabilities $B(E2)$ in neutron-deficient 198-208Po isotopes.
These nuclei can be studied through fusion-evaporation reactions with heavy
ions [5, 6], Coulomb excitation experiments [7, 8] or $\beta^{+}/$EC decay
[9]. These studies generally agree that the yrast structure (including the
$8^{+}$ seniority isomers) in the neutron-deficient polonium isotopes is
dominated by the two-proton configuration $\pi(h_{9/2})^{2}$ coupled to the
quadrupole vibrations of the underlying even-even 208Pb core [10, 11, 12].

For isotopes heavier than 210Po, the high-$j$ neutron orbital $\nu(g_{9/2})$
starts to be filled. It is therefore expected to compete with the
$\pi(h_{9/2})^{2}$ in determining the configuration of $0^{+}-8^{+}$ yrast
states. The experimental studies of neutron-rich nuclei with $Z>82$ and
$N>126$ are more challenging, as only very specific production techniques can
be employed. One of the methods to populate excited states in such nuclei is
through multi-nucleon transfer reactions such as 18O + 208Pb, used to study
210,212,214Po [13, 14, 15, 16].

A recent study of 214Po [14] reported its yrast $8^{+}$ state half-life of
$13(1)$ ns and $B(E2;8_{1}^{+}\to 6_{1}^{+})=0.54(4)$ W.u., underlining its
isomeric character. Based on the resemblance between the excitation energies
of yrast states in 214,216,218Po, the authors proposed that similar isomers
might exist in 216,218Po, their main excitation mechanism being one-neutron-
pair breaking.

To reach isotopes beyond 214Po, high-energy fragmentation or spallation in
direct or inverse kinematics are the only methods of choice [17, 18]. In the
past, by using proton-induced spallation in 232Th and 238U at ISOLDE, CERN,
yrast states in 216,218Po up to $8^{+}$ have been populated through the
$\beta^{-}$ decay of 216,218Bi [19, 20], however, half-life measurements of
excited states were not performed.

In this Letter, we report on the first fast-timing measurements of the yrast
$2^{+},4^{+},6^{+}$ and $8^{+}$ states in 214,216,218Po populated in the
$\beta^{-}$ decay of 214,216,218Bi at ISOLDE, CERN. The present results
strongly disagree with the previous literature value of the
$T_{1/2}(8_{1}^{+})$ in 214Po [14]. The extracted $B(E2)$ values are compared
to state-of-the-art shell-model calculations in order to understand the
underlying structure of the yrast bands of 214,216,218Po leading to a
comprehensive characterization of the single-particle excitations in the
“northeast” region of 208Pb in the nuclear chart. We note that the detailed
discussion of the decay schemes and newly-identified states in 214,216Po from
the same experiment has been recently reported for the $\beta^{-}$ decay of
214Bi [21] and 216Bi [22].

The 214,216,218Bi isotopes were produced at ISOLDE-CERN [23] by bombarding a
50 g/cm2 thick UCx target with the 1.4 GeV proton beam delivered by the Proton
Synchrotron Booster. The bismuth atoms released from the target were
selectively ionized by the Resonance Ionization Laser Ion Source (RILIS),
accelerated to 50 keV and mass-separated using the ISOLDE High Resolution
Separator.

The bismuth ions were implanted on an aluminised Mylar® tape at the center of
the ISOLDE Decay Station (IDS) detection setup equipped with a fast EJ232
plastic scintillator used as a $\beta$ detector, four HPGe Clover detectors
for the detection of $\gamma$ rays in the daughter nuclei, and two small-
volume (30.5 cm3) conic LaBr3(Ce) detectors. This detection setup enabled the
measurement of nuclear half-lives using the $\beta\gamma(t)$,
$\gamma\gamma(t)$ and $\beta\gamma\gamma(t)$ fast-timing method [24, 25], a
well established technique at IDS [26, 27, 28, 29]. It covers half-lives
between 10 ps and 100 ns by using fast-timing detectors such as LaBr3(Ce) and
plastic scintillators in coincidence with HPGe detectors. More details about
the calibration procedure of the detectors is provided in the Supplemental
Material [30].

Energy spectra recorded by the HPGe and LaBr3(Ce) detectors in coincidence
with $\beta^{-}$ events are shown in Fig. 1. The production yields of
214,216,218Bi ions were measured to be $>2\times 10^{4}$ (limited in order to
reduce the count rate in the detectors), $1.5\times 10^{3}$ and $2\times
10^{2}$ ions/$\mu$C, respectively.

Figure 1:  $\beta$-gated $\gamma$-ray spectra recorded by the HPGe (black) and
LaBr3(Ce) (red) detectors following the $\beta^{-}$ decay of 214Bi (a), 216Bi
(b) and 218Bi (c). The yrast transitions in 214,216,218Po are labeled.

The half-lives of the $8_{1}^{+}$ states in 214,216,218Po were extracted
through the $\beta\gamma\gamma(t)$ method by gating on $\beta$ particles in
the plastic scintillator as the START signal and the $8_{1}^{+}\rightarrow
6_{1}^{+}$ $\gamma$-rays in the LaBr3(Ce) detectors as the STOP signal for the
Time to Amplitude Converters. An extra gating condition was required on
$\gamma$ rays in the HPGe detectors originating from any of the three
transitions below the $6_{1}^{+}$ levels ($6_{1}^{+}\to 4_{1}^{+},4_{1}^{+}\to
2_{1}^{+}$ and $2_{1}^{+}\to 0_{1}^{+}$) in order to further reduce the
background of resulting time-difference distributions, shown in Fig. 2 (a-c).

The possible delayed contributions from high-lying long-lived levels feeding
the $8_{1}^{+}$ states were investigated by HPGe gating on the $\gamma$-ray
transitions from above. In the case of 214,216Po, we report half-lives of
73(7) and 155(14) ps for the $8_{2}^{+}$ 1824.5- [21] and 1802.6-keV [19]
levels, respectively. For 218Po, no long-lived states feeding the $8_{1}^{+}$
level were identified, either because of the lack of statistics or their
absence. The LaBr3(Ce) detectors cannot resolve the transitions feeding (240
keV) and de-exciting (244 keV) the yrast 8+ state in 214Po, see Fig. 1 (a),
however, the two contributions were taken into account for the half-life
analysis. The values of the half-lives reported in Table 1 are corrected for
the contributions mentioned above and result from the convolution fit between
a Gaussian (prompt response of the fast-timing detectors), an exponential
decay and a flat background.

The time distributions for the decays of $2^{+}_{1}$, $4^{+}_{1}$, $6^{+}_{1}$
states in 214,216,218Po are shown in Fig. 2 (d-l). The deduced half-lives are
reported in Table 1 and have been measured through the $\gamma\gamma(t)$
method by extracting the timing information from the LaBr3(Ce) detectors for
$\gamma$ rays directly populating and de-populating each level.

In the case of 214Po, a significant discrepancy is observed between the newly
measured half-life value of $T_{1/2}(8_{1}^{+})=607(14)$ ps, and the previous
one of 13(1) ns [14]. The latter was measured using the slope fitting of time-
difference distributions from HPGe detectors, known for having a much slower
and less accurate time response compared to the LaBr3(Ce) detectors.

Figure 2:  (a-c) Background-subtracted time difference distributions between
$\beta$ particles recorded in the fast plastic scintillator and $\gamma$ rays
directly depopulating the yrast $8^{+}$ states in 214,216,218Po, extracted
from $\beta$-$\gamma_{La}$-$\gamma_{Ge}$ events. (d-l) Delayed (black) and
anti-delayed (red) background-subtracted time distributions between $\gamma$
rays directly feeding and de-exciting the yrast 2+, 4+ and 6+ states in
214,216,218Po, extracted from $\gamma_{La}$-$\gamma_{La}$ events. The half-
lives, reported in Table 1, were measured using the decay slope method (a-c)
and the centroid-shift method (d-l) [24, 25]. The longer lived component
(tail) in the 214Po 6+ spectrum (d) originates from the 8+ half-life due to an
overlap in the LaBr3(Ce) energy spectra between the feeding (240 keV) and de-
exciting (245 keV) transitions. However, it was excluded when evaluating the
centroid shift of the distribution.

To link the newly measured data to the underlying nuclear structure, two
effective interactions were used in the present work for the calculation of
low-energy levels and $E2$ transition rates: (i) the recently-developed H208
effective interaction [31, 3] already used in the interpretation of
212,214,216Po new experimental data [21, 32, 33, 22]; (ii) the well-
established Kuo-Herling interaction as modified by Warburton and Brown [34],
denoted as KHPE. More details regarding the calculations are provided in the
Supplemental Material [30].

Additional modifications of the two interactions, denoted by H208-m and
KHPE-m, were implemented within this work in order to reconcile calculated and
experimental $B(E2)$ values. Specifically, the diagonal and off-diagonal
$1h_{9/2}$ proton-pairing matrix elements were reduced by 100 keV in H208-m.
The KHPE-m reduction was 100 keV for 218Po and 200 keV for 210-216Po (for a
similar 200-keV reduction in 218Po, the $8_{1}^{+}$ state excitation energy
becomes lower than the $6_{1}^{+}$ one). The $1h_{9/2}0i_{13/2}$ matrix
elements were not modified, as they had an insignificant effect on the studied
$E2$ transitions.

Figure 3: The calculated $0^{+}-8^{+}$ yrast energy levels of even-even
210-218Po isotopes using H208, H208-m, KHPE and KHPE-m interactions, compared
to the available experimental data [35, 36, 37, 38, 39, 14].

The calculated and experimental excitation energies of the $0^{+}-8^{+}$ yrast
levels in 210-218Po isotopes are shown in Fig. 3. The comparison indicates
that the H208 interaction is in better agreement with the adopted experimental
data [36, 37, 38, 39] than KHPE, except for 210Po [35]. A common shortcoming
is a compression by about 200 keV observed for the $8_{1}^{+}$ state of 218Po,
which we attribute to the adopted truncation for that nucleus.

In order to investigate in detail the wave function structure of the yrast
states, Table II in the Supplemental Material [30] reports their main
components (with $\geq$ 10% contribution), where we note the following points:
(i) all the yrast $2^{+}-8^{+}$ states in 210Po ($N=126$) have nearly pure
proton $1h_{9/2}$ contribution; (ii) for $N>126$ the contribution from the
pure $1h_{9/2}$ diminishes, and the composition is more fragmented. The
neutron $1g_{9/2}$ state starts to play a strong role in all the yrast states,
$e.g.$ nearly 54% in $0^{+}-6^{+}$ states in 212Po, and with a further
admixture of $1f_{7/2}$ proton for the $8_{1}^{+}$ state. An even more
fragmented wave function can be seen for 214-218Po, with the added
contribution from the $0i_{11/2}$ neutron orbital.

The calculated $B(E2;J^{\pi}\to(J-2)^{\pi})$ reduced transition probabilities
interconnecting the yrast $0^{+}-8^{+}$ states are compared to our new results
and to available data from literature [35, 36, 14, 15, 40] in Fig. 4, Table 1
here and Table I in the Supplemental Material [30]. The results of the H208
and KHPE calculations are close to each other, with the exception of 218Po,
where different truncations are implemented. Overall, they are within error
bars with the new measurements, with few exceptions mentioned below.

|  |  | $B(E2;J^{\pi}\to(J-2)^{\pi})$ (W.u.)
---|---|---|---
Nucl. | $J^{\pi}$ |  $T_{1/2}$ (ps) | Exp. | H208 | H208 | KHPE | KHPE
|  | Exp. |  |  | -m |  | -m
| $2_{1}^{+}$ | 13(5) | 7(3) | 13.6 | 14.3 | 12.3 | 13.2
| $4_{1}^{+}$ | 35(5) | 18(3) | 16.9 | 18.5 | 13.3 | 16.5
214Po | $6_{1}^{+}$ | 118(5) | 16(1) | 11.4 | 14.7 | 5.9 | 12.6
| $8_{1}^{+}$ | 607(14) | 11.3(3) | 1.2 | 9.2 | 0.1 | 6.6
|  | 13(1) ns | 0.54(4)[14] |  |  |  |
| $2_{1}^{+}$ | 11(5) | 13(6) | 18.1 | 18.9 | 17.7 | 18.7
| $4_{1}^{+}$ | 21(5) | 26(6) | 25.2 | 27.1 | 22.1 | 25.2
216Po | $6_{1}^{+}$ | 31(5) | 37(6) | 18.8 | 25.0 | 9.8 | 23.2
| $8_{1}^{+}$ | 409(16) | 24(1) | 16.2 | 17.8 | 3.2 | 15.5
| $2_{1}^{+}$ |  $<$15 |  $>$13 | 19.2 | 20.1 | 14.8 | 15.3
| $4_{1}^{+}$ |  $<$15 |  $>$33 | 29.9 | 31.5 | 21.5 | 22.9
218Po | $6_{1}^{+}$ | 20(8) | 40(16) | 28.3 | 35.0 | 2.8 | 3.2
| $8_{1}^{+}$ | 628(25) | 7.8(3) | 8.5 | 16.2 | 1.0 | 0.003
Table 1:  Experimental $T_{1/2}$ and $B(E2)$ values in 214-218Po measured in
the present work and Ref. [14] (bold), compared to calculated $B(E2)$s using
various effective interactions: H208, KHPE, and their pairing-modified
versions.  Figure 4:  $B(E2;J^{\pi}\rightarrow(J-2)^{\pi})$ values for
transitions between 214,216,218Po yrast states measured in the present work
(full circles with uncertainties) and from adopted experimental values [35,
36] or the most recent experimental results [14, 15, 40]. The comparison with
shell-model calculations using the H208 interaction initial and pairing-
modified versions (colored lines) is presented. Lower limits for the
experimental $B(E2)$ values are indicated with an arrow pointing up. All the
values are reported in Table 1.

We first note the over-estimation of the $B(E2,2_{1}^{+}\to 0_{1}^{+})$ values
in 210,212,214Po isotopes by a factor of two. However, the calculated
$B(E2;4_{1}^{+}\to 2_{1}^{+})$ and $B(E2;6_{1}^{+}\to 4_{1}^{+})$ values are
very well reproduced by the shell-model calculations, shown in Table 1, which
suggests that the correct structure of the $2_{1}^{+}$ state wave function was
taken into account. This $B(E2,2_{1}^{+}\to 0_{1}^{+})$ discrepancy was
already explained in Refs. [34, 41] by the absence of the contribution from
particle-hole core excitations in the $0_{1}^{+}$ wave function and later, the
$\alpha$-cluster model was proposed to address this issue [40].

The most striking discrepancy is observed for the $B(E2;8_{1}^{+}\to
6_{1}^{+})$ values in 214Po, where our newly measured value is at least 10
times larger than the original shell-model calculated one, questioning the
isomeric nature of the $8_{1}^{+}$ state in 214Po. To understand the origin of
this inconsistency between theory and our new measurements, the modified
effective interactions were employed. As a consequence, there was a remarkable
increase of the calculated $B(E2;8_{1}^{+}\to 6_{1}^{+})$ value in 214Po, from
1.2 to 9.2 W.u. in H208-m and from 0.1 to 6.6 W.u. in KHPE-m, close to the
measured value of 11.3(3) W.u., as shown in Fig. 4 and Table 1.

To provide a complementary insight into the effect of the aforementioned
pairing modifications on the wave functions structure of the lowest states in
214Po, the same calculations are performed in the seniority $(\nu)$ scheme
using the $j$-coupled Nathan code [42]. The results obtained using H208 and
H208-m are displayed in Fig. 5, where $T(x,y)$ represents the number of
neutron ($x$) or proton ($y$) pairs being broken ($e.g.$ the $0_{1}^{+}$ state
wavefunction is dominated by a 66% contribution from $T(0,0)$ which
corresponds to $\nu=0$, and the $2_{1}^{+}$ by $45\%$ of $T(1,0)$
corresponding to $\nu=2$). No important change between the two versions can be
observed in the seniority structure of the $0_{1}^{+},2_{1}^{+},4_{1}^{+}$ and
$6_{1}^{+}$ states. The effect of reducing the $1h_{9/2}$ proton pairing
matrix elements is apparent on the structure of the two lowest $8^{+}$ states.
In H208, the main wave function component of the $8_{1}^{+}$ state was
characterized by seniority $\nu=2$ (one neutron pair broken). The $8_{2}^{+}$
state was initially predicted very close in energy to the $8_{1}^{+}$ state,
connected by a strong $E2$ transition to the $6_{1}^{+}$ state and
characterized by a mixing of $\nu$=2, 4 and 6. Reducing the pairing in H208-m
leads to an inversion of the $8_{1,2}^{+}$ states, the $\pi(1h_{9/2}1f_{7/2})$
state comes below the $\nu g_{9/2}^{4}$ $8^{+}$ one, resulting in a good
agreement between experiment and calculations.

The increasing trend of the experimental $B(E2;J^{\pi}\to(J-2)^{\pi})$ values
for transitions from the $J=2,4,6,8$ states, shown in Fig. 4, was predicted
theoretically in neutron-rich Po isotopes up to $A=216$, with the exception of
a dip in 214Po, for the $8_{1}^{+}\to 6_{1}^{+}$ transition (see Figs. 6-9
from Ref. [3]). The trend was explained through the increase of the
collectivity and quadrupole correlations with respect to the neutron number,
leading to strong electromagnetic strengths.

The sudden decrease of the experimental $B(E2;8_{1}^{+}\to 6_{1}^{+})$ for
218Po is also reproduced theoretically, despite the adopted truncation in this
nucleus. This can be attributed to the contribution from $i_{13/2}$ and
$j_{15/2}$ proton and neutron orbitals in the structure of the $8_{1}^{+}$ and
$6_{1}^{+}$ states. A strong transition probability of 27 W.u. is calculated
when excluding the aforementioned orbitals. Their inclusion determines a
reduction of quadrupole correlations due to pairing, leading to the
$B(E2;8_{1}^{+}\to 6_{1}^{+})$ decrease observed in 218Po.

Figure 5: The wave-function components of 214Po states calculated in seniority
scheme where $T(x,y)$ represents the number of neutron ($x$) or proton ($y$)
pairs being broken. The initial and pairing-modified versions of the H208
effective interaction are employed within the Nathan code.

In conclusion, the $B(E2)$ transition probabilities of yrast states up to the
$8^{+}$ level in 214,216,218Po isotopes populated in the $\beta^{-}$ decay of
214,216,218Bi are established for the first time through fast-timing
measurements at the ISOLDE Decay Station. The comparison with shell-model
calculations using the standard versions of the H208 and KHPE effective
interactions shows reasonable agreement for the lower-lying transitions. A
significant disagreement was observed in the case of the $8_{1}^{+}\to
6_{1}^{+}$ transitions, especially for 214Po where a much higher transition
probability was measured ($B(E2)=11.3(3)$ W.u. instead of 0.54(4) W.u. [14]),
disproving its previously-claimed isomeric character. Similarly, based on the
large $B(E2)$ values, the $8_{1}^{+}$ states in 216-218Po should not be
considered isomeric. An improved agreement between theory and the present
measurement was achieved after reducing the pairing strength of the
interactions, effectively inverting the lowest predicted $8^{+}$ states and
confirming the two-proton configuration $\pi(1h_{9/2}1f_{7/2})$ of the yrast
$8^{+}$ state dominated by quadrupole correlations. The new findings can serve
as a guideline for future analysis, providing a clear indication of the nature
of $8_{1}^{+}$ states in neutron-rich polonium isotopes. Additionally, they
serve as an extremely valuable input for the different theoretical approaches
and constitute a stringent test for the effective interactions.

We would like to thank B.A. Brown, A. Gargano and the late H. Grawe for the
inspiring discussions and the ISOLDE staff for the support provided during the
experiment. We acknowledge support for this work by the Romanian IFA grant
CERN/ISOLDE and Nucleu project No. PN 23 21 01 02, the Slovak Research and
Development Agency (Contract No. APVV-22-0282) and the Slovak grant agency
VEGA (Contract No. 1/0651/21), the Research Foundation Flanders (FWO,
Belgium), BOF KU Leuven (C14/22/104), fWO and F.R.S.-FNRS under the Excellence
of Science (EOS 40007501) programme, the Spanish MCIN/AEI/10.13039/
501100011033 under grants RTI2018-098868-B-I00 and PID2021-126998OB-I00, the
Academy of Finland project No. 354968, the German BMBF under contract
05P21PKCI1 and Verbundprojekt 05P2021, and the United Kingdom Science and
Technology Facilities Council (STFC) through the grants ST/P004598/1 and
ST/V001027/1.

## References

  * Garg  _et al._ [2023] S. Garg _et al._ , Atomic Data and Nuclear Data Tables 150, 101546 (2023).
  * Walker and Dracoulis [1999] P. Walker and G. Dracoulis, Nature 52, 35 (1999).
  * Naïdja [2021] H. Naïdja, Phys. Rev. C 103, 054303 (2021).
  * Isacker [2011] P. V. Isacker, J. Phys.: Conf. Ser. 322, 012003 (2011).
  * Maj  _et al._ [1990] A. Maj _et al._ , Nuclear Physics A 509, 413 (1990).
  * Stoychev  _et al._ [2023] K. Stoychev _et al._ , Phys. Rev. C 108, 014316 (2023).
  * Grahn  _et al._ [2016] T. Grahn _et al._ , Eur. Phys. J. A 52, 340 (2016).
  * Wrzosek-Lipska and Gaffney [2016] K. Wrzosek-Lipska and L. P. Gaffney, J. Phys. G: Nucl. Phys. 43, 024012 (2016).
  * Bijnens  _et al._ [1995] N. Bijnens _et al._ , Phys. Rev. Lett. 75, 4571 (1995).
  * Häusser  _et al._ [1976] O. Häusser _et al._ , Nuclear Physics A 273, 253 (1976).
  * Neyens  _et al._ [1997] G. Neyens _et al._ , Nuclear Physics A 625, 668 (1997).
  * Dracoulis  _et al._ [2016] G. D. Dracoulis, P. M. Walker, and F. G. Kondev, Reports on Progress in Physics 79, 076301 (2016).
  * Astier  _et al._ [2010] A. Astier _et al._ , Eur. Phys. J. A 46, 165 (2010).
  * Astier and Porquet [2011] A. Astier and M.-G. Porquet, Phys. Rev. C 83, 014311 (2011).
  * Kocheva  _et al._ [2017a] D. Kocheva _et al._ , Eur. Phys. J. A 53, 175 (2017a).
  * Kocheva  _et al._ [2017b] D. Kocheva _et al._ , Phys. Rev. C 96, 044305 (2017b).
  * Benzoni  _et al._ [2012] G. Benzoni _et al._ , Phys. Lett. B 715, 293 (2012).
  * Morales  _et al._ [2014] A. I. Morales _et al._ , Phys. Rev. C 89, 014324 (2014).
  * Kurpeta  _et al._ [2000] J. Kurpeta _et al._ , Eur. Phys. J. A 7, 49 (2000).
  * De Witte  _et al._ [2004] H. De Witte _et al._ , Phys. Rev. C 69, 044305 (2004).
  * Andel  _et al._ [2021] B. Andel _et al._ (IDS Collaboration), Phys. Rev. C 104, 054301 (2021).
  * Andel  _et al._ [2024] B. Andel _et al._ (IDS Collaboration), Phys. Rev. C 109, 064321 (2024).
  * Borge and Jonson [2017] M. J. G. Borge and B. Jonson, J. Phys. G: Nucl. Phys. 44, 044011 (2017).
  * Moszynski and Mach [1989] M. Moszynski and H. Mach, NIM A 277,  (1989).
  * Mach  _et al._ [1989] H. Mach, R. L. Gill, and M. Moszynski, NIM A 280, 49 (1989).
  * Fraile [2017] L. M. Fraile, J. Phys. G: Nucl. Phys. 44, 094004 (2017).
  * Lică  _et al._ [2016] R. Lică _et al._ (IDS Collaboration), Phys. Rev. C 93, 044303 (2016).
  * Lică  _et al._ [2017] R. Lică _et al._ , J. Phys. G: Nucl. Phys. 44, 054002 (2017).
  * Lică  _et al._ [2018] R. Lică _et al._ (IDS Collaboration), Phys. Rev. C 97, 024305 (2018).
  * [30] See Supplemental Material at [URL will be inserted by publisher] for additional details about the detector calibration procedure and shell-model calculations.
  * Naïdja [2018] H. Naïdja, Physica Scripta 94, 014005 (2018).
  * Fernández  _et al._ [2021] A. Fernández _et al._ , Phys. Rev. C 104, 054316 (2021).
  * Karayonchev  _et al._ [2022] V. Karayonchev _et al._ , Phys. Rev. C 106, 064305 (2022).
  * Warburton and Brown [1991] E. K. Warburton and B. A. Brown, Phys. Rev. C 43, 602 (1991).
  * Shamsuzzoha Basunia [2014] M. Shamsuzzoha Basunia, NDS 121, 561 (2014).
  * Auranen and McCutchan [2020] K. Auranen and E. McCutchan, NDS 168, 117 (2020).
  * Zhu and McCutchan [2021] S. Zhu and E. McCutchan, NDS 175, 1 (2021).
  * Wu [2007] S.-C. Wu, NDS 108, 1057 (2007).
  * Singh  _et al._ [2019] B. Singh _et al._ , NDS 160, 405 (2019).
  * von Tresckow  _et al._ [2021] M. von Tresckow _et al._ , Phys. Lett. B 821, 136624 (2021).
  * Karayonchev  _et al._ [2019] V. Karayonchev _et al._ , Phys. Rev. C 99, 024326 (2019).
  * [42] E. Caurier and F. Nowacki, unpublished .