# Highly depleted alkali metals in Jupiter’s deep Atmosphere

Ananyo Bhattacharya Ananyo Bhattacharya<EMAIL_ADDRESS>University of
Michigan, Ann Arbor  Cheng Li  University of Michigan, Ann Arbor  Sushil K.
Atreya  University of Michigan, Ann Arbor  Paul G. Steffes  Georgia Institute
of Technology  Steven M. Levin  NASA Jet Propulsion Laboratory  Scott J.
Bolton  Southwest Research Institute  Tristan Guillot  Universite C´ote d’Azur
Pranika Gupta  University of Michigan  Andrew P. Ingersoll  California
Institute of Technology  Jonathan I. Lunine  Cornell University  Glenn S.
Orton  NASA Jet Propulsion Laboratory  Fabiano A. Oyafuso  NASA Jet Propulsion
Laboratory  J. Hunter Waite  Southwest Research Institute  Waite Science LLC
Amadeo Belloti  Michael H. Wong  University of California, Berkeley

###### Abstract

Water and ammonia vapors are known to be the major sources of spectral
absorption at pressure levels observed by the microwave radiometer (MWR) on
Juno. However, the brightness temperatures and limb darkening observed by the
MWR at its longest wavelength channel of 50 cm (600 MHz) in the first 9
perijove passes indicate the existence of an additional source of opacity in
the deep atmosphere of Jupiter (pressures beyond 100 bar). The absorption
properties of ammonia and water vapor, and their relative abundances in
Jupiter’s atmosphere do not provide sufficient opacity in the deep atmosphere
to explain the 600 MHz channel observation. Here we show that free electrons
due to the ionization of alkali metals, i.e. sodium, and potassium, with sub-
solar metallicity, [M/H] (log based 10 relative concentration to solar) in the
range of [M/H] = -2 to [M/H] = -5 can provide the missing source of opacity in
the deep atmosphere. If the alkali metals are not the source of additional
opacity in the MWR data, then their metallicity at 1000 bars can only be even
lower. This upper bound of -2 on the metallicity of the alkali metals
contrasts with the other heavy elements – C, N, S, Ar, Kr, and Xe – which are
all enriched relative to their solar abundances having a metallicity of
approximately +0.5.

Solar System (1528) – Chemical abundances(224) – Jupiter(873) – Extrasolar
gaseous giant planets(509)

††software: High-performance Atmospheric Radiation Package (HARP) [Li et al.
(2018b), Bhattacharya et al. (2023)]

##  1 Introduction

The alkali metals sodium and potassium have been previously detected in the
atmospheres of hot Jupiters and a super-Neptune together with lithium [Chen et
al. (2018)] in the latter. The detections show a large range of abundances
from highly substellar to super-stellar values [Welbanks et al. (2019), Demory
et al. (2011)]. Alkali metal abundances are important in understanding the
formation of hot Jupiters and represent a bridge between the refractory and
volatile elements, which in molecular form seed the growth of planets.
Obtaining the abundance of alkali metals in Jupiter can potentially serve as a
first constraint on the ratio of rocky to icy material in the interior of the
solar system’s largest planet when combined with the elemental and molecular
abundances provided by the Galileo Probe Mass Spectrometer (GPMS) [Atreya et
al. (1999), Wong et al. (2004), Atreya et al. (2019)] and Juno constraints on
water [Li et al. (2020)]. Here we derive observationally based abundances of
alkali metals in Jupiter’s atmosphere to determine whether they are enriched
relative to solar like the other heavy elements or depleted.

To obtain these abundances requires knowing the deep structure of Jupiter’s
atmosphere. The shallower part of Jupiter’s atmosphere has been previously
investigated at microwave frequencies by the Very Large Array (VLA) telescope
[de Pater & Dunn (2003), de Pater et al. (2019)]. VLA probes Jupiter at
frequencies in the range of 74 MHz to 50 GHz [de Pater et al. (2019)].
However, confusion from Jupiter’s powerful synchrotron radiation does not
allow VLA to observe Jupiter’s atmosphere below 5 GHz [de Pater & Dunn
(2003)], limiting its reach to less than 5 bars, leaving the deep atmosphere
of Jupiter inaccessible from microwave and radio frequency observatories from
Earth. The orbit of Juno and the spin of the spacecraft allow the spacecraft
to make observations at low frequencies, i.e. 0.6 GHz and 1.2 GHz, by avoiding
the energetic electron belts around Jupiter from its field of view. Access to
greater depths allows for the investigation of bulk elemental abundances of N
and O in Jupiter [Janssen et al. (2017), Bolton et al. (2017), Steffes et al.
(2017)].

The Microwave Radiometer (MWR) instrument onboard the Juno orbiter is a
passive radiometer that is designed to measure the internal heat emitted by
Jupiter’s atmosphere at six different frequencies ranging from 0.6 GHz to 22
GHz [Janssen et al. (2017)]. The brightness temperature measured by MWR at
these frequencies sounds different levels of Jupiter’s atmosphere
corresponding to pressures from 0.3 bar to 250 bar [Janssen et al. (2017)]. In
addition, the highly inclined polar orbit and rotation of the Juno spacecraft
aided in the high spatial resolution necessary for probing Jupiter’s
atmosphere at various latitudes [Bolton et al. (2017)].

Previous analysis of the MWR data at the 0.6 GHz found an unanticipated limb-
darkening signal, which cannot be explained by nominal absorbers such as
ammonia and water [Li et al. (2020)]. Based on investigation of thermodynamic
models of Jupiter’s deep atmosphere between 50 bar and 1 kbar [Fegley Jr &
Lodders (1994), Weidenschilling & Lewis (1973)] we conjecture that the free
electrons from thermally ionzied alkali metals may provide the missing
opacity. Alkali metals are expected to undergo condensation to form clouds in
the deep atmosphere [Visscher et al. (2006), Morley et al. (2012)]. Na2S and
KCl are the first chemical species to condense in the above pressure range and
thereby act as a sink for atomic sodium and potassium [Fegley Jr & Lodders
(1994)]. Furthermore, high-temperature environments cause alkali metals to
undergo ionization due to their low ionization energies [Bagenal et al.
(2007)]. Density and temperature play a role in governing the electron
densities according to the Saha ionization equation (Eq. 2). Electrons
generated from alkali metal ionization act as a source of absorption at
microwave frequencies that could affect the brightness temperatures at the 0.6
GHz frequency channel. Therefore, the objective of this study is to determine
the alkali metal abundance in the deep atmosphere of Jupiter.

To facilitate comparison of our results on alkali metals with those of the
extrasolar planets we express the abundances of non- hydrogen and helium
elements using astronomical terminology, e.g., metallicity. The metallicity
(_[M/H]_) of an element is the logarithm of the ratio of elemental abundance
in a system to the stellar (or solar, for the solar system) elemental
abundance. Generally, the metallicity of a star is defined in terms of the
ratio of the number of Fe atoms to the number of hydrogen atoms. Here we
define the metallicity in terms of alkali metal abundance in Jupiter to that
of Sun e.g. for potassium, _[K/H]_ = log10(_N K/NH_)Jupiter \- log10(_N
K/NH_)Sun. For the giant planets, iron and silicon is not measurable,
emphasizing the importance of proxy indicators such as the alkali metals along
with other elements measured by Galileo probe.

##  2 Methods

Brightness temperatures from 9 perijoves i.e. PJ 1,3-9, 12 have been taken
into consideration for this article. Variations in brightness temperatures
have been observed across the planetocentric latitudes from pole-to-pole at
0.6 and 1.2 GHz channels. These variations can be attributed to various
sources of origin from the atmosphere and space environment. The most
important sources of the observed variability are (i) changes in atmospheric
structure and composition, (ii) Jupiter’s synchrotron radiation in the
microwave band, and (iii) variation in acceleration due to gravity due to the
non-spherical shape of Jupiter. The latter sources, i.e. synchrotron and
gravity need to be taken into account for proper interpretation of MWR
observations. It will aid in investigating the true variability in Jupiter’s
deep atmosphere.

The contribution of Jupiter’s gravity can be corrected by taking into account
the non-spherical shape of Jupiter. Brightness temperatures are corrected
using a gravity correction factor defined as the ratio of theoretical _T b_ at
a given latitude to that at the equator of Jupiter taking into consideration
the acceleration due to gravity at the latitude. Thereby, it transforms the
Juno observations at each latitude for equatorial gravity, which effectively
removes variation in _T b_ due to changes in Jupiter’s gravity from the
equator to the poles.

Energetic electrons in Jupiter’s space environment contribute to the
synchrotron radiation [de Pater & Dunn (2003), Levin et al. (2001), Santos-
Costa et al. (2017)]. The signature of the emission is observed in MWR data
across all the perijoves which leads to anomalous changes in _T b_. Data at
extremely high latitudes are polluted by synchrotron emission and thus, remain
of no use for investigating Jupiter’s deep atmosphere. Therefore, we only
consider the MWR data between -60 to 60 deg. latitude. The correction for
synchrotron and other sources of anomalous _T b_ is done by filtering the data
at 0.6 and 1.2 GHz for each perijove. The process is carried out by sorting
the deviations of _T b_ from the least value of Tb in a group and removing the
values greater than a filter cutoff temperature of the order of 2 K.

##  3 Results

###  3.1 Sources of Microwave Opacity

The weighting function of Jupiter’s atmospheric absorption and emission at a
given microwave frequency determines the contribution of each region in the
atmosphere to the observed brightness temperature at the given frequency. The
peak structure of the weighting function gives the range of pressure levels
corresponding to the measurements. The weighting function can be expressed as
a function of microwave opacity of the atmosphere (1). Here, _T b_ is the
brightness temperature, _W(p)_ is the weighting function as a function of
pressure, and _T(p)_ is the physical temperature profile of the atmosphere.

$T_{b}=\int_{-\infty}^{\infty}W(p)T(p)dlnP$ (1)

Fig. 1 shows the relative weighting functions, i.e. weighting function divided
by the maximum value of the function, at 0.6 GHz and 1.2 GHz with and without
alkali metals. In the absence of alkali metals, the relative weighting
functions peak at 100 bar and 30 bar, respectively[Janssen et al. (2017)]. At
0.6 GHz, the relative weighting function extends to the deeper atmosphere
below the 100 bar level, and therefore, the _T b_ derived using this channel
is sensitive to the sources of microwave opacity present in the deep
atmosphere at _p_ greater than 100 bar. The relative weighting function at 0.6
GHz channel shows a broad shape with a second maxima at kbar pressure levels
which is attributed to the increase in mass absorption coefficients of water
vapor with pressure. The mass absorption coefficient of ammonia decreases
after a maximum near 1000 bar and eventually water vapor dominates the opacity
in the deep atmosphere. Moreover, the inclusion of free electrons as sources
of opacity due to alkali metal ionization causes a decrease in the value of
the relative weighting function at 0.6 GHz around 100 bar, and a global
maximum in the relative weighting function emerges at $\sim$ 1 kbar pressure
(magenta line). The shift of the global maximum can be attributed to the
increase in opacity from free electrons with pressure as the ionization
fraction of alkali metals increases with temperature under thermal equilibrium
conditions [Saha (1920)] (described later in this section). Inclusion of lower
amounts of alkali metals ([M/H] = -5) will lead to a peak at deeper levels
(Fig. 1). However as the metallicity is increased to solar, the maximum drifts
towards lower pressures around 1 kbar level. This could be attributed to the
fact that higher abundance of alkali metals can produce higher amount of
electrons at relatively lower pressures (magenta line), whereas low abundance
of alkali metals in Jupiter would need to reach higher pressure ($>$ 1 kbar)
to produce equivalent opacity (blue line). Thereby the abundance of alkali
metals directly affects the shape of weighting function.

Figure 1: Relative weighting functions at 0.6 GHz (black) and 1.2 GHz (gray)
for a Jupiter adiabat considering the Hanley model [Hanley et al. (2009)] for
NH3 absorption. The functions peak at 100 bar and 30 bar at 0.6 GHz and 1.2
GHz respectively without the inclusion of alkali metals. The inclusion of
alkali metals (orange, magenta and blue) decreases the relative weighting
function at $\sim$ 100 bar and produces a second peak that is observed at
$\sim$ 1 kbar pressure due to the opacity contributed by free electrons from
alkali metal ionization. As the metallicity of alkali metals increase, the
global maximum of weighting function shifts towards lower pressure.

The main sources of microwave opacity at 0.6 GHz and 1.2 GHz are ammonia,
water vapor, free electrons, and collision-induced absorption by hydrogen and
helium. Hydrogen-hydrogen and hydrogen-helium collisions are the dominant
sources of collision-induced absorption processes in Jupiter. Their magnitude
is well constrained due to the invariance of hydrogen and helium abundances in
Jupiter’s deep atmosphere. The microwave absorption behavior of water and
ammonia vapor has been investigated by laboratory experiments that show the
pressure and temperature dependence of mass absorption coefficients [Devaraj
et al. (2014), Karpowicz & Steffes (2011), Bellotti et al. (2016)]. In
addition, hydrogen, methane, and water vapor contribute to line broadening in
the ammonia vapor absorption. The models based on laboratory experiments show
significant divergent behavior when extrapolated to pressures greater than 50
bar and 550 K [Bellotti et al. (2016)]. In order to obtain a robust estimate
of the range of absorption coefficients at higher temperatures, we test a grid
model describing a power scaling relationship with temperature based on the
Hanley et al. (2009) model of ammonia absorption. For water vapor absorption
at microwave frequencies, the laboratory models show divergence by orders of
magnitude. However, recent laboratory measurements [Steffes et al. (2023)] at
high pressure show that water vapor absorption can be explained by the
Bellotti et al. (2016) model. Therefore, Bellotti et al. (2016) model is
chosen to compute the water vapor opacity which incorporates water opacity
measurements at high temperatures above 500 K.

Free electrons in the atmosphere can act as a source of opacity at microwave
wavelengths through the process of free-free absorption in which electrons
absorb photons during collisions with other ions and electrons. Electrons can
be generated by the ionization of various elemental and molecular species in
the atmosphere. Due to their low ionization energies, alkali metals i.e. Na, K
are expected to be the major sources of free electrons in the atmosphere
[Heays et al. (2017)]. In Jupiter’s atmosphere, the pressure and temperatures
corresponding to the transition between the alkali metals and their compounds
are calculated using an equilibrium cloud condensation model (ECCM) [Atreya et
al. (1999), Weidenschilling & Lewis (1973)] for Jupiter’s adiabat with
saturation vapor pressures of Na2S and KCl [Visscher et al. (2006), Morley et
al. (2012)]. The condensation of alkali metals at solar abundance [Figure 2]
takes place at 352 bar for KCl and 796 bar for Na2S, with corresponding
temperatures of 967 K and 1234 K, respectively, assuming thermodynamic
equilibrium. The condensation of Na2S at deeper levels, and a higher solar
abundance of Na compared to K [Asplund et al. (2009)] will cause Na2S clouds
to be significantly more massive than KCl clouds. Thermochemical equilibrium
models indicate formation of metal hydrides and hydroxides in gas phase,
however they are much lower in abundance [Fegley Jr & Lodders (1994)] as
compared to the condensates, thereby they will not act as the primary sink of
alkali metals in Jupiter. Condensation of the alkali metal compounds occurs
when the partial pressure of a compound exceeds its saturation vapor pressure.
If condensation occurs, it causes depletion in the alkali metal abundances at
altitudes above the condensation level.

Figure 2: Condensation curves of NH3, H2O, H2S and alkali metals Na2S and KCl
at 1X solar abundance. Our calculations are based on the equilibrium cloud
condensation model [Atreya et al. (1999)], and saturation vapor pressure
corresponding to Na2S and KCl [Visscher et al. (2006), Morley et al. (2012)].
The cloud bases are at the levels where the condensation curves cross the
adiabat considering T1bar = 166.1 K [Seiff et al. (1998)].

At high pressures 100 bar and beyond, alkali metals would undergo ionization
to form cold plasma, and the electrons generated in the process would act as
an additional source of opacity at microwave frequencies. The number density
of free electrons due to the ionization of alkali metal atoms in the gas phase
and is calculated using the Saha ionization [Saha (1920)] (Eq. 2) equation
assuming Jupiter’s atmosphere to be in a state of thermal equilibrium. The
ionization equation itself assumes a single component gas phase system.
Thereby, we add the electron densities from ionization of sodium and potassium
to determine total number density of free electrons. Here, _N e_ is the
electron density, _N_ is number density, _$\epsilon$_ is ionization energy,
_$\lambda$_ is De Broglie wavelength, _g 0_ and _g 1_ are statistical weights,
_k B_ is Boltzmann constant, _m e_ is mass of the electron and _h_ is Planck’s
constant.

$\frac{N_{e}^{2}}{N-N_{e}}=\frac{2}{\lambda^{3}}\frac{g_{1}}{g_{0}}e^{-\epsilon/k_{B}T}$
(2)

$\lambda=\sqrt{\frac{h^{2}}{2\pi m_{e}k_{B}T}}$ (3)

The brightness temperatures correspond to electromagnetic radiation traveling
from the interior of Jupiter radially outwards through the atmospheric layers.
Thus, the transmission through the deep atmosphere is similar to the
transmission through a cold plasma medium. The refractive index of microwaves
propagating through a cold plasma media can be described by the Appleton-
Hartree equation [Helliwell (2014)]. The formulation is applicable to low-
temperature plasma medium both in the presence or absence of magnetic fields.
At 100-1000 bar pressure levels, the contribution of the magnetic field is
insignificant in the Appleton-Hartree formulation [Helliwell (2014)].
Therefore, a simplified version of the Appleton-Hartree equation (Eq. 4) is
used to calculate the complex refractive index of the deep atmosphere using
the electron number density calculated from the Saha ionization equation. For
an unmagnetized cold plasma medium i.e. Jupiter’s deep atmosphere, the
Appleton- Hartree equation is simplified to:

$n^{2}=1-\frac{X}{1-iZ}$ (4)

$\alpha=\frac{2\pi}{\lambda_{ch}Q}$ (5)

Here, _X_ = $\frac{\omega_{0}^{2}}{\omega^{2}}$, _Z_ = $\frac{\nu}{\omega}$ ,
$\omega_{0}$ is electron plasma frequency, $\omega$ is the angular frequency
of microwave radiation, $\omega_{h}$ is electron gyro frequency, $\nu$ is
electron- neutral collision frequency, $\lambda_{ch}$ is the frequency of a
given MWR channel, _n_ is the refractive index, $\alpha$ is the extinction
coefficient and _Q_ is the quality factor i.e. the ratio of squares of real
and imaginary parts of the refractive index.

###  3.2 Radiative Transfer Modeling

In order to draw a comparison between the MWR observations and theoretical
knowledge of Jupiter’s atmosphere, a benchmark model for the ideal Jupiter
atmosphere is constructed using a moist hydrostatic adiabat following the
ideal gas law [Li et al. (2018b), Li et al. (2018a)]. The specific heat of
hydrogen is estimated from the mixing ratio of ortho and para hydrogen
assuming thermal equilibrium between the ortho and para states. Moreover, the
temperature profile of Jupiter’s atmosphere is constructed for two cases of
reference temperatures: (i) _T_ = 166.1 K at the 1-bar pressure level from the
Galileo probe [Seiff et al. (1998)] and (ii) _T_ = 168.8 K at the 1-bar
pressure level based on the reanalysis of the Voyager radio occultation
experiment at Jupiter [Gupta et al. (2022)]. Ammonia and water vapor are
considered vapors for the moist adiabat and their partial pressure is
controlled by the cloud condensation process by forcing the partial pressures
to be equal to their saturation vapor pressures. In the deep atmosphere of
Jupiter, water and ammonia are not expected to form clouds; however, alkali
metals are expected to undergo condensation. Therefore, a similar approach is
applied to alkali metals to estimate the concentration of alkali metals
present in the gas phase available for the ionization process.

Spectral radiance is proportional to the physical temperature of the
atmosphere in the Rayleigh-Jeans limit. For microwave frequencies, we compute
the brightness temperature (_T b_) from the physical temperature using Eq.
(1). The opacity of Jupiter’s atmosphere is the sum of opacities from
individual sources discussed in the previous section i.e. ammonia, water, free
electrons, and collision-induced absorption. The abundances of ammonia and
water vapor have been assumed to be 2.7 and 5 times the solar abundance [Li et
al. (2020), Li et al. (2017)]. Because there is no a priori information on the
alkali metal abundance in Jupiter, Therefore, we compare two cases, one
without alkali metals (baseline) and another with alkali metals (treatment) in
order to provide a comparison between our current knowledge of Jupiter and MWR
data.

The spatial resolution of MWR data also provides the limb darkening
coefficient at six microwave frequencies. Limb darkening (_L d_) is defined as
the percent change in _T b_ at a given viewing angle relative to _T b_ at a
position looking vertically down to the planet center i.e. nadir. For our
simulations, we compute the limb darkening at a 45- degree angle from the
nadir. The MWR channels at 0.6 GHz and 1.2 GHz are chosen to provide a
comparison between theory and observations at higher pressures using _T b_ and
_L d_ as the observables for comparison. The benchmark case of the ideal
Jupiter atmosphere is compared with MWR observations as a function of latitude
between -40 and 40 degrees planetocentric latitude. Data from higher latitudes
are neglected due to the presence of signatures from synchrotron radiation
that is inseparable from the atmospheric contribution.

Figure 3: Limb darkening and brightness temperature MWR observations compared
with simulation results at 0.6 GHz and 1.2 GHz corresponding to Jovian
adiabats at (i) _T 1bar_ = 166.1 K and (ii) _T 1bar_ = 168.8 K , (a) _T b_ vs.
latitude at 0.6 GHz, (b) _L d_ vs. latitude at 0.6 GHz, (c) _T b_ vs. latitude
at 1.2 GHz, (d) _L d_ vs. latitude at 1.2 GHz.

A latitudinal variation in brightness temperatures is observed at both 0.6 and
1.2 GHz (Figure 3, panels (a) and (c)). The small- scale variations in _T b_
and _L d_ in all the panels can be attributed to variations in the atmospheric
temperature structure and composition. It is important to note that the
baseline case (without alkali metals) corresponds to two different temperature
profiles of Jupiter’s atmosphere for two different _T 1bar_. There is an
agreement between the baseline case and observations at 1.2 GHz in the
equatorial region (panel (c)). On the other hand, brightness temperatures at
0.6 GHz are lower than the baseline case by 40-60 K at all latitudes (panel
(a)) indicating the possibility of an additional source of opacity. Such a
source is also supported by a depressed _L d_ observed by MWR; it is 4 percent
less than the _L d_ magnitude of the ideal Jupiter atmosphere across all
latitudes (panel (b)). The mismatch between the baseline and observations at
0.6 GHz is much greater than the uncertainty in measurements and variations in
_T b_ and _L d_. Since the brightness temperatures correspond to different
pressure regions in the atmosphere, the anomalous observations at 0.6 GHz must
be attributed to the presence of an additional opacity source in the deep
atmosphere or to a different opacity source that absorbs more effectively at
0.6 GHz than at 1.2 GHz. We test four confounding factors: (1) the
distribution of ammonia, (2) the ammonia opacity at temperatures exceeding the
range of laboratory measurements, (3) the opacity of water at high
temperatures and (4) the contribution of alkali metals. The theoretical
brightness temperature and limb darkening at 0.6 GHz and 1.2 GHz is shown in
Fig. 3.

The latitudinal distribution of brightness temperatures and limb darkening
from the forward model indicates the decrease in limb darkening from the
equator to the pole at 0.6 GHz. It is opposite to the variation of limb
darkening at 1.2 GHz across the latitudes. This effect could be attributed to
the free electrons in the deep atmosphere which could be inferred from the
shift in the contribution functions toward higher pressures in presence of
alkali metals (Fig. 1). Alkali metals greatly affect the absorption behavior
at 0.6 GHz which dominates the effect of gravitation on limb darkening.

###  3.3 Ammonia, Water and Alkali Metals

Brightness temperature variations with latitude and the spectral inversion of
brightness temperatures show a non-uniform distribution of ammonia vapor in
Jupiter’s atmosphere in the deep atmosphere region [Li et al. (2017),
Ingersoll et al. (2017)]. Therefore, the non-uniform distribution of ammonia
could contribute to variations in microwave opacity of the deep atmosphere. In
order to estimate the effect of ammonia concentration variations, we perturb
the ammonia profile in the model and use a scaling factor to vary the
magnitude of ammonia vapor concentration in the model as described in Eq. (6).

$q_{NH_{3}}(P)=q_{NH_{3},0}(P)-(q_{NH_{3},0}(P)-q_{NH_{3},MWR}(P))s$ (6)

Here, _q ${}_{NH_{3}}$_ is the ammonia mass mixing ratio at a given pressure
_P_ , _q ${}_{NH_{3},0}(P)$_ is the homogeneous ammonia mixing ratio which is
set to 2.7 times solar abundance for NH3 $\sim$ 360 ppm [Li et al. (2017)]
from the deep atmosphere till the NH3 vapor saturation point. Above the
saturation point, the mixing ratio follows the NH3 saturation vapor pressure
curve. _q ${}_{NH_{3},MWR}$(P)_is the mixing ratio retrieved from MWR
inversion. We use a scaling factor to vary the ammonia mixing ratio between
the homogeneous case to MWR derived profiles. The scaling factor, s ranges
from 0 to 1.5 where 0 is the case for homogeneous mixing ratio. Increasing s
to 1 will change the ammonia profile to MWR inversion case for equator and
mid-latitude regions. We also extend the scaling factor to 1.5, in order to
take into account the low ammonia mixing ratio observed at the North
Equatorial Belt (NEB) of Jupiter [Li et al. (2017)].

NH3 opacity measurements are currently not available for high temperatures (
550 K-3000 K) corresponding to Jupiter’s deep atmosphere and there is a
decrease in the magnitude of absorption of NH3 at high pressures. Thereby, we
invoke a scaling factor to the NH3 absorption coefficient to provide an
estimation of the opacity at high temperatures. The mass absorption
coefficient of ammonia is estimated by multiplying the temperature-scaling law
to the absorption coefficient based on Hanley et al. (2009) (Eq. 7). In this
equation, $\alpha$ is the absorption coefficient of NH3, h is the opacity
factor, _T_ is temperature and _T c_ is reference temperature equal to 750 K.
The NH3 opacity models show that the absorption coefficient peaks at 750 K and
decreases at temperatures beyond 750 K. In the simulations, the scaling factor
is multiplied to the NH3 opacity at temperatures higher than _T c_. The power
law index (h) is varied from 1 to 5 keeping the ammonia concentration
constant, i.e., 2.7 times solar abundance. We also keep the water vapor
constant at 5 times solar abundance as the laboratory measurements demonstrate
that water vapor absorption does not show a significant increase with pressure
and can be said to be relatively transparent when compared to the previous
model of microwave absorption [Steffes et al. (2023)].

$\alpha(NH_{3})\sim\left(\frac{T_{c}}{T}\right)^{h}$ (7)

Changing the ammonia profile and introducing the additional temperature-
dependent scaling factor produce brightness temperature and limb darkening
divergent from MWR data at 0.6 GHz as shown in Figure 4a. The difference
between _T b_ from the model and observations is in the range of 50-200 K at
0.6 GHz. Reducing the ammonia concentration causes a monotonic increase in _T
b_ and a decrease in _L d_. Further, reducing the ammonia opacity shows a
similar trend in _T b_, while a saturation in _L d_ is expected at a power law
factor of 5. Changing the ammonia profile and ammonia opacity has a similar
effect on _T b_ and _L d_ at 1.2 GHz. However, overall, the variation in the
MWR observations at 1.2 GHz can be explained by these two factors and does not
require the inclusion of alkali metals. The 1.2 GHz observations correspond to
$\sim$ 20 bar (Fig. 1), much above the cloud base of alkali metals and at
relatively lower pressure levels. Therefore, the contribution of free
electrons to opacity is expected to be less due to lower temperatures, and the
opacity contribution of ammonia vapor dominates at 1.2 GHz. However, a
comparison of MWR observations at both frequencies clearly implies that the
variation in ammonia vapor opacity cannot solely explain the anomalous
observations at the 0.6 GHz channel.//

Figure 4: Comparison is drawn between the Juno MWR observations and the
results of the radiative transfer model for _T b_ and _L d_ at 0.6 GHz and 1.2
GHz, keeping the water abundance constant $\sim$ 5 times solar abundance. (a,
b) Jupiter’s atmosphere in the absence of alkali metals with only variations
in the NH3 vapor profile and the NH3 opacity, (c, d) Jupiter’s atmosphere in
the presence of alkali metals with variations in the NH3 vapor profile and the
NH3 opacity. The NH3 profile of Jupiter’s atmosphere is varied using a scale
from 0 to 1.5 to take into account the contribution of non-uniform
distribution of NH3 vapor observed by MWR [Li et al. (2017)]. NH3 opacity at
temperatures above 750 K undergoes power law scaling as a function of
atmospheric temperature (Eq. 7). In the absence of alkali metals, the changes
in NH3 vapor profile and the scaling in NH3 vapor opacity deviate
significantly from Juno MWR observations at 0.6 GHz. However, in the presence
of alkali metals of low metallicity, i.e., in the range of -2 to -5, there is
an agreement between model results and MWR observations. Observations at 1.2
GHz can be explained by variations in the NH3 vapor profile and the NH3
opacity independent of opacity contributions from alkali metals.

Fig. 4c, 4d examines the overall effect of alkali metals and ammonia vapor on
the _T b_ and _L d_ at 0.6 GHz and 1.2 GHz. We vary the alkali metal
metallicities in a range from 0 to -7 (solar abundance of Na and K according
to Asplund et al. (2009)) for each condition of NH3 profile scaling and NH3
opacity scaling. The volume mixing ratio of Na and K corresponding to
abundance in solar photosphere [Asplund et al. (2009)] are 3.46 x 10-6
(_[Na/H]_ = -5.76) and 2.14 x 10-7 (_[K/H]_ = -6.97), respectively. Therefore,
we simulate a wide range of ammonia opacity conditions for a given alkali
metal abundance (colored dots). Both NH3 profile and opacity scaling cause a
change in _T b_ and _L d_, which is shown by the annotation in the figure. The
variation in _T b_ and _L d_ is similar to the pattern in Fig. 4a. NH3 profile
scaling causes a decrease in _L d_, while the scaling in NH3 vapor opacity
causes _L d_ to increase at 0.6 GHz. For each case of metallicity, we then
perform a scaling in ammonia vapor and ammonia opacity as described previously
in this section. This provides us with a matrix of _T b_ and _L d_ to take
into account all possible sources of opacity, i.e., collision-induced
absorption, ammonia, water vapor, and free electrons from alkali metals. The
free electron opacity is calculated from the Hartree-Appleton equation
explained in the previous section.

Figure 5: Latitudinal variation of brightness temperature and limb darkening
of Jupiter’s atmosphere at 0.6 GHz and 1.2 GHz at [M/H] = -3

When we compare the new model result with MWR observations (Fig. 4b), we
observe that model matches with observations at 0.6 GHz for free electrons
corresponding to alkali metal metallicities in the range of -2 to -5
(chocolate colored patches), i.e. 10-2 to 10-5 times the solar abundance.
There is an agreement between the model and observations at 1.2 GHz for the
same range of metallicities. The addition of free electrons from alkali metals
dominates the effect of gravity (Fig. 5) and we expect the limb darkening to
decrease from equator to the poles assuming uniform mixing ratio of water and
ammonia vapor. It serves as a baseline to understand the sole effect of free
electrons on latitudinal variation of microwave radiation from Jupiter’s deep
atmosphere.

##  4 Discussions

We infer metallicity of the alkali metals in Jupiter to be much lower than the
solar value. A possible indication of low metallicity of the alkali metals in
a hot Jupiter exoplanet was first proposed by Demory et al. (2011) as one
plausible explanation for the high albedo of Kepler-7b. They derived an alkali
metal abundance 10-100 times lower than the solar value. Since then the
abundance of alkali metals has been derived for several other giant
exoplanets, with abundances ranging from $\sim$ 100 times below solar to
$\sim$ 100 times above solar, although the uncertainties are large. Recent
observations of two hot Jupiters or Saturns with clear or mostly clear
atmospheres were made. The alkali metal abundance for one such hot Jupiter
(HAT-P-1b) is found to be sub-solar [Chen et al. (2022)], while it was found
to be solar to greatly super-solar for the other (WASP-96b)[Nikolov et al.
(2022)]. Considering the relatively small sample size of hot Jupiters with
clear atmospheres, it is premature to make a meaningful comparison between
their alkali metal metallicity and the metallicity in Jupiter presented in
this paper. On the other hand, it is instructive to compare the abundance of
alkali metals in Jupiter from this work with the abundance of the other heavy
elements. While the opacity contribution from alkali metals suggest that Na
and K are strongly depleted relative to solar at the level probed by MWR at
0.6 GHz, all other heavy elements are enriched by a factor of approximately
three to five; while nitrogen is highly variable but enriched, and the water
abundance remains uncertain [Atreya et al. (2019), Li et al. (2020), Li et al.
(2017), Mahaffy et al. (2000)]. The comparison to other heavy metal
measurements from the Galileo probe corresponds to much lower pressures i.e.
$<$ 22 bars. The estimation of alkali metal metallicity from MWR implies lower
metallicity at much higher pressures. The results (Fig. 4b) provide an
important constraint on alkali metal abundance at pressures sensitive to 0.6
GHz channel. A [M/H] = -1 for alkali metals provides too much opacity while
too little abundance or absence of alkali metals does not provide sufficient
opacity to match the MWR observations at 0.6 GHz.

The low abundance of alkali metals indicated by MWR observations could be
attributed to any of the following scenarios. (i) Initially enriched alkali
metals, consistent with the other heavy elements in the atmosphere, are
depleted by chemical reactions with other constituents deep in the atmosphere,
resulting in a low abundance of Na and K at $\sim$ 1 kilobar level sufficient
to provide the free electrons to explain the MWR data at 0.6 GHz. Fegley and
Lodders [Fegley Jr & Lodders (1994)] predict, for example, the formation of
gas-phase species of Na and K in the atmosphere i.e. NaCl, NaOH, and KOH.
Should there be chemical mechanisms that could selectively deplete K in the
atmosphere, leaving Na to be the most significant contributor to free
electrons in the deep atmosphere, the metallicity of Na would be expected to
be in the range of 0 to -2 i.e. solar to highly sub- solar abundance (Appendix
B). (ii) Unconventional planet formation processes, whereby Jupiter did not
accrete a solar complement of alkali metals, or that the alkali metals are not
well mixed at greater depths. If the depletion of alkali metals at $\sim$ 1
kbar inferred in this paper is representative of their bulk abundance, it
could be indicative of the depletion of all rock- forming elements, with
significant implications for the formation and evolution of Jupiter. Our
conclusion of depletion is based on the data of the 0.6 GHz channel, whose
weighting function peaks at 1 kilobar level with the inclusion of alkali
metals. Thus, we are confident about the result only at this level. Alkali
metals could well be more abundant deeper in the atmosphere and they could
have been depleted by some as yet unknown mechanism before reaching the 1
kilobar level though the degree of depletion would have to be huge. Barshay &
Lewis (1978) considered one such possibility, where silicates were found to be
a way of sequestration of gas phase alkali metals. However, a later study by
Fegley Jr & Lodders (1994) found it to be an ineffective mechanism. Further
modeling and laboratory studies are needed to cover the full parameter space
of combined thermochemistry of alkali metal and rock cloud forming species
corresponding to the very high temperature and high pressure conditions of the
deep atmosphere of Jupiter, together with any dynamical effects, before
drawing any firm conclusions about depletion of alkali metals in bulk Jupiter
below the level to which the MWR data of this paper are sensitive.

The new constraints on the abundance of alkalis are linked to their low
ionisation potential, and the fact that the electrons that they provide
directly affect opacities at 0.6 and 1.2 GHz (see Eq. 4). But when present,
they are strong absorbers at visible wavelengths (e.g., Burrows et al. (2000))
and therefore directly affect the planetary radiative flux. The low abundances
that we derive imply that a radiative zone may be present in Jupiter [Guillot
et al. (1994), Guillot et al. (2004)]. Interestingly, this could explain at
the same time the relatively low abundance of CO observed in Jupiter’s
atmosphere compared to expectations for a fully convective deep atmosphere
[Cavalié et al. (2023)].

##  5 Software and third party data repository citations

The software for the radiative transfer package will be available at zenodo
archive (https://doi.org/10.5281/zenodo.7893914) and the MWR data used in this
work, and associated files for data visualization are available at archive (
https://doi.org/10.5281/zenodo.7893817). They can be made available upon
request

##  Appendix A Appendix A: Electron Density and Conductivity

The electron density of Jupiter’s atmosphere is governed by two fundamental
processes: (i) condensation of alkali metal condensates i.e. Na2S and KCl, and
(ii) ionization of alkali metals in thermal equilibrium. Fig. 2 shows the
pressure levels corresponding to the cloud base of Na2S and KCl based on their
saturation vapor pressures. Cloud condensation reduces the amount of alkali
metals available in gas phase that act as a source of free electrons, and
restricts the abundance of Na and K corresponding to their respective
saturation vapor pressure. In the cloud region, electron density is controlled
by saturation vapor pressure of alkali metals whereas below the cloud base,
electron densities are governed by metallicity of alkali metals. Thereby, it
is evident that condensation controls the electron density and thereby,
conductivity at low pressure levels. Condensation limited ionization is
observed at low pressure (below 1 kbar) irrespective of the alkali metal
abundance as the electron density lines converge (Fig. A.1 (a)). Fig. A.1 (a)
and (b) show the presence of a kink in electron density and their respective
conductivity at the cloud base corresponding to different alkali metal
abundances. However, condensation does not play a significant role in
governing the electron densities at $\sim$ 1 kbar pressure level corresponding
to the global maxima in the weighting function at 0.6 GHz (Figure 1).

(a)

(b)

Figure A.1: (a) Electron density of Jupiter’s deep atmosphere at the solar
abundance and _[M/H]_ = -3 and -4, (b) electrical conductivity of Jupiter’s
deep atmosphere at the solar abundance and [M/H] = -3 and -4.

The electron density of the deep atmosphere is much lower than in the case of
alkali metals at solar abundance. It is the true representation of the
electron density of the deep atmosphere. At greater pressures, hydrogen
behaves as a semiconductor and becomes the major contributor to the electron
density [Liu et al. (2008)]. The electrical conductivity of the atmosphere is
calculated using Drude’s equation. It provides an estimate of the conductivity
due to the free electrons provided by alkali metal ionization.

##  Appendix B Appendix B: Selective Depletion of Alkali Metals

Even though Na2S has a deeper condensation level compared to KCl, the cloud
condensation is governed by atmospheric temperature, and does not reflect the
chemical reactivity of alkali metals. K is more electropositive than Na and
thereby, is expected to be more reactive as compared to Na. Therefore, it is
possible that there can be a chemical mechanism that could selectively deplete
K into other compounds, leaving Na as the only source of free electrons in
Jupiter. Under such conditions, we find that Na metallicity should be in the
range of 0 to -3 to match the MWR observations. The increase in alkali metal
metallicities can be attributed to two factors: (i) low ionization energy of
K, and (ii) Na2S condenses much below KCl (Figure 2). Thereby, a larger amount
of Na is required to produce enough free electrons to match the MWR brightness
temperatures and limb darkening.

Figure B.1: Limb darkening and brightness temperature comparison of MWR
observations and forward model results at 600 MHz and 1.2 GHz for
metallicities ranging from 0 to -7 at different ammonia vapor concentration
profiles and opacities. It showcases the sole effect of free electrons due to
the ionization of Na, without considering any contribution from K.

The elimination of K from the atmosphere highlights the role of the elemental
abundance of Na required to match the MWR observations. The results of the
forward model in Fig. B.1 indicate the possible solutions of Na metallicity
under different conditions of ammonia vapor concentration profiles and
microwave opacities. It is observed that the range of Na metallicity is
expected to be from 0 to -3 i.e. solar abundance to highly sub-solar
abundance. Thus, metallicity of Na required is expected to be higher than
those considering both Na and K to be sources of free electrons.

##  Appendix C Appendix C: Jovian Adiabats and Comparison of MWR with High
Temperature Adiabat

Fig. 2 shows that brightness temperatures at 600 MHz from two adiabats differ
by approximately 15 K. The relative weighting function for the adiabats is
that of the ideal Jupiter’s atmosphere without the inclusion of opacity due to
free electrons from alkali metals. It shows a peak at $\sim$ 100 bar. From the
difference in physical temperature of the atmosphere of the two adiabats, it
is seen that the difference reaches $\sim$ 10-15 K at 100 bar level (Fig.
C.1). The weighting function at 600 MHz also extends below 100 bar which could
explain the difference in brightness temperatures. An interesting observation
is that the difference in adiabat temperatures increases with increase in
atmospheric pressure. This increase can be attributed to the temperature
dependent specific heat of the atmospheric constituents.

Figure C.1: Pressure v. temperature difference in temperatures of Jovian
adiabats constructed using T1bar = 168.8 K and T1bar = 166.1 K

The interior models of Jupiter generally use a high temperature in the range
of 170-180 K at the outer boundary (1 bar pressure level) [Gupta et al.
(2022), Miguel et al. (2022)]. These temperatures are about 10-15 K higher
than the measurements from the Galileo probe (166.1 K)[Seiff et al. (1998)]
and Voyager radio occultation reanalysis (168.8 K)[Gupta et al. (2022)]. A
simulation of brightness temperatures and limb darkening at 0.6 GHz and 1.2
GHz is carried out for all cases of alkali metal metallicities, ammonia
concentration and opacity variation assuming _T 1bar_ = 175 K. It can be
clearly seen in Fig. C.2 that high temperature at 1 bar doesn’t match with
entire range of MWR observations for both the frequencies. Some alternate
possibilities could be the presence of a non-adiabatic gradient or a radiative
layer in Jupiter’s deep atmosphere that can possibly account for a higher
temperature at 1 bar level. However, the mismatch with MWR at 1.2 GHz poses a
serious question on the assumption. The current measurements of temperature at
1 bar level are from limited radio occultation experiments. There is a need
for radio science experiments from equator to the poles, in order to estimate
the true variability in temperatures at 1 bar.

Figure C.2: Limb darkening and brightness temperature comparison of MWR
observations and forward model results at 600 MHz and 1.2 GHz for
metallicities ranging from 0 to -7 at different ammonia vapor concentration
profiles and opacities considering T1bar = 175 K

## References

  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual review of astronomy and astrophysics, 47, 481
  * Atreya et al. (1999) Atreya, S., Wong, M., Owen, T., et al. 1999, Planetary and Space Science, 47, 1243
  * Atreya et al. (2019) Atreya, S. K., Crida, A., Guillot, T., et al. 2019, The origin and evolution of Saturn, with exoplanet perspective, in ”Saturn in the 21st Century”, (Eds. Baines K. H. et al.), pp. 5-43, Cambridge University Press Cambridge, UK
  * Bagenal et al. (2007) Bagenal, F., Dowling, T. E., & McKinnon, W. B. 2007, Jupiter: the planet, satellites and magnetosphere, Vol. 1 (Cambridge University Press)
  * Barshay & Lewis (1978) Barshay, S. S., & Lewis, J. S. 1978, Icarus, 33, 593
  * Bellotti et al. (2016) Bellotti, A., Steffes, P. G., & Chinsomboom, G. 2016, Icarus, 280, 255
  * Bhattacharya et al. (2023) Bhattacharya, A., Li, C., Atreya, S. K., et al. 2023, HARP Radiation Package for Planetary Atmospheres, Zenodo, doi: 10.5281/zenodo.7893914
  * Bolton et al. (2017) Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821
  * Burrows et al. (2000) Burrows, A., Marley, M., & Sharp, C. 2000, The Astrophysical Journal, 531, 438
  * Cavalié et al. (2023) Cavalié, T., Lunine, J., & Mousis, O. 2023, Nature Astronomy, 1
  * Chen et al. (2022) Chen, G., Wang, H., van Boekel, R., & Pallé, E. 2022, The Astronomical Journal, 164, 173
  * Chen et al. (2018) Chen, G., Pallé, E., Welbanks, L., et al. 2018, Astronomy & Astrophysics, 616, A145
  * de Pater & Dunn (2003) de Pater, I., & Dunn, D. E. 2003, Icarus, 163, 449
  * de Pater et al. (2019) de Pater, I., Sault, R., Wong, M. H., et al. 2019, Icarus, 322, 168
  * Demory et al. (2011) Demory, B.-O., Seager, S., Madhusudhan, N., et al. 2011, The Astrophysical Journal Letters, 735, L12
  * Devaraj et al. (2014) Devaraj, K., Steffes, P. G., & Duong, D. 2014, Icarus, 241, 165
  * Fegley Jr & Lodders (1994) Fegley Jr, B., & Lodders, K. 1994, Icarus, 110, 117
  * Guillot et al. (1994) Guillot, T., Chabrier, G., Morel, P., & Gautier, D. 1994, Icarus, 112, 354
  * Guillot et al. (2004) Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, The interior of Jupiter, in ”Jupiter: The planet, satellites and magnetosphere”, (Eds. Bagenal F. et al.), Cambridge University Press New York
  * Gupta et al. (2022) Gupta, P., Atreya, S. K., Steffes, P. G., et al. 2022, The Planetary Science Journal, 3, 159
  * Hanley et al. (2009) Hanley, T. R., Steffes, P. G., & Karpowicz, B. M. 2009, Icarus, 202, 316
  * Heays et al. (2017) Heays, A., Bosman, AD, v., & Van Dishoeck, E. 2017, Astronomy & Astrophysics, 602, A105
  * Helliwell (2014) Helliwell, R. A. 2014, Whistlers and related ionospheric phenomena (Courier Corporation)
  * Ingersoll et al. (2017) Ingersoll, A. P., Adumitroaie, V., Allison, M. D., et al. 2017, Geophysical Research Letters, 44, 7676
  * Janssen et al. (2017) Janssen, M., Oswald, J., Brown, S., et al. 2017, Space Science Reviews, 213, 139
  * Karpowicz & Steffes (2011) Karpowicz, B. M., & Steffes, P. G. 2011, Icarus, 212, 210
  * Levin et al. (2001) Levin, S. M., Bolton, S. J., Gulkis, S. L., et al. 2001, Geophysical research letters, 28, 903
  * Li et al. (2018a) Li, C., Ingersoll, A. P., & Oyafuso, F. 2018a, Journal of the Atmospheric Sciences, 75, 1063
  * Li et al. (2018b) Li, C., Le, T., Zhang, X., & Yung, Y. L. 2018b, Journal of Quantitative Spectroscopy and Radiative Transfer, 217, 353
  * Li et al. (2017) Li, C., Ingersoll, A., Janssen, M., et al. 2017, Geophysical Research Letters, 44, 5317
  * Li et al. (2020) Li, C., Ingersoll, A., Bolton, S., et al. 2020, Nature Astronomy, 4, 609
  * Liu et al. (2008) Liu, J., Goldreich, P. M., & Stevenson, D. J. 2008, Icarus, 196, 653, doi: https://doi.org/10.1016/j.icarus.2007.11.036
  * Mahaffy et al. (2000) Mahaffy, P., Niemann, H., Alpert, A., et al. 2000, Journal of Geophysical Research: Planets, 105, 15061
  * Miguel et al. (2022) Miguel, Y., Bazot, M., Guillot, T., et al. 2022, Astronomy and Astrophysics-A&A, 662
  * Morley et al. (2012) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, The Astrophysical Journal, 756, 172
  * Nikolov et al. (2022) Nikolov, N. K., Sing, D. K., Spake, J. J., et al. 2022, Monthly Notices of the Royal Astronomical Society, 515, 3037
  * Saha (1920) Saha, M. N. 1920, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 40, 472
  * Santos-Costa et al. (2017) Santos-Costa, D., Adumitroaie, V., Ingersoll, A., et al. 2017, Geophysical research letters, 44, 8676
  * Seiff et al. (1998) Seiff, A., Kirk, D. B., Knight, T. C., et al. 1998, Journal of Geophysical Research: Planets, 103, 22857
  * Steffes et al. (2023) Steffes, P. G., Bellotti, A., Cembellin, P., Levin, S. M., & Bolton, S. J. 2023, The Planetary Science Journal
  * Steffes et al. (2017) Steffes, P. G., Hanley, T. R., Karpowicz, B. M., et al. 2017, Space Science Reviews, 213, 187
  * Visscher et al. (2006) Visscher, C., Lodders, K., & Fegley Jr, B. 2006, The Astrophysical Journal, 648, 1181
  * Weidenschilling & Lewis (1973) Weidenschilling, S., & Lewis, J. 1973, Icarus, 20, 465
  * Welbanks et al. (2019) Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, The Astrophysical Journal, 887, L20
  * Wong et al. (2004) Wong, M. H., Mahaffy, P. R., Atreya, S. K., Niemann, H. B., & Owen, T. C. 2004, Icarus, 171, 153


# survivalContour: Visualizing predicted survival via colored contour plots

Yushu Shi,
Department of Population Health Sciences, Weill Cornell Medicine
Liangliang Zhang,
Department of Population and Quantitative Health Sciences,
Case Western Reserve University
Kim-Anh Do,
Department of Biostatistics,
The University of Texas MD Anderson Cancer Center
Robert R. Jenq,
Department of Genomic Medicine,
The University of Texas MD Anderson Cancer Center
and
Christine B. Peterson
Department of Biostatistics,
The University of Texas MD Anderson Cancer Center

###### Abstract

Advances in survival analysis have facilitated unprecedented flexibility in
data modeling, yet there remains a lack of tools for graphically illustrating
the influence of continuous covariates on predicted survival outcomes. We
propose the utilization of a colored contour plot to depict the predicted
survival probabilities over time, and provide a Shiny app and R package as
implementations of this tool. Our approach is capable of supporting
conventional models, including the Cox and Fine-Gray models. However, its
capability shines when coupled with cutting-edge machine learning models such
as random survival forests and deep neural networks.

Running title: Survival contour plots

Additional information: Christine B. Peterson, The University of Texas MD
Anderson Cancer Center, 7007 Bertner Ave, Houston, TX, 77030.
<EMAIL_ADDRESS>
The authors declare no potential conflicts of interest.

##  1 Introduction

Over the past few decades, there have been major advances in survival
modeling, including improvements in approaches for handling competing risks
and the development of machine learning methods for survival prediction.
However, there is a lack of appropriate visualization tools for translating
these cutting-edge models into effective graphical representations that can
illustrate the relationships between continuous predictors and survival
probabilities in an intuitive manner. Even for classical models, such as the
Cox or Fine-Gray methods, there is a gap between the model outputs and the
presentation of results.

A typical visualization approach in current practice is to split the covariate
of interest at its median, and present Kaplan-Meier curves for subjects with
low vs. high values. Figures 1(a) and 1(b) display two examples from the
biomedical literature where Kaplan-Meier curves were obtained from binarizing
a continuous predictor at the median value. This approach has significant
shortcomings, as precise information from the continuous predictor gets lost
during binarization (Altman and Royston, 2006; Fernandes et al., 2019;
MacCallum et al., 2002). In addition, the $p$-value for the coefficient in the
Cox model using the continuous version of the predictor does not necessarily
match the $p$-value that one would obtain from a log-rank test comparing the
two Kaplan-Meier curves. This mismatch creates a divergence between the
analysis and the visual presentation. This inconsistency is clear in Figure
1(b), where a hazard ratio derived from a Cox model is displayed alongside a
two-group representation of patient survival.

To highlight interesting trends for extremes of the predictor distribution,
narrower categories may be constructed using quantiles of the predictor.
However, this step introduces an additional element of subjective data
exploration in choosing how to represent the association between the
continuous predictor and survival. Figures 1(c) and 1(d) depict examples from
the literature where tail quantiles, such as the bottom 10% or top 10%, were
identified to highlight differences in the survival curves for patients with
extreme predictor values.

Binary representations with continuous predictor split at the median

(a) HRD score (Abkevich et al., 2012)

(b) Microbiome diversity (Peled et al., 2020)

Multi-group representations with split at extreme quantiles

(c) Metabolomic state (Buergel et al., 2022)

(d) TMB (Samstein et al., 2019)

Spaghetti plots with predicted survival curves for individual subjects

(e) Random survival forests predictions colored by group (Kim et al., 2019)

(f) Deep learning predictions with extreme patients highlighted (Vale-Silva
and Rohr, 2021)

Figure 1: Examples of existing visualization methods

The lack of flexible and informative visualization approaches becomes even
more prominent when considering advanced machine learning-based methods, which
relax the assumption of proportional hazards or accelerated failure times.
Foundational work on the application of machine learning methods to survival
prediction includes the random survival forests (Ishwaran et al., 2008). More
recently, there has been a flurry of proposals on the use of deep learning for
survival prediction, including DeepHit (Lee et al., 2018), DeepSurv (Katzman
et al., 2018), Nnet-survival (Gensheimer and Narasimhan, 2019), and pycox
(Kvamme et al., 2019; Kvamme and Borgan, 2021). The challenge of presenting
survival predictions generated by these approaches is particularly pressing
given the complex non-linear relationships they encode.

To visualize model outputs, spaghetti plots can be used to display predicted
survival or cumulative incidence curves for individual patients. In spaghetti
plots, which have traditionally been applied to depict trends in longitudinal
data, a distinct trajectory, or “noodle”, is plotted for each subject (Szabo,
2009). Figures 1(e) and 1(f) highlight the use of spaghetti plots to display
predicted survival from random survival forests and deep learning models.
Although the curve for an individual subject may be informative, when applied
to a patient cohort spaghetti plots tend to result in a tangle of overlapping
curves without shedding light on the link between covariate values and
predicted survival.

Alternative approaches to visualization have been proposed in the statistical
literature. In particular, the idea of visualizing survival outcomes through a
contour plot appeared in Lumley and Heagerty (2000). However, their method was
implemented in XLisp-Stat, which is not widely used today, and does not
reflect developments in the field of survival analysis over the past two
decades. Currently available tools for survival visualization include the
survminer package (Kassambara et al., 2021), which can produce publication-
ready plots and risk tables, but still focuses on the display of Kaplan-Meier
and cumulative incidence curves. Most recently, Denz and Timmesfeld (2023)
proposed the use of a survival area plot for depicting the effect of a
predictor within the causal inference framework. However, there remains a lack
of user-friendly visualization tools for presenting the results of popular
survival models, especially machine learning based models. To address this
gap, we have designed a Shiny app for interactive visualization and also
provide an R package for more advanced users. Both tools enable the production
of visually effective colored contour plots and offer flexible options for
handling competing risks and interval censoring. Our software captures a broad
set of models for survival association and offers more advanced users seamless
integration to deep learning models.

Figure 2: A snapshot of the survivalContour Shiny app

##  2 Materials and Methods

###  2.1 Shiny app

We have designed a user-friendly Shiny app for automatically generating
survival contour plots, available at
https://biostatistics.mdanderson.org/shinyapps/survivalContour/. A snapshot of
the survivalContour Shiny app is shown above in Figure 2. In the Shiny app, we
guide users to select an appropriate survival model for their data. There are
options for settings with competing risks, where the event of interest may be
precluded by a different event, and for interval-censored data, where the
event time is only known to occur within a given window. Competing risks and
interval censoring are both common in biomedical studies. For example,
competing risks occur when a patient is no longer at risk for cancer
progression because they have died due to an infection. Interval censoring
occurs whenever there are fixed time points for screening or follow-up. For
example, in a cancer screening study, if a scan reveals that the patient
developed a tumor since their last visit, the true time-to-event lies
somewhere during the intervening time window but is not exactly observed. For
settings without competing risks, we offer parametric, nonparametric, semi-
parametric, and machine learning based survival models, as these may be
preferred for different applications. Parametric models, which assume a known
distribution for the survival times, offer an advantage for settings where
out-of-sample prediction is a priority (Jackson, 2016). Nonparametric survival
models, such as spline-based methods, offer greater flexibility on the shape
of the hazard (Royston and Parmar, 2002). Semi-parametric methods, in
particular the Cox model, which avoids assumptions on the baseline hazard
function, are a mainstay in applied survival analysis. We offer both the
classic and stratified Cox model, which allows the baseline hazard function to
differ across levels of a covariate. Finally, machine learning-based methods
offer the most flexible modeling of survival data. Although they do not allow
for inference on the predictors, those methods can provide predicted survival
probabilities for given covariates at specific times. These options are
illustrated as a decision tree in Figure 3.

Figure 3: Decision tree implemented in the Shiny app

Once the desired model has been selected, we ask users to upload their data,
designate the time and censoring indicator, and select the continuous
predictor for visualization. We then produce a 2D contour plot: in this plot,
the x-axis is time, the y-axis is the continuous predictor of interest, and
the intensity of regions in the plot corresponds to the predicted survival
probability from the fitted model. A histogram showing the distribution of the
continuous predictor is displayed alongside the contour plot. In addition to
the major continuous covariate of interest for visualization, users can
provide other covariates to be included as adjusters in the predictive model.
The survival predictions presented in the contour plots are conditional on the
covariate values: these may be specified to a desired level by the user, or
taken as the default of the median value for continuous covariates and the
most frequent value for categorical covariates.

###  2.2 R package

In addition to the Shiny app, we provide an R package for more advanced users,
which allows the display of predictions from deep neural network-based models.
Because of the manual tuning involved with deep neural network-based models,
we restrict the choice of these models to the R package only. The
survivalContour R package can be downloaded from the GitHub page
https://github.com/YushuShi/survivalContour/. We list the models that are
compatible with our tools in Table 1.

Model |  Package
---|---
Cox model and stratified Cox model (Cox, 1972; Therneau and Grambsch, 2000) |  survival (Therneau, 2023)
Cox model and stratified Cox model for interval-censored data |  mets (Scheike et al., 2014; Holst et al., 2016)
Parametric and spline models (Royston and Parmar, 2002) |  flexsurv (Jackson, 2016)
Fine-Gray model for competing risks data (Fine and Gray, 1999) |  riskRegression (Gerds and Kattan, 2021; Gerds et al., 2023)
Fine-Gray model with interval censoring (Bakoyannis et al., 2017) |  intccr (Park et al., 2019)
random survival forests (Ishwaran et al., 2008) |  randomForestSRC (Ishwaran and Kogalur, 2023)
Cox-Time (Kvamme et al., 2019) |  survivalmodels (Sonabend, 2022)
DeepSurv (Katzman et al., 2018) |
DeepHit (Lee et al., 2018) |
Nnet-Survival (Gensheimer and Narasimhan, 2019) |
PC-Hazard (Kvamme and Borgan, 2021) |
Table 1: Models and packages supported by survivalContour

###  2.3 3D plots

Both the Shiny app and the R package support visualization of the survival
contours in three dimensions. This provides a more comprehensive understanding
of predicted survival, since a two-dimensional (2D) colored plot essentially
represents a top-down view of the three-dimensional (3D) structure.
Additionally, the three dimensional representation allows for the illustration
of confidence intervals as semi-transparent layers.

##  3 Results

In this section, we illustrate the capacity of our Shiny app and R package to
present predictions from standard Cox models, stratified Cox models, competing
risks regression, random survival forests, and deep learning methods.

###  3.1 Standard Cox model

To highlight the utility of the Shiny app for presenting the results from a
standard Cox model, we reanalyzed the data presented in Figure 1(b), which was
originally published in Peled et al. (2020). This study established that lower
microbiome diversity at baseline was associated with increased risk of
mortality following allogeneic hematopoietic-cell transplant across multiple
independent patient cohorts. Here, we obtained the data from Cohort 1, which
was composed of patients treated at Memorial Sloan Kettering Cancer Center.
Microbiome diversity was computed using the inverse Simpson index. In the
presentation of the data shown in Figure 1(b), patients were separated into
low vs. high diversity groups based on the median value of diversity for
samples collected between Days 7 and 21. The Kaplan-Meier curves for each
group were presented alongside the hazard ratio for death from a Cox model
with the log10 of diversity as the primary predictor of interest.

A richer representation of this association is presented in Figure 4(a), which
depicts a contour plot of the predicted survival from the Cox model. In this
plot, the x-axis represents time in months after Day 21, the y-axis represents
microbiome diversity on the log10 scale, and the color intensity represents
the predicted survival probability. As an example of how to interpret the
plot, at 10 months, a subject with microbiome diversity $<0.2$ on the log10
scale would have a predicted survival probability in the range 70-75%, while a
subject with a diversity value of 1.2 on the log10 scale would have a
predicted survival of 80-85%. This plot corresponds more closely to the
analysis performed than the Kaplan-Meier curves and offers an intuitive
representation of the strength of association between the predictor and the
outcome. The orange histogram at right provides information on the
distribution of the predictor values as well as the number of subjects
included in the analysis (via the bar heights). Here, the bar heights provide
a key insight that very low diversity values (which confer worse expected
survival) are in fact quite common in the observed data for this cohort, as
reflected by the histogram mode occurring at the lowest bin.

In Figure 4(b), we show the predicted survival curves for 5 quantiles of the
predictor values. These plots are automatically generated by the app, and
provide a snapshot of the predicted survival at evenly spaced quantile values.
Importantly, these plots include confidence intervals to highlight the degree
of statistical uncertainty in the predictions.

Here, we presented the results from a univariate Cox model for the sake of
simplicity. However, the survivalContour app allows for the inclusion of
covariates in the model.

(a) Contour plot. Diversity corresponds to microbiome diversity on the log10
scale.

(b) Predicted survival curves. Diversity corresponds to baseline microbiome
diversity on the log10 scale.

Figure 4: Plots created by survivalContour Shiny app using Peled et al.
(2020)’s data

###  3.2 Stratified Cox model

To demonstrate the capabilities of the survivalContour Shiny app in tandem
with the stratified Cox model, we provide an illustrative example using the
Veterans’ Administration Lung Cancer study (Kalbfleisch and Prentice, 2011),
available in the survival R package. The data set includes covariate data and
survival outcomes for 137 subjects. Our focus is on elucidating the influence
of the Karnofsky performance score (Karnofsky, 1949) on survival times. Since
this relationship varies by cancer subtype, we fit a stratified Cox model,
with strata defined by the four subtypes in the dataset: adenocarcinoma, large
cell carcinoma, small cell lung cancer, and squamous cell carcinoma. We
adjusted for the additional covariates age, prior therapy status, treatment
arm, and months from diagnosis to randomization. The predicted survival
probabilities are shown in Figure 5, where the projected survival for each
cancer type is presented in a subpanel. Higher Karnofsky scores are associated
with improved survival predictions across all four subtypes: however, the
strength of this association and overall hazard rates vary by cancer subtype.

Figure 5: Plot generated by the survivalContour R package for stratified
analysis. Each subpanel corresponds to a lung cancer subtype: adenocarcinoma
(orange), large cell carcinoma (red), small cell carcinoma (brown), and
squamous cell carcinoma (gray).

###  3.3 Competing risks

The survivalContour Shiny app, along with its accompanying R package, can
handle data with competing risks. To demonstrate this capacity, we use the
Paquid dataset, which was collected as part of a study on brain ageing
(Dartigues et al., 1992) and is included in the riskRegression package (Gerds
et al., 2023). The data set includes a total of 2561 participants, with
observed event times for the onset of dementia for 449 subjects and for death
without dementia for 634 subjects. We focused on predicting the onset of
dementia while considering death without dementia as a competing risk.

As our primary predictor variable of interest, we considered the Digit Symbol
Substitution Score Test (DSST), which explores attention and psychomotor
speed. Concurrently, we adjust for the Mini-Mental State Examination score
(MMSE), which is a measure of overall cognitive performance. Employing the
Fine and Gray model, we illustrate the influence of DSST in Figure 6 via both
2D and 3D contour representations. Here, lower scores on the DSST are
associated with an increased probability of developing dementia over time,
accounting for the fact that subjects that died without dementia were removed
from the risk pool.

Figure 6: 2D plot (left) and 3D plot (right) generated by the survivalContour
R package for the Paquid competing risks data.

###  3.4 Machine learning methods

The proposed contour plot visualization shines when integrated with advanced
machine learning models like random survival forests and deep neural networks.
In this section, we rely on data collected by the SUPPORT III study (Knaus et
al., 1995), which aimed to characterize survival outcomes of seriously ill
hospitalized adults. This dataset, which contains 8873 observations, was also
analyzed in the paper introducing the DeepSurv method (Katzman et al., 2018).
In this case study, we are interested in depicting the influence of
respiratory rate on survival. Additional covariates in the data set include
age, sex, race, the number of comorbidities, the presence of cancer, diabetes,
and dementia, mean arterial blood pressure, temperature, white blood count,
heart rate, and serum creatine.

Figure 7(a) depicts 2D and 3D survival contour plots showing the relationship
between respiratory rate and predicted survival from a random survival forest
model. The color variation in the contour plot corresponds to different
survival probabilities, with darker colors representing worse predicted
survival. The 2D and 3D visualizations from the DeepSurv method are shown in
Figure 7(b). What static Figures 7(a) and 7(b) cannot reveal is that users can
interactively obtain predicted survival rates at specific time points and
covariate values by gliding the mouse across the plot. Also, 3D interactive
graphics allow users to rotate the 3D surface to view them from multiple
angles.

Since the random survival forest and DeepSurv utilize different approaches for
prediction, the resulting contour plots are not exactly the same.
Nevertheless, the non-monotone relation between respiratory rate and survival
is captured by both models: survival rates deteriorate at extreme respiratory
rates, whereas the median respiratory rate forecasts better survival, aligning
with clinical intuition. This delicate relationship would not be captured by a
basic Cox model, which imposes a multiplicative hazards assumption. Moreover,
the common practice of dichotomizing the predictor would fail to identify this
relationship.

To illustrate the limitation of existing methods, Figure 7(c) displays Kaplan-
Meier curves for both the high and the low respiratory rate groups, using a
cutoff set at the median value. These curves appear perfectly smooth due to
the large number of subjects in the patient cohort ($>$ 8000). Additionally,
we randomly selected 20 samples from each group and plotted their predicted
survival from the DeepSurv model. The non-linear association evident in Figure
7(a) and 7(b) becomes lost in these presentations, which reflect common
practice in presenting survival predictions.

(a) 2D plot (left) and 3D plot (right) generated by the survivalContour R
package for SUPPORT III study using random survival forests model.

(b) 2D plot (left) and 3D plot (right) generated by the survivalContour R
package to depict predictions from a DeepSurv model.

(c) Kaplan Meier curves of the high and the low respiratory rate groups (left)
and survival probabilities of randomly selected samples from each group
(right).

Figure 7: Survival predictions generated using machine learning methods

##  4 Discussion

Our survivalContour software opens a new avenue for visualizing the
relationship between a continuous predictor and survival outcomes. The Shiny
app offers a user-friendly guided experience and can accommodate many popular
survival models. The full potential of our plotting approach is revealed when
incorporated with state-of-the-art machine learning models through the
survivalContour R package.

Importantly, the survival contour plot focuses on displaying predicted
survival and does not aim to provide insight regarding model fit. We recommend
that users compute Harrell’s Concordance Index (C-index, Harrell et al., 1982)
and the integrated Brier score (Graf et al., 1999) to evaluate model fit
before sharing results from survivalContour. Complex deep neural network-based
models, like DeepHit or DeepSurv, often require significant manual tuning and
the best practices for fitting such complex models are beyond the scope of
this paper. Additionally, flexible models may be susceptible to overfitting,
so it is always advisable to use an independent dataset for validation.

### Data availability statement

The data presented in Section 3.1 were originally published by Peled et al.
(2020). Restrictions apply to the availability of these data, which were used
under license for this study. The Veterans’ Administration Lung Cancer data
analyzed in Sections 3.2 is available in the R package survival. The Paquid
data analyzed in Section 3.3 is available in the R package riskRegression. The
data from the SUPPORT III study analyzed in Section 3.4 is available from the
GitHub page https://github.com/havakv/pycox/.

### Acknowledgments

This work was partially funded by National Institutes of Health grants
P30CA016672, SPORE P50CA140388, R01 HL158796, R01 HL124112; CCTS TR000371,
Cancer Prevention and Research Institute of Texas RP160693, CCSG P30CA016672
and National Science Foundation grant 2310955. We would like to thank Jonathan
Peled and Keimya Sadeghi of Memorial Sloan Kettering Cancer Center for sharing
the MSKCC cohort data from Peled et al. (2020) presented in Section 3.1.

## References

  * (1)
  * Abkevich et al. (2012) Abkevich, V., Timms, K. M., Hennessy, B. T., Potter, J., Carey, M. S., Meyer, L. A. et al. (2012), ‘Patterns of genomic loss of heterozygosity predict homologous recombination repair defects in epithelial ovarian cancer’, British Journal of Cancer 107(10), 1776–1782.
  * Altman and Royston (2006) Altman, D. G. and Royston, P. (2006), ‘The cost of dichotomising continuous variables’, BMJ 332(7549), 1080.
  * Bakoyannis et al. (2017) Bakoyannis, G., Yu, M. and Yiannoutsos, C. T. (2017), ‘Semiparametric regression on cumulative incidence function with interval-censored competing risks data’, Statistics in Medicine 36(23), 3683–3707.
  * Buergel et al. (2022) Buergel, T., Steinfeldt, J., Ruyoga, G., Pietzner, M., Bizzarri, D., Vojinovic, D. et al. (2022), ‘Metabolomic profiles predict individual multidisease outcomes’, Nature Medicine 28(11), 2309–2320.
  * Cox (1972) Cox, D. R. (1972), ‘Regression models and life-tables’, Journal of the Royal Statistical Society: Series B (Methodological) 34(2), 187–202.
  * Dartigues et al. (1992) Dartigues, J.-F., Gagnon, M., Barberger-Gateau, P., Letenneur, L., Commenges, D., Sauvel, C., Michel, P. and Salamon, R. (1992), ‘The Paquid epidemiological program on brain ageing’, Neuroepidemiology 11(Suppl. 1), 14–18.
  * Denz and Timmesfeld (2023) Denz, R. and Timmesfeld, N. (2023), ‘Visualizing the (causal) effect of a continuous variable on a time-to-event outcome’, Epidemiology 34(5), 652.
  * Fernandes et al. (2019) Fernandes, A., Malaquias, C., Figueiredo, D., da Rocha, E. and Lins, R. (2019), ‘Why quantitative variables should not be recoded as categorical’, Journal of Applied Mathematics and Physics 7, 1519–1530.
  * Fine and Gray (1999) Fine, J. P. and Gray, R. J. (1999), ‘A proportional hazards model for the subdistribution of a competing risk’, Journal of the American Statistical Association 94(446), 496–509.
  * Gensheimer and Narasimhan (2019) Gensheimer, M. F. and Narasimhan, B. (2019), ‘A scalable discrete-time survival model for neural networks’, PeerJ 7, e6257.
  * Gerds and Kattan (2021) Gerds, T. A. and Kattan, M. W. (2021), Medical Risk Prediction Models: With Ties to Machine Learning (1st ed.), Chapman and Hall/CRC.
  * Gerds et al. (2023) Gerds, T. A., Ohlendorff, J. S. and Ozenne, B. (2023), riskRegression: Risk Regression Models and Prediction Scores for Survival Analysis with Competing Risks.  R package version 2023.03.22.
  * Graf et al. (1999) Graf, E., Schmoor, C., Sauerbrei, W. and Schumacher, M. (1999), ‘Assessment and comparison of prognostic classification schemes for survival data’, Statistics in Medicine 18(17-18), 2529–2545.
  * Harrell et al. (1982) Harrell, F. E. J., Califf, R. M., Pryor, D. B., Lee, K. L. and Rosati, R. A. (1982), ‘Evaluating the yield of medical tests’, Journal of American Medical Association 247(18), 2543–2546.
  * Holst et al. (2016) Holst, K. K., Scheike, T. H. and Hjelmborg, J. B. (2016), ‘The liability threshold model for censored twin data’, Computational Statistics and Data Analysis 93, 324–335.
  * Ishwaran and Kogalur (2023) Ishwaran, H. and Kogalur, U. (2023), Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC).  R package version 3.2.2.
  * Ishwaran et al. (2008) Ishwaran, H., Kogalur, U., Blackstone, E. and Lauer, M. (2008), ‘Random survival forests’, Annals of Applied Statistics 2(3), 841–860.
  * Jackson (2016) Jackson, C. (2016), ‘flexsurv: A platform for parametric survival modeling in R’, Journal of Statistical Software 70(8), 1–33.
  * Kalbfleisch and Prentice (2011) Kalbfleisch, J. D. and Prentice, R. L. (2011), The statistical analysis of failure time data, John Wiley & Sons.
  * Karnofsky (1949) Karnofsky, D. A. (1949), ‘The clinical evaluation of chemotherapeutic agents in cancer’, Evaluation of Chemotherapeutic Agents pp. 191–205.
  * Kassambara et al. (2021) Kassambara, A., Kosinski, M. and Biecek, P. (2021), survminer: Drawing Survival Curves using ’ggplot2’.  R package version 0.4.9.
https://CRAN.R-project.org/package=survminer

  * Katzman et al. (2018) Katzman, J. L., Shaham, U., Cloninger, A., Bates, J., Jiang, T. and Kluger, Y. (2018), ‘DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network’, BMC Medical Research Methodology 18(1), 24.
  * Kim et al. (2019) Kim, D. W., Lee, S., Kwon, S., Nam, W., Cha, I.-H. and Kim, H. J. (2019), ‘Deep learning-based survival prediction of oral cancer patients’, Scientific Reports 9(1), 6994.
  * Knaus et al. (1995) Knaus, W. A., Harrell, F. E., Lynn, J., Goldman, L., Phillips, R. S., Connors, A. F. et al. (1995), ‘The SUPPORT prognostic model: Objective estimates of survival for seriously ill hospitalized adults’, Annals of Internal Medicine 122(3), 191–203.
  * Kvamme and Borgan (2021) Kvamme, H. and Borgan, Ø. (2021), ‘Continuous and discrete-time survival prediction with neural networks’, Lifetime Data Analysis 27(4), 710–736.
  * Kvamme et al. (2019) Kvamme, H., Borgan, Ø. and Scheel, I. (2019), ‘Time-to-event prediction with neural networks and cox regression’, Journal of Machine Learning Research 20, 129:1–129:30.
  * Lee et al. (2018) Lee, C., Zame, W., Yoon, J. and van der Schaar, M. (2018), ‘DeepHit: A deep learning approach to survival analysis with competing risks’, Proceedings of the AAAI Conference on Artificial Intelligence 32(1).
  * Lumley and Heagerty (2000) Lumley, T. and Heagerty, P. (2000), ‘Graphical exploratory analysis of survival data’, Journal of Computational and Graphical Statistics 9(4), 738–749.
  * MacCallum et al. (2002) MacCallum, R. C., Zhang, S., Preacher, K. J. and Rucker, D. D. (2002), ‘On the practice of dichotomization of quantitative variables’, Psychological methods 7(1), 19–40.
  * Park et al. (2019) Park, J., Bakoyannis, G. and Yiannoutsos, C. (2019), ‘Semiparametric competing risks regression under interval censoring using the r package intccr’, Computer Methods and Programs in Biomedicine 173, 167–176.
  * Peled et al. (2020) Peled, J. U., Gomes, A. L., Devlin, S. M., Littmann, E. R., Taur, Y., Sung, A. D. et al. (2020), ‘Microbiota as predictor of mortality in allogeneic hematopoietic-cell transplantation’, New England Journal of Medicine 382(9), 822–834.
  * Royston and Parmar (2002) Royston, P. and Parmar, M. K. (2002), ‘Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects’, Statistics in Medicine 21(15), 2175–2197.
  * Samstein et al. (2019) Samstein, R. M., Lee, C.-H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y. et al. (2019), ‘Tumor mutational load predicts survival after immunotherapy across multiple cancer types’, Nature Genetics 51(2), 202–206.
  * Scheike et al. (2014) Scheike, T. H., Holst, K. K. and Hjelmborg, J. B. (2014), ‘Estimating heritability for cause specific mortality based on twin studies’, Lifetime Data Analysis 20(2), 210–233.
  * Sonabend (2022) Sonabend, R. (2022), survivalmodels: Models for Survival Analysis.  R package version 0.1.13.
  * Szabo (2009) Szabo, A. (2009), ‘Spaghetti is good for more than just eating’, Datum: Newsletter of the Biostatistics & Epidemiology Key Function of the CTSI 15(4).
  * Therneau (2023) Therneau, T. M. (2023), A Package for Survival Analysis in R.  R package version 3.5-5.
https://CRAN.R-project.org/package=survival

  * Therneau and Grambsch (2000) Therneau, T. M. and Grambsch, P. M. (2000), Modeling Survival Data: Extending the Cox Model, Springer, New York.
  * Vale-Silva and Rohr (2021) Vale-Silva, L. A. and Rohr, K. (2021), ‘Long-term cancer survival prediction using multimodal deep learning’, Scientific Reports 11(1), 13505.